{"cells":[{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Introduction\n\nTo maximize our search I think it might be advantageous to perform some segmentation on the dataset before applying and algorithm to identify bounding boxes.  I know this might be tedious but coming from the classical ML approach where by features are extracted from the data domain and classification is done on the feature space.  I think this approach is far more intuitve for me at the moment.  Therefore, I propose the following steps:\n\n1. Pre-Processing\n    - Masking the image based on the distribution of bounding boxes\n2. Feature Extraction\n    - Extract histogram as a feature for each of the test cases\n3. Classification\n    - Use classification techniques to classify normal vs. ( opaque, non-opaque) cases\n \n **This notebook will only focus on the pre-processing aspects of this workflow using another dataset to generate the lung segmentation portion.  We can further retrain the model on our model and see if it has improved. **"},{"metadata":{"_uuid":"88b6889047da86b866b4d4a1df9514a4a2ada163"},"cell_type":"markdown","source":"# Imports "},{"metadata":{"trusted":true,"_uuid":"bc033c6e9ad62dafd41b113060ee36651fbd9e46"},"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,sys\n\nimport pydicom\nimport tensorflow as tf\n\nfrom sklearn.utils import shuffle\nfrom scipy.ndimage import imread\nfrom scipy.misc import imresize\n\nfrom glob import glob\n\n%matplotlib inline\nimport matplotlib.pyplot as plt","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5900613426ce145b38dae36f232b4226ae8c6d9c"},"cell_type":"markdown","source":"# Pre-Processing Step"},{"metadata":{"_uuid":"ef292f56eaed82b1cb07f5d821eb5d684f6af9d6"},"cell_type":"markdown","source":"## Reading the Data\n\nLet us use the data uploaded by Kevin Mader [kernel](https://www.kaggle.com/kmader/gaussian-mixture-lung-segmentation) to train our segmentation U-Net algorithm.  It is noted that he has also done this but uses code that I cannot follow at the moment due to my lack of inexperiance but please refer to his guide."},{"metadata":{"trusted":true,"_uuid":"cef763a4f4b8596b3dba29c8d85f103dc4271bbb","scrolled":true},"cell_type":"code","source":"cxr_paths = glob(os.path.join('..', 'input', 'pulmonary-chest-xray-abnormalities','Montgomery', 'MontgomerySet', '*', '*.png'))\nprint(\"There are {} images in the dataset\".format(len(cxr_paths)))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a744761bea8a7c4e151ac4c2639a088f8658c97b"},"cell_type":"markdown","source":"Obtain the masks for each 138 patients. "},{"metadata":{"trusted":true,"_uuid":"42dbee18ef621864015cee55be55092e72b8aa5b"},"cell_type":"code","source":"cxr_images = [(c_path, \n               [os.path.join('/'.join(c_path.split('/')[:-2]),'ManualMask','leftMask', os.path.basename(c_path)),\n               os.path.join('/'.join(c_path.split('/')[:-2]),'ManualMask','rightMask', os.path.basename(c_path))]\n              ) for c_path in cxr_paths]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"867d1f92aa728fc90cdad0b807c74c745e8b60d3"},"cell_type":"markdown","source":"An example patient in the list has the structures below.  \n\ncxr_images[ **#ofPatients**][**0 = Location of the Patient Original Scan**][** 0/1 = Mask for the left and right lung respectively**]"},{"metadata":{"trusted":true,"_uuid":"fea0a158cbb83afa224eb08c7125f0e0f65c1a4e"},"cell_type":"code","source":"cxr_images[0]","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"0a0ab5057773b9dd1f06652f67627505c4f484d5"},"cell_type":"markdown","source":"## Reading Imports\n\nWe need more imports for reading the actual data"},{"metadata":{"_uuid":"d73f4b6623b3f17eb1725804176972b8888bc09e","trusted":true},"cell_type":"code","source":"from skimage.io import imread as imread_raw\nfrom skimage.transform import resize\nimport warnings\nfrom tqdm import tqdm\nwarnings.filterwarnings('ignore', category=UserWarning, module='skimage') # skimage is really annoying","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1631173c92b7044baafeb0e18e3de24ab6b9b71c"},"cell_type":"markdown","source":"These images are 512 x 512 images so let us create a global variable so that we can follow it."},{"metadata":{"_uuid":"b023989553c4e9d076376689bf522f62783e44e9"},"cell_type":"markdown","source":"## Function to read the images"},{"metadata":{"trusted":true,"_uuid":"a5da48046b696eeb0ef28e288f351f491c8d5862"},"cell_type":"code","source":"def imread(in_path):\n    OUT_DIM = (512, 512)\n    \n    # use the skimge function to read the file specified in the path\n    img_data = imread_raw(in_path)\n    \n    # make sure the image data is between the range 0-255 and convert the variable into uint8\n    n_img = (255*resize(img_data, OUT_DIM, mode = 'constant')).clip(0,255).astype(np.uint8)\n        \n    return np.expand_dims(n_img, -1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":true,"_uuid":"bdf33234370949b033275509e2b6f260ffec6a9d"},"cell_type":"code","source":"# init empty array for images and masks or in this case segmentations\nimg_vol, seg_vol = [], []\n\nfor img_path, s_paths in tqdm(cxr_images):\n    # first read the image paths\n    img_vol += [imread(img_path)]    \n    \n    # read both images, stack them up, then store them    \n    seg_vol += [np.max(np.stack([imread(s_path) for s_path in s_paths],0),0)]\n    \nimg_vol = np.stack(img_vol,0)\nseg_vol = np.stack(seg_vol,0)\n\nprint('Images', img_vol.shape, 'Segmentations', seg_vol.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"cc7abdff12df2e60990406a791651f415147b3ea"},"cell_type":"markdown","source":"As you can see from the above line there are 138 patients consisting of images with a 512x512 image and segmentation corresponding to the labeled lung segmentation.  Let us look at them side by side to get a better intuitive understading of the above images."},{"metadata":{"trusted":true,"_uuid":"4a47657e532818ef7e99efc0574dae82f02d112b"},"cell_type":"code","source":"# Get a random patient\nnp.random.seed(64)\n\nrandomPatient = int(np.random.rand()*138)\nt_img, m_img = img_vol[randomPatient], seg_vol[randomPatient]\n\n#plot it\nscan = t_img[:,:,0]\nmask = m_img[:,:,0]\nsegmented = scan*mask\n\ndef drawImage(ax_img,img,label):\n    ax_img.imshow(img,interpolation='none',cmap='bone')\n    ax_img.set_title(label)\n    ax_img.set_axis_off()\n    \n\nfig, (ax_img, ax_mask,ax_segmentedImage) = plt.subplots(1,3, figsize = (12, 6))\ndrawImage(ax_img,scan,label='CXR')\ndrawImage(ax_mask,mask,label='Labeled Mask')\ndrawImage(ax_segmentedImage,segmented,label='Segmented Image')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"6c06bea2346124768b1f03b65199e9ceae4b5893"},"cell_type":"markdown","source":"## Mask Generation"},{"metadata":{"_uuid":"f6839a3e6dd78f66f7ce295e0ae868ab3784beb2"},"cell_type":"markdown","source":"In this section we will use the [U-Net algorithm](https://arxiv.org/abs/1505.04597) for segmenting the lungs.\n\nFor some guidance we will use the following diagram from the article and also use the working sample by Jae Duk Seo in his [article](https://towardsdatascience.com/medical-image-segmentation-part-1-unet-convolutional-networks-with-interactive-code-70f0f17f46c6).\n\nNote the colouration as they will be useful when generating the layers.\n\n![img2](https://cdn-images-1.medium.com/max/1600/1*aRMefObpm7AMVOZYYiQAMQ.png)  \n\n"},{"metadata":{"_uuid":"8f5f689cce44e2351e254fe688c5d081bbde6ce8"},"cell_type":"markdown","source":"## Set up environment and utility functions"},{"metadata":{"trusted":true,"_uuid":"06ec5d31edcf0f02f97db5f4a26f7dd3f675b6dc"},"cell_type":"code","source":"tf.set_random_seed(6464)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"205d675c474d424c4d10e9adc54f0ce042ea147d"},"cell_type":"markdown","source":"### Activation Functions"},{"metadata":{"trusted":true,"_uuid":"d93202cfed55f1e041252ccfef7aaac4744fb998"},"cell_type":"code","source":"def tf_relu(x): \n    return tf.nn.relu(x)\n\ndef d_tf_relu(s): \n    return tf.cast(tf.greater(s,0),dtype=tf.float32)\n\ndef tf_softmax(x): \n    return tf.nn.softmax(x)\n\ndef np_sigmoid(x): \n    1/(1 + np.exp(-1 *x))","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5ac0520c84c33b41f6aa4f863fe0a80dbdea1cd5"},"cell_type":"markdown","source":"### Generate Classes for Convolution Layers"},{"metadata":{"trusted":true,"_uuid":"7971cf3d9475e08ac6719d97d632bcd08cdc38ca"},"cell_type":"code","source":"class conlayer_left():\n    \n    def __init__(self,ker,in_c,out_c):\n        self.w = tf.Variable(tf.random_normal([ker,ker,in_c,out_c],stddev=0.05))\n\n    def feedforward(self,input,stride=1,dilate=1):\n        self.input  = input\n        self.layer  = tf.nn.conv2d(input,self.w,strides = [1,stride,stride,1],padding='SAME')\n        self.layerA = tf_relu(self.layer)\n        return self.layerA\n\nclass conlayer_right():\n    \n    def __init__(self,ker,in_c,out_c):\n        self.w = tf.Variable(tf.random_normal([ker,ker,in_c,out_c],stddev=0.05))\n\n    def feedforward(self,input,stride=1,dilate=1,output=1):\n        self.input  = input\n\n        current_shape_size = input.shape\n\n        self.layer = tf.nn.conv2d_transpose(input,self.w,\n        output_shape=[batch_size] + [int(current_shape_size[1].value*2),int(current_shape_size[2].value*2),int(current_shape_size[3].value/2)],strides=[1,2,2,1],padding='SAME')\n        self.layerA = tf_relu(self.layer)\n        return self.layerA","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b3deec99efe380e629db3fd58f9e6f8d9b1d1948"},"cell_type":"markdown","source":"## Normalize the data"},{"metadata":{"trusted":true,"_uuid":"b31b2aa1379eed3a21239d69faa8844fa1bb1fe0","scrolled":true},"cell_type":"code","source":"train_images = (img_vol - img_vol.min()) / (img_vol.max() - img_vol.min())\ntrain_labels = (seg_vol - seg_vol.min()) / (seg_vol.max() - seg_vol.min())\n\nprint('Images', img_vol.shape, 'Segmentations', seg_vol.shape)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9293cf2df9bffcd5f6a0f55bd612bf17e833ff45"},"cell_type":"markdown","source":"## Hyper-Parameters"},{"metadata":{"trusted":true,"_uuid":"1e63f4dc0b523aba296db61f5979be0dc276398a"},"cell_type":"code","source":"num_epoch = 100\ninit_lr = 0.0001\nbatch_size = 2","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"9b28bca3b85450c459322eb993459e69b06c7502"},"cell_type":"markdown","source":"## Generate Layers"},{"metadata":{"_uuid":"0d7bfa884637aa860d6bbc75efbbb4592d0af20f"},"cell_type":"markdown","source":"### Left Side of th U-Net architecutre or Red in the above image"},{"metadata":{"trusted":true,"_uuid":"f7be77397a39980a2e6da0e70f3b0b13b42f3057"},"cell_type":"code","source":"l1_1 = conlayer_left(3,1,3)\nl1_2 = conlayer_left(3,3,3)\nl1_3 = conlayer_left(3,3,3)\n\nl2_1 = conlayer_left(3,3,6)\nl2_2 = conlayer_left(3,6,6)\nl2_3 = conlayer_left(3,6,6)\n\nl3_1 = conlayer_left(3,6,12)\nl3_2 = conlayer_left(3,12,12)\nl3_3 = conlayer_left(3,12,12)\n\nl4_1 = conlayer_left(3,12,24)\nl4_2 = conlayer_left(3,24,24)\nl4_3 = conlayer_left(3,24,24)\n\nl5_1 = conlayer_left(3,24,48)\nl5_2 = conlayer_left(3,48,48)\nl5_3 = conlayer_left(3,48,24)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"f4ee264aff89a82b7c4a49dc32ffbe37eeaf9f63"},"cell_type":"markdown","source":"### Right Side of the U-Net Architecture or Blue in the above image\n"},{"metadata":{"trusted":true,"_uuid":"8a85ec30728d2eb040ad33c907d50bdd83bf6836"},"cell_type":"code","source":"# right\nl6_1 = conlayer_right(3,24,48)\nl6_2 = conlayer_left(3,24,24)\nl6_3 = conlayer_left(3,24,12)\n\nl7_1 = conlayer_right(3,12,24)\nl7_2 = conlayer_left(3,12,12)\nl7_3 = conlayer_left(3,12,6)\n\nl8_1 = conlayer_right(3,6,12)\nl8_2 = conlayer_left(3,6,6)\nl8_3 = conlayer_left(3,6,3)\n\nl9_1 = conlayer_right(3,3,6)\nl9_2 = conlayer_left(3,3,3)\nl9_3 = conlayer_left(3,3,3)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"749d3ac1a539bcecd9706fb7312fa102e68e3d74"},"cell_type":"markdown","source":"### Combining layer or bottle neck layer or Green in the above image"},{"metadata":{"trusted":true,"_uuid":"d320e96db756a5f6e92ff94b0969717968e3b905"},"cell_type":"code","source":"l10_final = conlayer_left(3,3,1)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"fecda321794c956371b66f7bbc94d363942bcf00"},"cell_type":"markdown","source":"## Generate architecture"},{"metadata":{"trusted":true,"_uuid":"fb121ba6217f411681f00cce3593a36c64c76213"},"cell_type":"code","source":"x = tf.placeholder(shape=[None,512,512,1],dtype=tf.float32)\ny = tf.placeholder(shape=[None,512,512,1],dtype=tf.float32)\n\nlayer1_1 = l1_1.feedforward(x)\nlayer1_2 = l1_2.feedforward(layer1_1)\nlayer1_3 = l1_3.feedforward(layer1_2)\n\nlayer2_Input = tf.nn.max_pool(layer1_3,ksize=[1,2,2,1],strides=[1,2,2,1],padding='VALID')\nlayer2_1 = l2_1.feedforward(layer2_Input)\nlayer2_2 = l2_2.feedforward(layer2_1)\nlayer2_3 = l2_3.feedforward(layer2_2)\n\nlayer3_Input = tf.nn.max_pool(layer2_3,ksize=[1,2,2,1],strides=[1,2,2,1],padding='VALID')\nlayer3_1 = l3_1.feedforward(layer3_Input)\nlayer3_2 = l3_2.feedforward(layer3_1)\nlayer3_3 = l3_3.feedforward(layer3_2)\n\nlayer4_Input = tf.nn.max_pool(layer3_3,ksize=[1,2,2,1],strides=[1,2,2,1],padding='VALID')\nlayer4_1 = l4_1.feedforward(layer4_Input)\nlayer4_2 = l4_2.feedforward(layer4_1)\nlayer4_3 = l4_3.feedforward(layer4_2)\n\nlayer5_Input = tf.nn.max_pool(layer4_3,ksize=[1,2,2,1],strides=[1,2,2,1],padding='VALID')\nlayer5_1 = l5_1.feedforward(layer5_Input)\nlayer5_2 = l5_2.feedforward(layer5_1)\nlayer5_3 = l5_3.feedforward(layer5_2)\n\nlayer6_Input = tf.concat([layer5_3,layer5_Input],axis=3)\nlayer6_1 = l6_1.feedforward(layer6_Input)\nlayer6_2 = l6_2.feedforward(layer6_1)\nlayer6_3 = l6_3.feedforward(layer6_2)\n\nlayer7_Input = tf.concat([layer6_3,layer4_Input],axis=3)\nlayer7_1 = l7_1.feedforward(layer7_Input)\nlayer7_2 = l7_2.feedforward(layer7_1)\nlayer7_3 = l7_3.feedforward(layer7_2)\n\nlayer8_Input = tf.concat([layer7_3,layer3_Input],axis=3)\nlayer8_1 = l8_1.feedforward(layer8_Input)\nlayer8_2 = l8_2.feedforward(layer8_1)\nlayer8_3 = l8_3.feedforward(layer8_2)\n\nlayer9_Input = tf.concat([layer8_3,layer2_Input],axis=3)\nlayer9_1 = l9_1.feedforward(layer9_Input)\nlayer9_2 = l9_2.feedforward(layer9_1)\nlayer9_3 = l9_3.feedforward(layer9_2)\n\nlayer10 = l10_final.feedforward(layer9_3)\n\ncost = tf.reduce_mean(tf.square(layer10-y))\nauto_train = tf.train.AdamOptimizer(learning_rate=init_lr).minimize(cost)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"scrolled":false,"_uuid":"d0c48cb8162b6c176b4cc7a5c2397f79131d8902"},"cell_type":"code","source":"with tf.Session() as sess:\n    sess.run(tf.global_variables_initializer())\n\n    for iter in range(num_epoch):\n        \n        # train\n        for current_batch_index in range(0,len(train_images),batch_size):\n            current_batch = train_images[current_batch_index:current_batch_index+batch_size,:,:,:]\n            current_label = train_labels[current_batch_index:current_batch_index+batch_size,:,:,:]\n            sess_results = sess.run([cost,auto_train],feed_dict={x:current_batch,y:current_label})\n            print(' Iter: ', iter, \" Cost:  %.32f\"% sess_results[0],end='\\r')\n        print('\\n-----------------------')\n        train_images,train_labels = shuffle(train_images,train_labels)\n\n        if iter % 10 == 0:\n            test_example =   train_images[:2,:,:,:]\n            test_example_gt = train_labels[:2,:,:,:]\n            sess_results = sess.run([layer10],feed_dict={x:test_example})\n\n            sess_results = sess_results[0][0,:,:,:]\n            test_example = test_example[0,:,:,:]\n            test_example_gt = test_example_gt[0,:,:,:]\n\n            plt.figure()\n            plt.imshow(np.squeeze(test_example),cmap='gray')\n            plt.axis('off')\n            plt.title('Original Image')  \n            plt.savefig(str(iter)+\"a_Original_Image.png\")\n\n            plt.figure()\n            plt.imshow(np.squeeze(test_example_gt),cmap='gray')\n            plt.axis('off')\n            plt.title('Ground Truth Mask')  \n            plt.savefig(str(iter)+\"b_Original_Mask.png\")            \n\n            plt.figure()\n            plt.imshow(np.squeeze(sess_results),cmap='gray')\n            plt.axis('off')\n            plt.title('Generated Mask') \n            plt.savefig(str(iter)+\"c_Generated_Mask.png\")\n\n            plt.figure()\n            plt.imshow(np.multiply(np.squeeze(test_example),np.squeeze(test_example_gt)),cmap='gray')\n            plt.axis('off')\n            plt.title(\"Ground Truth Overlay\")   \n            plt.savefig(str(iter)+\"d_Original_Image_Overlay.png\")\n\n            plt.figure()\n            plt.axis('off')\n            plt.imshow(np.multiply(np.squeeze(test_example),np.squeeze(sess_results)),cmap='gray')\n            plt.title(\"Generated Overlay\")                 \n\n            plt.close('all')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"dc14eafdf51bb8d44c3d9a505fe6772c3b6fec9b"},"cell_type":"markdown","source":"Look at the ouputs to view the results of the U-Net alogirthm"}],"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}