{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport cv2 # computer vision\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport matplotlib\nmatplotlib.rcParams['image.cmap'] = 'bone'\nimport pydicom\nfrom glob import glob\nfrom tqdm import tqdm\nfrom sklearn.metrics import precision_score, recall_score\nfrom sklearn.model_selection import train_test_split\nfrom skimage.exposure import equalize_adapthist","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"seed=42\nversion=37\nBASE_PATH = '../input/rsna-intracranial-hemorrhage-detection/'\nTRAIN_DIR = 'rsna-intracranial-hemorrhage-detection/stage_2_train/'\nTEST_DIR = 'rsna-intracranial-hemorrhage-detection/stage_2_test/'\nDICOM_PATH = '../input/dicom-dataframe/df_trn.fth' # Obtained from https://www.kaggle.com/jhoward/creating-a-metadata-dataframe-fastai/output","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"df = pd.read_csv(BASE_PATH + '/rsna-intracranial-hemorrhage-detection/stage_2_train.csv').rename(columns={'Label': 'label'})\ndf[['id', 'img', 'subtype']] = df['ID'].str.split('_', n=3, expand=True)\ndf['filename'] = df['id'] + '_' + df['img'] + '.dcm'\ndf['path'] = BASE_PATH + TRAIN_DIR + df['filename']\n\nbdf = pd.read_feather(DICOM_PATH).rename(columns={'SOPInstanceUID': 'filename'}).assign(filename=lambda x: x['filename']+'.dcm')[['filename', 'img_pct_window']]\ndf = df.merge(bdf, on='filename', how='left').drop(['ID'], axis=1)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['label'].value_counts().plot(kind='bar', grid=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['label'] == 1]['subtype'].value_counts().plot(kind='bar', grid=True, colors='blue')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Set a threshold to the minimum amount of pixels containing brain tissues."},{"metadata":{"trusted":true},"cell_type":"code","source":"brain_threshold = 0.15\ndf = df[(df['img_pct_window'] > brain_threshold)].reset_index(drop=True)\nprint(df.shape)\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Preprocessing the images before viewing, code taken from [here](https://www.kaggle.com/akensert/inceptionv3-prev-resnet50-keras-baseline-model) further reading [here](http://uwmsk.org/jupyter/Jupyter_DICOM_toolbox.html)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def window_correction(dcm, window_center, window_width):\n    if (dcm.BitsStored == 12) and (dcm.PixelRepresentation == 0) and (int(dcm.RescaleIntercept) > -100):\n        dcm = fix_error_img(dcm)\n    img = dcm.pixel_array * dcm.RescaleSlope + dcm.RescaleIntercept\n    img_min = window_center - window_width // 2\n    img_max = window_center + window_width // 2\n    img = np.clip(img, img_min, img_max)\n    return img\n\ndef fix_error_img(dcm):\n    if dcm.PixelRepresentation != 0 or dcm.RescaleIntercept<-100: \n        return dcm\n    x = dcm.pixel_array + 1000\n    px_mode = 4096\n    x[x>=px_mode] = x[x>=px_mode] - px_mode\n    dcm.PixelData = x.tobytes()\n    dcm.RescaleIntercept = -1000\n    return dcm\n\ndef image_preprocessing(img, target_size=(256, 256)):\n    brain_img = window_correction(img, 40, 80)\n    subdural_img = window_correction(img, 80, 200)\n    soft_img = window_correction(img, 40, 380)\n    \n    brain_img = (brain_img - 0) / 80\n    subdural_img = (subdural_img - (-20)) / 200\n    soft_img = (soft_img - (-150)) / 380\n    \n    bsb_img = np.array([brain_img, subdural_img, soft_img]).transpose(1,2,0)\n    bsb_img = cv2.resize(bsb_img, target_size[:2], interpolation=cv2.INTER_LINEAR)\n    return bsb_img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Analyze label 1\nLabel 1 are image that contain no damage, a healthy brain"},{"metadata":{"trusted":true},"cell_type":"code","source":"one_df = df[df['label'] == 1]\nrow=1; col=3;\nfig = plt.figure(figsize=(10, 10)) #ID_39e5d2a39.dcm ID_f14b31779.dcm ID_330d2fd8c ID_ff776bf6b.dcm ID_8e373a29e\nsmp1 = []\nfor i in range(3):\n    sample=np.random.randint(one_df.shape[0])\n    smp1.append(sample)\n    fig.add_subplot(row, col, i+1)\n    data = one_df.iloc[sample]\n    img = image_preprocessing(pydicom.dcmread(data['path']))\n    plt.imshow(img, cmap=plt.cm.bone)\n    plt.title(data['filename'] + '_' + data['label'].astype(str))\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Analyze label 0\nLabel 0 corresponds to images contains no injuries"},{"metadata":{"trusted":true},"cell_type":"code","source":"zero_df = df[df['label'] == 0]\nrow=1; col=3;\nfig = plt.figure(figsize=(10, 10))\nsmp0 = []\nfor i in range(3):\n    sample=np.random.randint(zero_df.shape[0])\n    smp0.append(sample)\n    fig.add_subplot(row, col, i+1)\n    data = zero_df.iloc[sample]\n    img = image_preprocessing(pydicom.dcmread(data['path']))\n    plt.title(data['filename'] + data['label'].astype(str))\n    plt.imshow(img, cmap=plt.cm.bone)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_subtypes(subtype, salt=0, figsize=(10, 5)):\n    sample_df = positive_df[positive_df['subtype'] == subtype].sample(5, random_state=seed+salt)['path'].values\n    \n    fig, ax = plt.subplots(1, 5, figsize=figsize)\n    for idx, val in enumerate(sample_df):\n        img = image_preprocessing(pydicom.dcmread(val))\n        ax[idx].imshow(img)\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"positive_df = df[df['label'] == 1]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Epidural\nEpidural hematoma (EDH) is a traumatic accumulation of blood between the inner table of the skull and the stripped-off dural membrane. EDH results from traumatic head injury, usually with an associated skull fracture and arterial laceration. [Source](https://emedicine.medscape.com/article/824029-overview)\n\nThe images below show the blood at the skull perimeter"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_subtypes('epidural', figsize=(20, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Intraparenchymal\nIntraparenchymal hemorrhage appearance is generally as bright white acutely. The size may vary from punctate to catastrophically large, with associated mass effect and midline shift. For intraparenchmyal hemorrhage, mass effect such as midline shift or ventricular effacement should be assessed. Intraparenchymal hemorrhage is similar to subarachnoid hemorrhage. [Source](https://www.sciencedirect.com/topics/neuroscience/intraparenchymal-hemorrhage)\n\nThe images below (3, 4) show a bright white acutely while the others show the mass in varying spots"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_subtypes('intraparenchymal', figsize=(20, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Intraventricular\nIntraventricular hemorrhage (IVH) is bleeding inside or around the ventricles in the brain. The ventricles are the spaces in the brain that contain the cerebral spinal fluid. Bleeding in the brain can put pressure on the nerve cells and damage them. Severe damage to cells can lead to brain injury. [Source](https://www.urmc.rochester.edu/encyclopedia/content.aspx?contenttypeid=90&contentid=P02608)\n\nThe second image below shows blood in the cerebral fluid but the other images are difficult to tell."},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_subtypes('intraventricular', figsize=(20, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Subarachnoid\nThe subarachnoid space is the area between the brain and the skull. It is filled with cerebrospinal fluid (CSF), which acts as a floating cushion to protect the brain. When blood is released into the subarachnoid space, it irritates the lining of the brain, increases pressure on the brain, and damages brain cells. At the same time, the area of brain that previously received oxygen-rich blood from the affected artery is now deprived of blood, resulting in a stroke. [Source](https://mayfieldclinic.com/pe-sah.htm)\n\nAll the images below show blood in the subarachnoid space. The last image shows "},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_subtypes('subarachnoid', figsize=(20, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Any\nAs the label states these images could be any of the 5 above. I am no doctor but from what I learned above my guess is:\n - Epidural\n - Subarachnoid\n - Intraparenchymal\n - Intraventricular\n - Epidual"},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_subtypes('any', figsize=(20, 10))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"df = df.loc[:, [\"label\", \"subtype\", \"filename\"]].drop_duplicates(['filename', 'subtype'])\ndf = df.set_index(['filename', 'subtype']).unstack(level=-1).droplevel(0, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.utils import class_weight\n\ntdf = df[df.sum(axis=1) != 0]\nfor idx, i in enumerate(df.columns):\n    tdf[i] = tdf[i].map({1: idx})\n    \ntdf = tdf.fillna(0)\ntdf.head()\n\nytdf = tdf.values.flatten('f')\ndel tdf\n\ncw = class_weight.compute_class_weight('balanced', classes=np.unique(ytdf), y=ytdf)\nprint(cw)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(df.sum() / df.shape[0])\nhealthy_images = (df.sum(axis=1) == 0)\nprint(f'\\nHealthy images: {healthy_images.sum()}')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df = pd.read_csv(BASE_PATH + '/rsna-intracranial-hemorrhage-detection/stage_2_sample_submission.csv')\ntest_df[\"Image\"] = test_df[\"ID\"].str.slice(stop=12)\ntest_df[\"Diagnosis\"] = test_df[\"ID\"].str.slice(start=13)\n\ntest_df = test_df.loc[:, [\"Label\", \"Diagnosis\", \"Image\"]]\ntest_df = test_df.set_index(['Image', 'Diagnosis']).unstack(level=-1)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Create model"},{"metadata":{"trusted":true},"cell_type":"code","source":"import keras.applications as ka\nfrom keras.models import Sequential, Model, load_model, Input\nfrom keras.optimizers import Adam\nfrom keras import layers\nfrom keras.layers import Input, Add, Dense, Activation, ZeroPadding2D, BatchNormalization, Flatten, Conv2D, AveragePooling2D, MaxPooling2D, GlobalMaxPooling2D,Dropout, GlobalAveragePooling2D\nfrom keras.models import Model, load_model\nfrom keras.preprocessing import image\nfrom keras.utils import layer_utils, Sequence\nfrom keras.wrappers.scikit_learn import KerasClassifier\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.initializers import glorot_uniform\nimport keras.backend as K\nfrom keras.callbacks import ModelCheckpoint, LearningRateScheduler, TensorBoard, EarlyStopping, ReduceLROnPlateau, Callback\nfrom imgaug import augmenters as iaa","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DENSE = 128\nDROPOUT = 0.4\nLEARNING_RATE = 2e-3\nBATCH_SIZE = 32\nEPOCHS = 6\nN_CLASS = 6\n\ntarget_size=(224, 224, 3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def weighted_categorical_loss(y_true, y_pred):\n    class_weights = np.array([2., 1., 1., 1., 1., 1.])\n    eps = K.epsilon()\n    y_pred = K.clip(y_pred, eps, 1.0-eps)\n\n    out = -(y_true  * K.log(y_pred) * class_weights + (1.0 - y_true) * K.log(1.0 - y_pred) * class_weights)\n    return K.mean(out, axis=-1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class DataGenerator(Sequence):\n\n    def __init__(self, list_IDs, labels=None, batch_size=1, img_size=target_size, img_dir='rsna-intracranial-hemorrhage-detection/stage_2_train/'):\n\n        self.list_IDs = list_IDs\n        self.labels = labels\n        self.batch_size = batch_size\n        self.img_size = target_size\n        self.img_dir= BASE_PATH + img_dir\n        self.on_epoch_end()\n\n    def __len__(self):\n        return int(np.ceil(len(self.indices) / self.batch_size))\n\n    def __getitem__(self, index):\n        indices = self.indices[index*self.batch_size:(index+1)*self.batch_size]\n        list_IDs_temp = [self.list_IDs[k] for k in indices]\n        \n        if self.labels is not None:\n            X, Y = self.__data_generation(list_IDs_temp)\n            return X, Y\n        else:\n            X = self.__data_generation(list_IDs_temp)\n            return X\n        \n    def on_epoch_end(self):\n        if self.labels is not None: # for training phase we undersample and shuffle\n            # keep probability of any=0 and any=1\n            keep_prob = self.labels.iloc[:, 0].map({0: 0.35, 1: 0.5})\n            keep = (keep_prob > np.random.rand(len(keep_prob)))\n            self.indices = np.arange(len(self.list_IDs))[keep]\n            np.random.shuffle(self.indices)\n        else:\n            self.indices = np.arange(len(self.list_IDs))\n\n    def read_dicom(self, path, target_size):\n        dcm = pydicom.dcmread(path)\n\n        try:\n            img = image_preprocessing(dcm, target_size=target_size)\n        except:\n            img = np.zeros(target_size)\n\n        return img\n            \n    def augment_img(self, image): \n        augment_img = iaa.Sequential([\n            iaa.Sometimes(0.25, iaa.GaussianBlur(sigma=(0, 0.5))),\n            iaa.Fliplr(0.25),\n            iaa.Flipud(0.25)])\n        image_aug = augment_img.augment_image(image)\n        return image_aug\n\n    def __data_generation(self, list_IDs_temp):\n        X = np.empty((self.batch_size, *self.img_size))\n        \n        if self.labels is not None: # training phase\n            Y = np.empty((self.batch_size, 6), dtype=np.float32)\n        \n            for i, ID in enumerate(list_IDs_temp):\n                #X[i,] = self.augment_img(self.read_dicom(self.img_dir+ID, self.img_size))\n                X[i,] = self.read_dicom(self.img_dir+ID, self.img_size)\n                Y[i,] = self.labels.loc[ID].values\n        \n            return X, Y\n        \n        else:\n            for i, ID in enumerate(list_IDs_temp):\n                X[i,] = self.read_dicom(self.img_dir+ID+'.dcm', self.img_size)\n            \n            return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"class PredictionCheckpoint(Callback):\n    \n    def __init__(self, test_df, batch_size=BATCH_SIZE, input_size=target_size):\n        \n        self.test_df = test_df\n        self.test_images_dir = BASE_PATH + TEST_DIR\n        self.batch_size = batch_size\n        self.input_size = input_size\n        \n    def on_train_begin(self, logs={}):\n        self.test_predictions = []\n        \n    def on_epoch_end(self,batch, logs={}):\n        self.test_predictions.append(\n            self.model.predict_generator(\n                DataGenerator(self.test_df.index, None, self.batch_size, self.input_size, img_dir=TEST_DIR), verbose=2)[:len(self.test_df)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def build_model(input_shape, pretrained_model=None):\n    net = pretrained_model(include_top=False, input_shape=input_shape, weights='imagenet')\n    \n    model = Sequential()\n    model.add(net)\n    model.add(GlobalAveragePooling2D())\n    #model.add(BatchNormalization())\n    #model.add(Dropout(DROPOUT))\n    \n    #model.add(Dense(DENSE, activation='relu'))\n    #model.add(BatchNormalization())\n    #model.add(Dropout(DROPOUT-0.3))\n    \n    model.add(Dense(N_CLASS, activation='sigmoid'))\n    model.summary()\n    \n    model.compile(loss='binary_crossentropy', optimizer=Adam(), metrics=[weighted_categorical_loss])\n    return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_idx, valid_idx = train_test_split(df.index, test_size=0.01, random_state=42)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def lr_scheduler(epoch, lr):\n    decay_rate = 0.85\n    decay_step = 1\n    if epoch % decay_step == 0 and epoch:\n        return lr * pow(decay_rate, np.floor(epoch / decay_step))\n    return lr","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_generator = DataGenerator(train_idx, labels=df.loc[train_idx], batch_size=BATCH_SIZE, img_size=target_size)\n#valid_generator = DataGenerator(valid_idx, labels=df.loc[valid_idx], batch_size=len(valid_idx), img_size=target_size)\n\nmodel_path = f'InceptionV3_{version}_{EPOCHS}_{DENSE}.h5'\npred_history = PredictionCheckpoint(test_df)\n\ncallbacks = [\n    LearningRateScheduler(lr_scheduler, verbose=1),\n    ModelCheckpoint(model_path, monitor='val_loss', verbose=1, save_best_only=True, save_weights_only=False, mode='auto'),\n    pred_history\n]\n    \nmodel = build_model(pretrained_model=ka.InceptionV3, input_shape=(224, 224, 3))\nhistory = model.fit_generator(train_generator, callbacks=callbacks, \n                              epochs=EPOCHS, verbose=1, use_multiprocessing=True, workers=4, class_weight=cw)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"hdf = pd.DataFrame(history.history)\nhdf[['loss', 'weighted_categorical_loss']].plot(grid=True, figsize=(15, 3), title='Loss and Accuracy Graphs')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.iloc[:, :] = np.average(pred_history.test_predictions, axis=0, weights=[0, 0, 1, 2, 4, 8])\ntest_df = test_df.stack().reset_index()\ntest_df.insert(loc=0, column='ID', value=test_df['Image'].astype(str) + \"_\" + test_df['Diagnosis'])\ntest_df = test_df.drop([\"Image\", \"Diagnosis\"], axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"test_df.to_csv('submission.csv', index=False)","execution_count":null,"outputs":[]}],"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":4,"nbformat_minor":1}