{"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":"### Pretrained UNET Lung Segmentation Dataset\n\nThe objective of this notebook is quite self explanatory :   \n<center><b><h2>Create a Lung Masks Segmented Dataset</h2></b></center>\n\n### Why Lung Segmentation?\n\nMost of the COVID-19 images are hastily obtained in this situation.  \nHence, more often then not, they contain **Medical Noise**. Obviously, we can't blame doctors for this situation as they are working hard day and night to treat patients.\n\nWell then, **What we as Machine Learning Engineers can do about this??**\n\nWe can apply segmentation models for cleaning the dataset first! In this notebook, I apply my favourite model, the UNET for Segmenting and separating out the Lung masks from the image. You can read more about UNets here : [Official Paper](https://arxiv.org/abs/1505.04597) or [Blog](https://towardsdatascience.com/understanding-semantic-segmentation-with-unet-6be4f42d4b47).\n\nAlso I have a pretrained model ready for use here, borrowed from [Nikhil Pandey](https://www.kaggle.com/nikhilpandey360/lung-segmentation-from-chest-x-ray-dataset/output).","metadata":{}},{"cell_type":"markdown","source":"### Import Helper Functions","metadata":{}},{"cell_type":"code","source":"! conda install -c conda-forge gdcm -y","metadata":{"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nfrom glob import glob\nfrom tqdm.notebook import tqdm\nimport matplotlib.pyplot as plt\nfrom skimage import exposure\nimport cv2\nimport warnings\nwarnings.filterwarnings('ignore')\nimport shutil \nimport tensorflow as tf\n%matplotlib inline\n\n\nimport matplotlib.pylab as pylab\nimport seaborn as sns\nimport pprint\nimport pydicom as dicom\nfrom pydicom.pixel_data_handlers.util import apply_voi_lut\nimport wandb\n\nimport PIL\nfrom PIL import Image\nfrom colorama import Fore, Back, Style\nviz_counter=0\n\ndef create_dir(dir, v=1):\n    \"\"\"\n    Creates a directory without throwing an error if directory already exists.\n    dir : The directory to be created.\n    v : Verbosity\n    \"\"\"\n    if not os.path.exists(dir):\n        os.makedirs(dir)\n        if v:\n            print(\"Created Directory : \", dir)\n        return 1\n    else:\n        if v:\n            print(\"Directory already existed : \", dir)\n        return 0\n\nvoi_lut=True\nfix_monochrome=True\n\ndef dicom_dataset_to_dict(filename):\n    \"\"\"Credit: https://github.com/pydicom/pydicom/issues/319\n               https://www.kaggle.com/raddar/convert-dicom-to-np-array-the-correct-way\n    \"\"\"\n    \n    dicom_header = dicom.dcmread(filename) \n    \n    #====== DICOM FILE DATA ======\n    dicom_dict = {}\n    repr(dicom_header)\n    for dicom_value in dicom_header.values():\n        if dicom_value.tag == (0x7fe0, 0x0010):\n            #discard pixel data\n            continue\n        if type(dicom_value.value) == dicom.dataset.Dataset:\n            dicom_dict[dicom_value.name] = dicom_dataset_to_dict(dicom_value.value)\n        else:\n            v = _convert_value(dicom_value.value)\n            dicom_dict[dicom_value.name] = v\n      \n    del dicom_dict['Pixel Representation']\n    \n    #====== DICOM IMAGE DATA ======\n    # VOI LUT (if available by DICOM device) is used to transform raw DICOM data to \"human-friendly\" view\n    if voi_lut:\n        data = apply_voi_lut(dicom_header.pixel_array, dicom_header)\n    else:\n        data = dicom_header.pixel_array\n    # depending on this value, X-ray may look inverted - fix that:\n    if fix_monochrome and dicom_header.PhotometricInterpretation == \"MONOCHROME1\":\n        data = np.amax(data) - data\n    data = data - np.min(data)\n    data = data / np.max(data)\n    modified_image_data = (data * 255).astype(np.uint8)\n    \n    return dicom_dict, modified_image_data\n\ndef _sanitise_unicode(s):\n    return s.replace(u\"\\u0000\", \"\").strip()\n\ndef _convert_value(v):\n    t = type(v)\n    if t in (list, int, float):\n        cv = v\n    elif t == str:\n        cv = _sanitise_unicode(v)\n    elif t == bytes:\n        s = v.decode('ascii', 'replace')\n        cv = _sanitise_unicode(s)\n    elif t == dicom.valuerep.DSfloat:\n        cv = float(v)\n    elif t == dicom.valuerep.IS:\n        cv = int(v)\n    else:\n        cv = repr(v)\n    return cv\n\n\nimport os, fnmatch\ndef find(pattern, path):\n    \"\"\"Utility to find files wrt a regex search\"\"\"\n    result = []\n    for root, dirs, files in os.walk(path):\n        for name in files:\n            if fnmatch.fnmatch(name, pattern):\n                result.append(os.path.join(root, name))\n    return result\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets get all the .dcm files","metadata":{}},{"cell_type":"code","source":"FIND_FOLDER=\"/kaggle/input/siim-covid19-detection\"\ndcm_files=find('*.dcm', FIND_FOLDER)\nprint(len(dcm_files),\"Files Found.\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import random\nrandom.seed(52) # 42 # 2021\nsubset_dcm_files =  random.choices(dcm_files, k=3) # dcm_files[:3]\nsubset_dcm_files","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"viz_counter=0","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def props(arr):\n    print(\"Shape :\",arr.shape,\"Maximum :\",arr.max(),\"Minimum :\",arr.min(),\"Data Type :\",arr.dtype)\nfor path in subset_dcm_files:\n    dicom_dict, modified_image_data = dicom_dataset_to_dict(path)\n    props(modified_image_data)\n    # print(dicom_dict)\n    fig, ax = plt.subplots(1, 2, figsize=(20, 12))\n    ax[0].imshow(modified_image_data, cmap=\"gray\")\n    ax[0].axis('off')\n    ax[1].imshow(modified_image_data, cmap=\"viridis\")    \n    ax[1].axis('off')\n    plt.savefig(str(viz_counter)+\".png\",dpi=300)\n    viz_counter+=1\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define UNET Model","metadata":{}},{"cell_type":"code","source":"from keras.models import *\nfrom keras.layers import *\nfrom keras.optimizers import *\nfrom keras import backend as keras\nfrom keras.preprocessing.image import ImageDataGenerator\nfrom keras.callbacks import ModelCheckpoint, LearningRateScheduler\n\n\ndef dice_coef(y_true, y_pred):\n    y_true_f = keras.flatten(y_true)\n    y_pred_f = keras.flatten(y_pred)\n    intersection = keras.sum(y_true_f * y_pred_f)\n    return (2. * intersection + 1) / (keras.sum(y_true_f) + keras.sum(y_pred_f) + 1)\n\ndef dice_coef_loss(y_true, y_pred):\n    return -dice_coef(y_true, y_pred)\n\ndef unet(input_size=(256,256,1)):\n    inputs = Input(input_size)\n    \n    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(inputs)\n    conv1 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv1)\n    pool1 = MaxPooling2D(pool_size=(2, 2))(conv1)\n\n    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(pool1)\n    conv2 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv2)\n    pool2 = MaxPooling2D(pool_size=(2, 2))(conv2)\n\n    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(pool2)\n    conv3 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv3)\n    pool3 = MaxPooling2D(pool_size=(2, 2))(conv3)\n\n    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(pool3)\n    conv4 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv4)\n    pool4 = MaxPooling2D(pool_size=(2, 2))(conv4)\n\n    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(pool4)\n    conv5 = Conv2D(512, (3, 3), activation='relu', padding='same')(conv5)\n\n    up6 = concatenate([Conv2DTranspose(256, (2, 2), strides=(2, 2), padding='same')(conv5), conv4], axis=3)\n    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(up6)\n    conv6 = Conv2D(256, (3, 3), activation='relu', padding='same')(conv6)\n\n    up7 = concatenate([Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same')(conv6), conv3], axis=3)\n    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(up7)\n    conv7 = Conv2D(128, (3, 3), activation='relu', padding='same')(conv7)\n\n    up8 = concatenate([Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same')(conv7), conv2], axis=3)\n    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(up8)\n    conv8 = Conv2D(64, (3, 3), activation='relu', padding='same')(conv8)\n\n    up9 = concatenate([Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same')(conv8), conv1], axis=3)\n    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(up9)\n    conv9 = Conv2D(32, (3, 3), activation='relu', padding='same')(conv9)\n\n    conv10 = Conv2D(1, (1, 1), activation='sigmoid')(conv9)\n\n    return Model(inputs=[inputs], outputs=[conv10])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = unet(input_size=(512,512,1))\nmodel.compile(optimizer=Adam(lr=1e-5), loss=dice_coef_loss,\n                  metrics=[dice_coef, 'binary_accuracy'])\nmodel.summary()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load the Pretrained UNet Model","metadata":{}},{"cell_type":"code","source":"model_weights_path = \"/kaggle/input/unet-lung-segmentation-weights-for-chest-x-rays/cxr_reg_weights.best.hdf5\"\n\nmodel.load_weights(model_weights_path)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nShapes that you wish to resize to\n\"\"\"\n\nShape_X,Shape_Y=512,512\n\nfor path in subset_dcm_files:\n    dicom_dict, modified_image_data = dicom_dataset_to_dict(path)\n    resized_image_data = cv2.resize(modified_image_data,(Shape_Y,Shape_X)) # cv2 has this opposite\n    # props(resized_image_data)\n    prep_unet_input_img_1 = resized_image_data.reshape(1,Shape_X,Shape_Y,1)\n    prep_unet_input_img = (prep_unet_input_img_1-127.0)/127.0\n    pred_img = model.predict(prep_unet_input_img)\n    pred_img_preprocessed_1 = np.squeeze(pred_img)\n    pred_img_preprocessed = (pred_img_preprocessed_1*255>127).astype(np.int8)\n    # props(pred_img_preprocessed)\n    # print(\"Unique Values :\",np.unique(pred_img_preprocessed))\n    res = cv2.bitwise_and(resized_image_data,resized_image_data,mask = pred_img_preprocessed)\n    fig, ax = plt.subplots(1, 3, figsize=(20, 12))\n    ax[0].imshow(resized_image_data, cmap=\"viridis\")\n    ax[0].axis('off')\n    ax[1].imshow(pred_img_preprocessed, cmap=\"viridis\")    \n    ax[1].axis('off')\n    ax[2].imshow(res, cmap=\"viridis\")    \n    ax[2].axis('off')\n    plt.savefig(str(viz_counter)+\".png\",dpi=300)\n    viz_counter+=1\n    cv2.imwrite(str(viz_counter)+\".png\",res)\n    viz_counter+=1\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_id = []\ndim0 = []\ndim1 = []\nsplits = []\n\nfor split in ['test', 'train']:\n    # save_dir = f'/kaggle/tmp/{split}/'\n    save_dir = f'/kaggle/working/segmented_data/{split}/'\n    print(split)\n    os.makedirs(save_dir, exist_ok=True)\n    \n    for dirname, _, filenames in tqdm(os.walk(f'../input/siim-covid19-detection/{split}')):\n        for file in filenames:\n            # set keep_ratio=True to have original aspect ratio\n            fpath = os.path.join(dirname, file)\n            dicom_dict, modified_image_data = dicom_dataset_to_dict(fpath)\n            resized_image_data = cv2.resize(modified_image_data,(Shape_Y,Shape_X)) # cv2 has this opposite\n            # props(resized_image_data)\n            prep_unet_input_img_1 = resized_image_data.reshape(1,Shape_X,Shape_Y,1)\n            prep_unet_input_img = (prep_unet_input_img_1-127.0)/127.0\n            pred_img = model.predict(prep_unet_input_img)\n            pred_img_preprocessed_1 = np.squeeze(pred_img)\n            pred_img_preprocessed = (pred_img_preprocessed_1*255>127).astype(np.int8)\n            # props(pred_img_preprocessed)\n            # print(\"Unique Values :\",np.unique(pred_img_preprocessed))\n            res = cv2.bitwise_and(resized_image_data,resized_image_data,mask = pred_img_preprocessed)\n            save_path = os.path.join(save_dir, file.replace('dcm', 'png'))\n            cv2.imwrite(save_path,res)\n\n            image_id.append(file.replace('.dcm', ''))\n            dim0.append(res.shape[0])\n            dim1.append(res.shape[1])\n            splits.append(split)\n\"\"\"\n2475/?\n12386/?\n07:34 | 5.38it/s\n36:51 | 8.13it/s\n\"\"\"\nprint(\"Generation Complete!\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame.from_dict({'image_id': image_id, 'dim0': dim0, 'dim1': dim1, 'split': splits})\ndf.to_csv('meta.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport zipfile\nimport shutil\n\n#taken from : https://www.kaggle.com/xhlulu/recursion-2019-load-resize-and-save-images\n\ndef zip_and_remove(path):\n    ziph = zipfile.ZipFile(f'{path}.zip', 'w', zipfile.ZIP_DEFLATED)\n    \n    for root, dirs, files in os.walk(path):\n        for file in files:\n            file_path = os.path.join(root, file)\n            ziph.write(file_path)\n            os.remove(file_path)\n    \n    ziph.close()\n    shutil.rmtree(path)\nsave_dir = 'segmented_data'\nzip_and_remove(save_dir)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}