{"cells":[{"metadata":{"trusted":true,"collapsed":true,"_uuid":"c7773c37ecf955e2b06f460403bac3ec864a6b57"},"cell_type":"markdown","source":"#  In this notebook I experiment with different image enhancement techniques that may be used to form features for cnns. \n# I also extract features from three popular cnn architectures namely VGG, Resnet and Xception to see if these features make sense in the context of Pneumonia Detection. \n"},{"metadata":{"_uuid":"e6e7f39dbcc5c2fcf94bce04f4c2d3005943fb89"},"cell_type":"markdown","source":"# Loading Libraries"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"collapsed":true},"cell_type":"code","source":"import os\nimport sys\nimport random\n\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n%matplotlib inline\n\nfrom skimage.transform import resize\nfrom skimage.morphology import label\nfrom skimage.feature import hog\nfrom skimage import exposure\nfrom keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img\nfrom skimage.feature import canny\nfrom skimage.filters import sobel\nfrom skimage.morphology import watershed\nfrom scipy import ndimage as ndi\nimport warnings\nwarnings.filterwarnings(\"ignore\")\nfrom skimage.segmentation import mark_boundaries\nfrom scipy import signal\nimport cv2\nimport glob, pylab, pandas as pd\nimport pydicom, numpy as np\nimport tqdm","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fbbd5ce0a95297597f388fad655cbfbb8bd8457f"},"cell_type":"markdown","source":"# Helper Functions"},{"metadata":{"trusted":true,"_uuid":"10feec03b44374c513d156b4a0b165adc32b945a","collapsed":true},"cell_type":"code","source":"# https://www.kaggle.com/peterchang77/exploratory-data-analysis\ndef parse_data(df):\n    \"\"\"\n    Method to read a CSV file (Pandas dataframe) and parse the \n    data into the following nested dictionary:\n\n      parsed = {\n        \n        'patientId-00': {\n            'dicom': path/to/dicom/file,\n            'label': either 0 or 1 for normal or pnuemonia, \n            'boxes': list of box(es)\n        },\n        'patientId-01': {\n            'dicom': path/to/dicom/file,\n            'label': either 0 or 1 for normal or pnuemonia, \n            'boxes': list of box(es)\n        }, ...\n\n      }\n\n    \"\"\"\n    # --- Define lambda to extract coords in list [y, x, height, width]\n    extract_box = lambda row: [row['y'], row['x'], row['height'], row['width']]\n\n    parsed = {}\n    for n, row in df.iterrows():\n        # --- Initialize patient entry into parsed \n        pid = row['patientId']\n        if pid not in parsed:\n            parsed[pid] = {\n                'dicom': '../input/stage_1_train_images/%s.dcm' % pid,\n                'label': row['Target'],\n                'boxes': []}\n\n        # --- Add box if opacity is present\n        if parsed[pid]['label'] == 1:\n            parsed[pid]['boxes'].append(extract_box(row))\n\n    return parsed","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8f5e52758c34c1d1e2b48f836edd8c6d6dcafe19","collapsed":true},"cell_type":"code","source":"# https://www.kaggle.com/peterchang77/exploratory-data-analysis\ndef draw(data,im):\n    \"\"\"\n    Method to draw single patient with bounding box(es) if present \n\n    \"\"\"\n\n    # --- Convert from single-channel grayscale to 3-channel RGB\n    im = np.stack([im] * 3, axis=2)\n\n    # --- Add boxes with random color if present\n    for box in data['boxes']:\n        rgb = np.floor(np.random.rand(3) * 256).astype('int')\n        im = overlay_box(im=im, box=box, rgb=rgb, stroke=6)\n        \n    return im\n\ndef overlay_box(im, box, rgb, stroke=1):\n    \"\"\"\n    Method to overlay single box on image\n\n    \"\"\"\n    # --- Convert coordinates to integers\n    box = [int(b) for b in box]\n    \n    # --- Extract coordinates\n    y1, x1, height, width = box\n    y2 = y1 + height\n    x2 = x1 + width\n\n    im[y1:y1 + stroke, x1:x2] = rgb\n    im[y2:y2 + stroke, x1:x2] = rgb\n    im[y1:y2, x1:x1 + stroke] = rgb\n    im[y1:y2, x2:x2 + stroke] = rgb\n\n    return im","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"abc3b9684e2451b7402ce159dc843e78d44e53ab"},"cell_type":"code","source":"df = pd.read_csv('../input/stage_1_train_labels.csv')\nparsed = parse_data(df)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"c3b9ffcebc7d64a4a6b441b72e8aa493ed58fc6c","collapsed":true},"cell_type":"code","source":"det_class_path = '../input/stage_1_detailed_class_info.csv'\ndet_class_df = pd.read_csv(det_class_path)\ndet_class_df.head()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0cfa8e69c78b9ab314ed963755df1613b2806d3b"},"cell_type":"markdown","source":"It can be seen that there are three classes namely \n- No Lung Opacity / Not Normal\n- Normal\n- Lung Opacity\n\nLet us visualize each class separately to get better idea."},{"metadata":{"_uuid":"204d3e72cdbf383ac5ade828bb393ad4cd95d11c"},"cell_type":"markdown","source":"# No Lung Opacity / Not Normal"},{"metadata":{"trusted":true,"_uuid":"480bacbe1e24dc42f480bfb62f7fa9cc4b5d6112","collapsed":true},"cell_type":"code","source":"# simple features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\n\nplt.figure(figsize=(30,15))\n# plt.subplots_adjust(bottom=0.2, top=0.8, hspace=0)  #adjust this to change vertical and horiz. spacings..\nnImg = 3  #no. of images to process\nj = -1\ndf = det_class_df[det_class_df['class']=='No Lung Opacity / Not Normal']\ndf = df.reset_index()\nwhile True:\n# for j in range(nImg):\n    if j == nImg-1:\n        break\n        \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    \n    data = parsed[patientId]\n    j += 1\n        \n    q = j+1\n    \n#     # Contrast stretching\n    p2, p98 = np.percentile(img, (2, 98))\n    img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))\n    \n    # Equalization\n    img_eq = exposure.equalize_hist(img)\n\n    # Adaptive Equalization\n    img_adapteq = exposure.equalize_adapthist(img)\n\n    plt.subplot(nImg,5,q*5-4)\n    plt.imshow(draw(data,img), cmap='binary')\n    plt.title('Original Image')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-3)    \n    plt.imshow(draw(data,img_rescale), cmap='binary')\n    plt.title('Contrast stretching')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-2)\n    plt.imshow(draw(data,img_eq), cmap='binary')\n    plt.title('Equalization')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-1)\n    plt.imshow(draw(data,img_adapteq), cmap='binary')\n    plt.title('Adaptive Equalization')\n    plt.axis('off')\n\nplt.show()\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"ad795a28989e10e8e1c1e0448bff42844386dd8a"},"cell_type":"markdown","source":"# Normal"},{"metadata":{"trusted":true,"_uuid":"e9b47e532dd011453ebb6e28a069e881a7624043","collapsed":true},"cell_type":"code","source":"# simple features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\n\nplt.figure(figsize=(30,15))\n# plt.subplots_adjust(bottom=0.2, top=0.8, hspace=0)  #adjust this to change vertical and horiz. spacings..\nnImg = 3  #no. of images to process\nj = -1\ndf = det_class_df[det_class_df['class']=='Normal']\ndf = df.reset_index()\nwhile True:\n# for j in range(nImg):\n    if j == nImg-1:\n        break\n        \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    \n    data = parsed[patientId]\n    j += 1\n        \n    q = j+1\n    \n#     # Contrast stretching\n    p2, p98 = np.percentile(img, (2, 98))\n    img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))\n    \n    # Equalization\n    img_eq = exposure.equalize_hist(img)\n\n    # Adaptive Equalization\n    img_adapteq = exposure.equalize_adapthist(img)\n\n    plt.subplot(nImg,5,q*5-4)\n    plt.imshow(draw(data,img), cmap='binary')\n    plt.title('Original Image')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-3)    \n    plt.imshow(draw(data,img_rescale), cmap='binary')\n    plt.title('Contrast stretching')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-2)\n    plt.imshow(draw(data,img_eq), cmap='binary')\n    plt.title('Equalization')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-1)\n    plt.imshow(draw(data,img_adapteq), cmap='binary')\n    plt.title('Adaptive Equalization')\n    plt.axis('off')\n\nplt.show()\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b00022249c6ac958c32b9294e3392e4dde8afc5e"},"cell_type":"markdown","source":"# Lung Opacity"},{"metadata":{"trusted":true,"_uuid":"7da75fea003b23ef97b4295d0f59b80cda4b0301","collapsed":true},"cell_type":"code","source":"# simple features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\n\nplt.figure(figsize=(30,15))\n# plt.subplots_adjust(bottom=0.2, top=0.8, hspace=0)  #adjust this to change vertical and horiz. spacings..\nnImg = 3  #no. of images to process\nj = -1\ndf = det_class_df[det_class_df['class']=='Lung Opacity']\ndf = df.reset_index()\nwhile True:\n# for j in range(nImg):\n    if j == nImg-1:\n        break\n        \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    \n    data = parsed[patientId]\n    j += 1\n        \n    q = j+1\n    \n#     # Contrast stretching\n    p2, p98 = np.percentile(img, (2, 98))\n    img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))\n    \n    # Equalization\n    img_eq = exposure.equalize_hist(img)\n\n    # Adaptive Equalization\n    img_adapteq = exposure.equalize_adapthist(img)\n\n    plt.subplot(nImg,5,q*5-4)\n    plt.imshow(draw(data,img), cmap='binary')\n    plt.title('Original Image')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-3)    \n    plt.imshow(draw(data,img_rescale), cmap='binary')\n    plt.title('Contrast stretching')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-2)\n    plt.imshow(draw(data,img_eq), cmap='binary')\n    plt.title('Equalization')\n    plt.axis('off')\n    \n    plt.subplot(nImg,5,q*5-1)\n    plt.imshow(draw(data,img_adapteq), cmap='binary')\n    plt.title('Adaptive Equalization')\n    plt.axis('off')\n\nplt.show()\nplt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"8083c9a67b592f0191a7a326b159ba446b7f4257","collapsed":true},"cell_type":"code","source":"from keras.applications.vgg16 import VGG16\nfrom keras.applications.xception import Xception\nfrom keras.applications.resnet50 import ResNet50\n\nfrom keras.preprocessing import image\nfrom keras.applications.vgg16 import preprocess_input\nfrom keras import backend as K\nfrom keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"acf931a63b8a5392e10a7cada54aad9a1ddde837"},"cell_type":"markdown","source":"# Extracting and Visualizing VGG features"},{"metadata":{"trusted":true,"_uuid":"5d50b0c47cd34e39b4d4242a468b403c51183d18","collapsed":true},"cell_type":"code","source":"model = VGG16(weights='imagenet', include_top=False)\ninp = model.input                                           # input placeholder\noutputs = [layer.output for layer in model.layers]          # all layer outputs\nfunctors = [K.function([inp]+ [K.learning_phase()], [out]) for out in outputs]  # evaluation functions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"2c2cde6b1ca2eb13b0140c2d0bf492150e1f595e","collapsed":true},"cell_type":"code","source":"# vgg features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\nrandom.seed(40)\nplt.figure(figsize=(30,15))\nplt.subplots_adjust(bottom=0.2, top=0.8, hspace=0.2)  #adjust this to change vertical and horiz. spacings..\nnImg = 3  #no. of images to process\nfor j in range(nImg):\n    q = j+1\n    \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    img = cv2.resize(img,(224, 224))\n    img = np.expand_dims(img, axis=-1)\n    img = np.repeat(img,3,axis=2)\n    x = image.img_to_array(img)\n    x = np.expand_dims(x, axis=0)\n    x = preprocess_input(x)\n\n    \n    layer_outs = [func([x, 0.]) for func in functors]\n    feat = np.reshape(layer_outs[4][0],(112,112,128))\n    layer4 = np.max(feat,axis=2)\n    \n    feat = np.reshape(layer_outs[6][0],(56,56,128))\n    layer6 = np.max(feat,axis=2)\n    \n    feat = np.reshape(layer_outs[10][0],(28,28,256))\n    layer10 = np.max(feat,axis=2)\n    \n    plt.subplot(nImg,6,q*6-5)\n    plt.imshow(img, cmap='binary')\n    plt.title('Original Image')\n    \n    plt.subplot(nImg,6,q*6-4)\n    plt.imshow(img, cmap='binary')\n    plt.title('Image Mask')\n    \n    plt.subplot(nImg,6,q*6-3)    \n    plt.imshow(layer4, cmap='binary')\n    plt.title('VGG Layer 4')\n    \n    plt.subplot(nImg,6,q*6-2)\n    plt.imshow(layer6, cmap='binary')\n    plt.title('VGG Layer 6')\n    \n    plt.subplot(nImg,6,q*6-1)\n    plt.imshow(layer10, cmap='binary')\n    plt.title('VGG Layer 10')\n\n\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"88865731dffa2d70ce83078ebdd44bae56961fdf","collapsed":true},"cell_type":"code","source":"# model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"67ddae4ff49c873166c75364717f81d597ccb75f","collapsed":true},"cell_type":"code","source":"model = ResNet50(weights='imagenet',input_shape=(224, 224, 3), include_top=False)\ninp = model.input                                           # input placeholder\noutputs = [layer.output for layer in model.layers]          # all layer outputs\nfunctors = [K.function([inp]+ [K.learning_phase()], [out]) for out in outputs]  # evaluation functions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7cee19be4908d037e763176f69297d7c01aed409","collapsed":true},"cell_type":"code","source":"# model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"9ee2bcb1b00498197cd53c282a32b6a3284ecfb1","collapsed":true},"cell_type":"code","source":"# resnet features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\nrandom.seed(40)\nplt.figure(figsize=(15,30))\n# plt.subplots_adjust(bottom=0.2, top=0.8, hspace=0.2)  #adjust this to change vertical and horiz. spacings..\nnImg = 5  #no. of images to process\nfor j in range(nImg):\n    q = j+1\n    \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    img = cv2.resize(img,(224, 224))\n    img = np.expand_dims(img, axis=-1)\n    img = np.repeat(img,3,axis=2)\n    x = image.img_to_array(img)\n    x = np.expand_dims(x, axis=0)\n    x = preprocess_input(x)\n\n    \n    layer_outs = [func([x, 0.]) for func in functors]\n    feat = np.reshape(layer_outs[4][0],(112,112,64))\n    layer4 = np.max(feat,axis=2)\n    \n    plt.subplot(nImg,3,q*3-2)\n    plt.imshow(img, cmap='binary')\n    plt.title('Original Image')\n    plt.axis('off')\n    \n    plt.subplot(nImg,3,q*3-1)    \n    plt.imshow(layer4, cmap='binary')\n    plt.title('ResNet activation_1 ')\n    plt.axis('off')\n\n\nplt.show()\n# plt.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0dd3a973919480ef103f0ff600493d7e2f9d9467"},"cell_type":"markdown","source":"# Visualizing Xception features "},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"b976fe9fc0f5258f199a8d2fc8302a47af4b22ba"},"cell_type":"code","source":"model = Xception(weights='imagenet', include_top=False)\ninp = model.input                                           # input placeholder\noutputs = [layer.output for layer in model.layers]          # all layer outputs\nfunctors = [K.function([inp]+ [K.learning_phase()], [out]) for out in outputs]  # evaluation functions","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"7c7bc582fb0846509a54b9f6e6722442a1dcb715","collapsed":true},"cell_type":"code","source":"# Xception features that can be easily extracted and used for training deep networks\n# these features may be used along with original image\nrandom.seed(40)\nplt.figure(figsize=(30,15))\nplt.subplots_adjust(bottom=0.2, top=0.8, hspace=0.2)  #adjust this to change vertical and horiz. spacings..\nnImg = 3  #no. of images to process\nfor j in range(nImg):\n    q = j+1\n    \n    ind = np.random.randint(df.shape[0])\n    patientId = df['patientId'][ind]\n    dcm_file = '../input/stage_1_train_images/%s.dcm' % patientId\n    dcm_data = pydicom.read_file(dcm_file)\n    img = dcm_data.pixel_array\n    img = cv2.resize(img,(224, 224))\n    img = np.expand_dims(img, axis=-1)\n    img = np.repeat(img,3,axis=2)\n    x = image.img_to_array(img)\n    x = np.expand_dims(x, axis=0)\n    x = preprocess_input(x)\n\n    \n    layer_outs = [func([x, 0.]) for func in functors]\n    feat = np.reshape(layer_outs[4][0],(109,109,64))\n    layer4 = np.max(feat,axis=2)\n    \n    feat = np.reshape(layer_outs[16][0],(55,55,128))\n    layer6 = np.max(feat,axis=2)\n    \n    feat = np.reshape(layer_outs[26][0],(28,28,256))\n    layer10 = np.max(feat,axis=2)\n    \n    plt.subplot(nImg,6,q*6-5)\n    plt.imshow(img, cmap='binary')\n    plt.title('Original Image')\n    \n    plt.subplot(nImg,6,q*6-4)\n    plt.imshow(img, cmap='binary')\n    plt.title('Image Mask')\n    \n    plt.subplot(nImg,6,q*6-3)    \n    plt.imshow(layer4, cmap='binary')\n    plt.title('Xception Block 1')\n    \n    plt.subplot(nImg,6,q*6-2)\n    plt.imshow(layer6, cmap='binary')\n    plt.title('Xception Block 2')\n    \n    plt.subplot(nImg,6,q*6-1)\n    plt.imshow(layer10, cmap='binary')\n    plt.title('Xception Block 3 ')\n\n\nplt.show()","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}