{"cells":[{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"#import modules\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport cv2\nfrom PIL import Image\nfrom collections import Counter\n\nimport os\nprint(os.listdir(\"../input\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#import training data\ntrain = pd.read_csv(\"../input/hpa-single-cell-image-classification/train.csv\")\nprint(train.head())\n\n#map of targets in a dictionary\nsubcell_locs = {\n0:  \"Nucleoplasm\", \n1:  \"Nuclear membrane\",   \n2:  \"Nucleoli\",   \n3:  \"Nucleoli fibrillar center\" ,  \n4:  \"Nuclear speckles\",\n5:  \"Nuclear bodies\",\n6:  \"Endoplasmic reticulum\",   \n7:  \"Golgi apparatus\",\n8:  \"Intermediate filaments\",\n9:  \"Actin filaments\", \n10: \"Microtubules\",\n11:  \"Mitotic spindle\",\n12:  \"Centrosome\",   \n13:  \"Plasma membrane\",\n14:  \"Mitochondria\",   \n15:  \"Aggresome\",\n16:  \"Cytosol\",   \n17:  \"Vesicles and punctate cytosolic patterns\",   \n18:  \"Negative\"\n}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Each image is a 4-channel image with the protein of interest in the green channel. It is the subcellular localization of this protein which is recorded in the Target column of the train.csv file. The red channel corresponds to microtubules, the blue channel to the nucleus and the yellow channel to the endoplasmid reticulum. Let's display the different channels of the image with ID == 1, since it contains several subcelullar locations for our protein of interest. Then we will overlay the green and yellow channel, as the yellow channel gives a good indication of the cell shape."},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"The image with ID == 1 has the following labels:\", train.loc[1, \"Label\"])\nprint(\"These labels correspond to:\")\nfor location in train.loc[1, \"Label\"].split('|')[0]:\n    print(\"-\", subcell_locs[int(location)])\n\n#reset seaborn style\nsns.reset_orig()\n\n#get image id\nim_id = train.loc[1, \"ID\"]\n\n#create custom color maps\ncdict1 = {'red':   ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0)),\n\n         'green': ((0.0,  0.0, 0.0),\n                   (0.75, 1.0, 1.0),\n                   (1.0,  1.0, 1.0)),\n\n         'blue':  ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0))}\n\ncdict2 = {'red':   ((0.0,  0.0, 0.0),\n                   (0.75, 1.0, 1.0),\n                   (1.0,  1.0, 1.0)),\n\n         'green': ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0)),\n\n         'blue':  ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0))}\n\ncdict3 = {'red':   ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0)),\n\n         'green': ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0)),\n\n         'blue':  ((0.0,  0.0, 0.0),\n                   (0.75, 1.0, 1.0),\n                   (1.0,  1.0, 1.0))}\n\ncdict4 = {'red': ((0.0,  0.0, 0.0),\n                   (0.75, 1.0, 1.0),\n                   (1.0,  1.0, 1.0)),\n\n         'green': ((0.0,  0.0, 0.0),\n                   (0.75, 1.0, 1.0),\n                   (1.0,  1.0, 1.0)),\n\n         'blue':  ((0.0,  0.0, 0.0),\n                   (1.0,  0.0, 0.0))}\n\nplt.register_cmap(name='greens', data=cdict1)\nplt.register_cmap(name='reds', data=cdict2)\nplt.register_cmap(name='blues', data=cdict3)\nplt.register_cmap(name='yellows', data=cdict4)\n\n#get each image channel as a greyscale image (second argument 0 in imread)\ngreen = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_green.png'.format(im_id), 0)\nred = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_red.png'.format(im_id), 0)\nblue = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_blue.png'.format(im_id), 0)\nyellow = cv2.imread('../input/hpa-single-cell-image-classification/train/{}_yellow.png'.format(im_id), 0)\n\n#display each channel separately\nfig, ax = plt.subplots(nrows = 2, ncols=2, figsize=(15, 15))\nax[0, 0].imshow(green, cmap=\"greens\")\nax[0, 0].set_title(\"Protein of interest\", fontsize=18)\nax[0, 1].imshow(red, cmap=\"reds\")\nax[0, 1].set_title(\"Microtubules\", fontsize=18)\nax[1, 0].imshow(blue, cmap=\"blues\")\nax[1, 0].set_title(\"Nucleus\", fontsize=18)\nax[1, 1].imshow(yellow, cmap=\"yellows\")\nax[1, 1].set_title(\"Endoplasmic reticulum\", fontsize=18)\nfor i in range(2):\n    for j in range(2):\n        ax[i, j].set_xticklabels([])\n        ax[i, j].set_yticklabels([])\n        ax[i, j].tick_params(left=False, bottom=False)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#stack nucleus and microtubules images\n#create blue nucleus and red microtubule images\nnuclei = cv2.merge((np.zeros((2048, 2048),dtype='uint8'), np.zeros((2048, 2048),dtype='uint8'), blue))\nmicrotub = cv2.merge((red, np.zeros((2048, 2048),dtype='uint8'), np.zeros((2048, 2048),dtype='uint8')))\n\n#create ROI\nrows, cols, _ = nuclei.shape\nroi = microtub[:rows, :cols]\n\n#create a mask of nuclei and invert mask\nnuclei_grey = cv2.cvtColor(nuclei, cv2.COLOR_BGR2GRAY)\nret, mask = cv2.threshold(nuclei_grey, 10, 255, cv2.THRESH_BINARY)\nmask_inv = cv2.bitwise_not(mask)\n\n#make area of nuclei in ROI black\nred_bg = cv2.bitwise_and(roi, roi, mask=mask_inv)\n\n#select only region with nuclei from blue\nblue_fg = cv2.bitwise_and(nuclei, nuclei, mask=mask)\n\n#put nuclei in ROI and modify red\ndst = cv2.add(red_bg, blue_fg)\nmicrotub[:rows, :cols] = dst\n\n#show result image\nfig, ax = plt.subplots(figsize=(8, 8))\nax.imshow(microtub)\nax.set_title(\"Nuclei (blue) + microtubules (red)\", fontsize=15)\nax.set_xticklabels([])\nax.set_yticklabels([])\nax.tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labels_num = [value.split('|') for value in train['Label']]\nlabels_num_flat = list(map(int, [item for sublist in labels_num for item in sublist]))\nlabels = [\"\" for _ in range(len(labels_num_flat))]\nfor i in range(len(labels_num_flat)):\n    labels[i] = subcell_locs[labels_num_flat[i]]\n\nfig, ax = plt.subplots(figsize=(15, 5))\npd.Series(labels).value_counts().plot(kind='bar', fontsize=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"According to Chen et al. 2007, if images are segmented into single cell regions, additional features that are not appropriate for whole fields can be calculated after seeded watershed segmentation. Nucleus images provide a means to identify each cell, so image segmentation may start by identification of nuclei in images. The function cv2.connectedComponents provides a simple and effective means to label nuclei in images. Conversely, as shown on the following notebook cell, identification of whole cells using cv2.connectedComponents is not as efficient, due to the less homogeneous signal in the yellow channel of the image."},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the nucleus image\nret, thresh = cv2.threshold(blue, 0, 255, cv2.THRESH_BINARY)\n#display threshold image\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(thresh, cmap=\"Greys\")\nax[0].set_title(\"Threshold\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise\nkernel = np.ones((5,5),np.uint8)\nopening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)\nax[1].imshow(opening, cmap=\"Greys\")\nax[1].set_title(\"Morphological opening\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n# Marker labelling\nret, markers = cv2.connectedComponents(opening)\n# Map component labels to hue val\nlabel_hue = np.uint8(179 * markers / np.max(markers))\nblank_ch = 255 * np.ones_like(label_hue)\nlabeled_img = cv2.merge([label_hue, blank_ch, blank_ch])\n# cvt to BGR for display\nlabeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img[label_hue==0] = 0\nax[2].imshow(labeled_img)\nax[2].set_title(\"Markers\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the endoplasmic reticulum image\nret, thresh = cv2.threshold(yellow, 4, 255, cv2.THRESH_BINARY)\n#display threshold image\nfig, ax = plt.subplots(ncols=4, figsize=(20, 20))\nax[0].imshow(thresh, cmap=\"Greys\")\nax[0].set_title(\"Threshold\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise\nkernel = np.ones((5,5),np.uint8)\nopening = cv2.morphologyEx(thresh, cv2.MORPH_OPEN, kernel)\nax[1].imshow(opening, cmap=\"Greys\")\nax[1].set_title(\"Morphological opening\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n#morphological closing\nclosing = cv2.morphologyEx(opening, cv2.MORPH_CLOSE, kernel)\nax[2].imshow(closing, cmap=\"Greys\")\nax[2].set_title(\"Morphological closing\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)\n\n# Marker labelling\nret, markers = cv2.connectedComponents(closing)\n# Map component labels to hue val\nlabel_hue = np.uint8(179 * markers / np.max(markers))\nblank_ch = 255 * np.ones_like(label_hue)\nlabeled_img = cv2.merge([label_hue, blank_ch, blank_ch])\n# cvt to BGR for display\nlabeled_img = cv2.cvtColor(labeled_img, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img[label_hue==0] = 0\nax[3].imshow(labeled_img)\nax[3].set_title(\"Markers\", fontsize=15)\nax[3].set_xticklabels([])\nax[3].set_yticklabels([])\nax[3].tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try different simple thresholding methods. Description of threshold types can be found here and here."},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply threshold on the endoplasmic reticulum image\nret, thresh1 = cv2.threshold(yellow, 4, 255, cv2.THRESH_BINARY)\nret, thresh2 = cv2.threshold(yellow, 4, 255, cv2.THRESH_TRUNC)\nret, thresh3 = cv2.threshold(yellow, 4, 255, cv2.THRESH_TOZERO)\n\n#display threshold images\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(thresh1, cmap=\"Greys\")\nax[0].set_title(\"Binary\", fontsize=15)\n\nax[1].imshow(thresh2, cmap=\"Greys\")\nax[1].set_title(\"Trunc\", fontsize=15)\n\nax[2].imshow(thresh3, cmap=\"Greys\")\nax[2].set_title(\"To zero\", fontsize=15)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"To zero simple thresholding is not adapted at all for identifying cell boundaries based on the yellow channel. Even after playing with the upper and lower parameter values, no satisfactory result is obtained. Binary and truncate methods work better. Let's see how connectedComponents work after both thresholding methods."},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, ax = plt.subplots(ncols=4, figsize=(20, 20))\n\n#morphological opening to remove noise after binary thresholding\nkernel = np.ones((5,5),np.uint8)\nopening1 = cv2.morphologyEx(thresh1, cv2.MORPH_OPEN, kernel)\nax[0].imshow(opening1, cmap=\"Greys\")\nax[0].set_title(\"Morphological opening (binary)\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n#morphological closing after binary thresholding\nclosing1 = cv2.morphologyEx(opening1, cv2.MORPH_CLOSE, kernel)\nax[1].imshow(closing1, cmap=\"Greys\")\nax[1].set_title(\"Morphological closing (binary)\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)\n\n#morphological opening to remove noise after truncate thresholding\nkernel = np.ones((5,5),np.uint8)\nopening2 = cv2.morphologyEx(thresh2, cv2.MORPH_OPEN, kernel)\nax[2].imshow(opening2, cmap=\"Greys\")\nax[2].set_title(\"Morphological opening (truncate)\", fontsize=15)\nax[2].set_xticklabels([])\nax[2].set_yticklabels([])\nax[2].tick_params(left=False, bottom=False)\n\n#morphological closing after truncate thresholding\nclosing2 = cv2.morphologyEx(opening2, cv2.MORPH_CLOSE, kernel)\nax[3].imshow(closing2, cmap=\"Greys\")\nax[3].set_title(\"Morphological closing (truncate)\", fontsize=15)\nax[3].set_xticklabels([])\nax[3].set_yticklabels([])\nax[3].tick_params(left=False, bottom=False)\n\nfig, ax = plt.subplots(ncols=2, figsize=(10, 10))\n# Marker labelling for binary thresholding\nret, markers1 = cv2.connectedComponents(closing1)\n# Map component labels to hue val\nlabel_hue1 = np.uint8(179 * markers1 / np.max(markers1))\nblank_ch1 = 255 * np.ones_like(label_hue1)\nlabeled_img1 = cv2.merge([label_hue1, blank_ch1, blank_ch1])\n# cvt to BGR for display\nlabeled_img1 = cv2.cvtColor(labeled_img1, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img1[label_hue1==0] = 0\nax[0].imshow(labeled_img1)\nax[0].set_title(\"Markers (binary)\", fontsize=15)\nax[0].set_xticklabels([])\nax[0].set_yticklabels([])\nax[0].tick_params(left=False, bottom=False)\n\n# Marker labelling for truncate thresholding\nret, markers2 = cv2.connectedComponents(closing2)\n# Map component labels to hue val\nlabel_hue2 = np.uint8(179 * markers2 / np.max(markers2))\nblank_ch2 = 255 * np.ones_like(label_hue2)\nlabeled_img2 = cv2.merge([label_hue2, blank_ch2, blank_ch2])\n# cvt to BGR for display\nlabeled_img2 = cv2.cvtColor(labeled_img2, cv2.COLOR_HSV2BGR)\n# set bg label to black\nlabeled_img2[label_hue2==0] = 0\nax[1].imshow(labeled_img2)\nax[1].set_title(\"Markers (truncate)\", fontsize=15)\nax[1].set_xticklabels([])\nax[1].set_yticklabels([])\nax[1].tick_params(left=False, bottom=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At this point it's not clear if truncate thresholding is an improvement compared to binary thresholding. Some cells are fused to each other while they should not be.\n\nOn the other hand. Adaptive thresholding methods apply a different threshold on different parts of the image, let's see how well it does on our images. See here for more explanations."},{"metadata":{"trusted":true},"cell_type":"code","source":"#apply adaptive threshold on endoplasmic reticulum image\ny_blur = cv2.medianBlur(yellow, 3)\n\n#apply adaptive thresholding\nret,th1 = cv2.threshold(y_blur, 5,255, cv2.THRESH_BINARY)\n\nth2 = cv2.adaptiveThreshold(y_blur, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 15, 3)\n\nth3 = cv2.adaptiveThreshold(y_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 15, 3)\n\n#display threshold images\nfig, ax = plt.subplots(ncols=3, figsize=(20, 20))\nax[0].imshow(th1, cmap=\"Greys\")\nax[0].set_title(\"Binary\", fontsize=15)\n\nax[1].imshow(th2, cmap=\"Greys_r\")\nax[1].set_title(\"Adaptive: mean\", fontsize=15)\n\nax[2].imshow(th3, cmap=\"Greys_r\")\nax[2].set_title(\"Adaptive: gaussian\", fontsize=15)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### I am totally refered to this great [kernel](https://www.kaggle.com/jschnab/exploring-the-human-protein-atlas-images) by [Jonathan Schnabel](https://www.kaggle.com/jschnab)\n\n# To be continued..."}],"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}