{"cells":[{"metadata":{},"cell_type":"markdown","source":"### Lung Segmentation from RANZCR Chest X-rays \n\n[radda](https://www.kaggle.com/c/ranzcr-clip-catheter-line-classification/discussion/207183) kindly provided lung masks for the RANZCR training data. Lung masks are believed to be critical in order to successfully detect intubation/catheter malpositions.\nAs suggested I build my own UNet model, that can map the x-ray chest of the competition data to lung masks.\n\nThis notebook illustrates a simple custom Keras model to learn the lung-mask and largely follows [Peter Grenholm's ](https://www.kaggle.com/toregil/a-lung-u-net-in-keras) structure. The network is trained from scratch and does not use imagenet weights. I am not quite satisfied with the performance, therefore next I will try using qubvel's segmentation-model-keras that leverage pretrained Unet-models and work nicely with the image augmentation library albumentation.\n\nUpdated and improved notebook using transfer learning [here](https://www.kaggle.com/philippschwarz/ranzcr-lung-mask-transfer-learning)"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd \nimport cv2\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom sklearn.model_selection import train_test_split\nfrom pathlib import Path\nimport ast\nfrom tqdm import tqdm_notebook, tqdm","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"from keras.models import Model\nfrom keras.layers import *\nfrom keras.optimizers import Adam\nfrom keras.regularizers import l2\nfrom keras.preprocessing.image import ImageDataGenerator\nimport keras.backend as K\nfrom keras.callbacks import  CSVLogger, ModelCheckpoint, LearningRateScheduler, EarlyStopping, ReduceLROnPlateau\nfrom keras.losses import binary_crossentropy","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"TRAIN_PATH = '../input/ranzcr-clip-catheter-line-classification/train/'\nIMAGE_LIB = TRAIN_PATH\nIMG_HEIGHT, IMG_WIDTH = 128, 128\nSEED=42\nNUM_SAMPLES = 4000\nBATCH_SIZE = 32\nEPOCHS = 20\nctr = pd.read_csv('../input/ranzcr-clip-lung-contours/RANZCR_CLiP_lung_contours.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Helper Function to read masks"},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_mask(StudyInstanceUID):\n    img = cv2.imread(IMAGE_LIB+StudyInstanceUID+'.jpg',-1)\n    ctr_left = ast.literal_eval(ctr.loc[ctr.StudyInstanceUID==StudyInstanceUID,'left_lung_contour'].values[0])\n    ctr_right = ast.literal_eval(ctr.loc[ctr.StudyInstanceUID==StudyInstanceUID,'right_lung_contour'].values[0])\n    img = cv2.drawContours(img, np.array([[np.array(x) for x in ctr_left]]), 0, (255), -1)\n    img = cv2.drawContours(img, np.array([[np.array(x) for x in ctr_right]]), 0, (255), -1)\n    img = np.where(img>=255, 1.0, 0.0)\n    return img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"radda  explained how to load the masks in this [notebook](https://www.kaggle.com/raddar/simple-lung-contour-visualization). I updated the script such that the background is encoded 0 and the mask is 1."},{"metadata":{"trusted":true},"cell_type":"code","source":"all_images = os.listdir(TRAIN_PATH)[:NUM_SAMPLES]\nall_images = [Path(e).stem for e in all_images]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load images and masks into memory"},{"metadata":{"trusted":true},"cell_type":"code","source":"x_data = np.empty((len(all_images), IMG_HEIGHT, IMG_WIDTH), dtype='float32')\nfor i, name in enumerate(tqdm(all_images)):\n    im = cv2.imread(IMAGE_LIB + name +'.jpg', cv2.IMREAD_UNCHANGED).astype(\"int16\").astype('float32')\n    im = cv2.resize(im, dsize=(IMG_WIDTH, IMG_HEIGHT), interpolation=cv2.INTER_LANCZOS4)\n    im = (im - np.min(im)) / (np.max(im) - np.min(im))\n    x_data[i] = im\n\ny_data = np.empty((len(all_images), IMG_HEIGHT, IMG_WIDTH), dtype='float32')\nfor i, name in enumerate(tqdm(all_images)):\n    im = load_mask(name)\n    im = cv2.resize(im, dsize=(IMG_WIDTH, IMG_HEIGHT), interpolation=cv2.INTER_NEAREST)\n    y_data[i] = im\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Verify visually that images and masks are correct"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,2, figsize = (8,4))\nax[0].imshow(x_data[0], cmap='gray')\nax[1].imshow(y_data[0], cmap='gray')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"x_data = x_data[:,:,:,np.newaxis]\ny_data = y_data[:,:,:,np.newaxis]\nx_train, x_val, y_train, y_val = train_test_split(x_data, y_data, test_size = 0.5)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def dice_coef(y_true, y_pred):\n    y_true_f = K.flatten(y_true)\n    y_pred = K.cast(y_pred, 'float32')\n    y_pred_f = K.cast(K.greater(K.flatten(y_pred), 0.5), 'float32')\n    intersection = y_true_f * y_pred_f\n    score = 2. * K.sum(intersection) / (K.sum(y_true_f) + K.sum(y_pred_f))\n    return score\n\ndef dice_loss(y_true, y_pred):\n    smooth = 1.\n    y_true_f = K.flatten(y_true)\n    y_pred_f = K.flatten(y_pred)\n    intersection = y_true_f * y_pred_f\n    score = (2. * K.sum(intersection) + smooth) / (K.sum(y_true_f) + K.sum(y_pred_f) + smooth)\n    return 1. - score\n\ndef bce_dice_loss(y_true, y_pred):\n    return binary_crossentropy(y_true, y_pred) + dice_loss(y_true, y_pred)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"input_layer = Input(shape=x_train.shape[1:])\nc1 = Conv2D(filters=8, kernel_size=(3,3), activation='relu', padding='same')(input_layer)\nl = MaxPool2D(strides=(2,2))(c1)\nc2 = Conv2D(filters=16, kernel_size=(3,3), activation='relu', padding='same')(l)\nl = MaxPool2D(strides=(2,2))(c2)\nc3 = Conv2D(filters=32, kernel_size=(3,3), activation='relu', padding='same')(l)\nl = MaxPool2D(strides=(2,2))(c3)\nc4 = Conv2D(filters=32, kernel_size=(1,1), activation='relu', padding='same')(l)\nl = concatenate([UpSampling2D(size=(2,2))(c4), c3], axis=-1)\nl = Conv2D(filters=32, kernel_size=(2,2), activation='relu', padding='same')(l)\nl = concatenate([UpSampling2D(size=(2,2))(l), c2], axis=-1)\nl = Conv2D(filters=24, kernel_size=(2,2), activation='relu', padding='same')(l)\nl = concatenate([UpSampling2D(size=(2,2))(l), c1], axis=-1)\nl = Conv2D(filters=16, kernel_size=(2,2), activation='relu', padding='same')(l)\nl = Conv2D(filters=64, kernel_size=(1,1), activation='relu')(l)\nl = Dropout(0.5)(l)\noutput_layer = Conv2D(filters=1, kernel_size=(1,1), activation='sigmoid')(l)\n                                                         \nmodel = Model(input_layer, output_layer)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def my_generator(x_train, y_train, batch_size):\n    data_generator = ImageDataGenerator(\n            width_shift_range=0.1,\n            height_shift_range=0.1,\n            rotation_range=10,\n            zoom_range=0.1).flow(x_train, x_train, batch_size, seed=SEED)\n    mask_generator = ImageDataGenerator(\n            width_shift_range=0.1,\n            height_shift_range=0.1,\n            rotation_range=10,\n            zoom_range=0.1).flow(y_train, y_train, batch_size, seed=SEED)\n    while True:\n        x_batch, _ = data_generator.next()\n        y_batch, _ = mask_generator.next()\n        yield x_batch, y_batch\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_batch, mask_batch = next(my_generator(x_train, y_train, 8))\nfix, ax = plt.subplots(8,2, figsize=(8,20))\nfor i in range(8):\n    ax[i,0].imshow(image_batch[i,:,:,0])\n    ax[i,1].imshow(mask_batch[i,:,:,0])\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.compile(optimizer=Adam(2e-4), loss=bce_dice_loss, metrics=[dice_coef, binary_crossentropy])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"\nearly_stopping = EarlyStopping(patience=10, verbose=1, monitor='val_dice_coeff', mode='max')\nmodel_checkpoint = ModelCheckpoint(\"unet_custom_128-128_{epoch:02d}-{val_loss:.3f}.hdf5\", \n#                                    save_best_only=True, \n                                   save_weights_only=True, \n                                   monitor='val_dice_coeff', verbose=1, mode='max', period=2)\nreduce_lr = ReduceLROnPlateau(factor=0.5, patience=5, min_lr=0.000001, verbose=1, monitor='val_dice_coeff', mode='max')\n\nhist = model.fit_generator(my_generator(x_train, y_train, batch_size = BATCH_SIZE),\n                           steps_per_epoch = NUM_SAMPLES//BATCH_SIZE,\n                           validation_data = (x_val, y_val),\n                           epochs=EPOCHS,  \n                           callbacks=[ reduce_lr, model_checkpoint], # early_stopping\n                           verbose=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import json \n\nclass MyJsonEncoder(json.JSONEncoder):\n    def default(self, obj):\n        #if isinstance(obj, np.integer):\n        #    return int(obj)\n        if isinstance(obj, np.floating):\n            return float(obj)\n        #if isinstance(obj, np.ndarray):\n        #    return obj.tolist()\n        return super(MyJsonEncoder, self).default(obj)\n\n\nwith open('history.json', 'w') as f:\n    json.dump(hist.history, f, cls=MyJsonEncoder)\n    \nhistory_df = pd.DataFrame(hist.history)\nhistory_df.head(2)\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(1,3,figsize=(20,4))\nhistory_df.val_loss.plot(ax=ax[0], color='red', title='Validation Loss',ylim=(0,5))\nhistory_df.val_dice_coef.plot(ax=ax[1], color='blue', title='Validation binary_crossentropy', )\nhistory_df.val_binary_crossentropy.plot(ax=ax[2], color='green', title='Validation Dice_Coef');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(model.predict(x_train[2].reshape(1,IMG_HEIGHT, IMG_WIDTH, 1))[0,:,:,0], cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"n = 10\ny_hat = model.predict(x_val)\nfig, ax = plt.subplots(1,3,figsize=(12,6))\nax[0].imshow(x_val[n,:,:,0], cmap='gray')\nax[1].imshow(y_val[n,:,:,0])\nax[2].imshow(y_hat[n,:,:,0]);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"1. ### Test quality of masks on test dataset"},{"metadata":{"trusted":true},"cell_type":"code","source":"TEST_PATH = '../input/ranzcr-clip-catheter-line-classification/test/'\nname = '1.2.826.0.1.3680043.8.498.10023042737818625910026668901358652653'\nim = cv2.imread(TEST_PATH + name +'.jpg', cv2.IMREAD_UNCHANGED).astype(\"int16\").astype('float32')\nim = cv2.resize(im, dsize=(IMG_WIDTH, IMG_HEIGHT), interpolation=cv2.INTER_LANCZOS4)\nim = (im - np.min(im)) / (np.max(im) - np.min(im))\nim = im.reshape(1,IMG_WIDTH, IMG_HEIGHT, 1 )\n\ny_hat = model.predict(im)\nfig, ax = plt.subplots(1,2,figsize=(12,4))\nax[0].imshow(im[0,:,:,0], cmap='gray')\nax[1].imshow(y_hat[0,:,:,0]);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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":4}