{"cells":[{"metadata":{},"cell_type":"markdown","source":"This notebook contains some basic EDA and notes about segmentation as it pertains to the HPA Single Cell Classification competition. I am a beginner in this field, and explaining my work as I go is as much for my benefit as anyone else's, so I appreciate any advice/corrections etc. "},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"!pip install selectivesearch\n!pip install \"../input/pycocotools/pycocotools-2.0-cp37-cp37m-linux_x86_64.whl\"\n!pip install \"../input/hpapytorchzoozip/pytorch_zoo-master\"\n!pip install \"../input/hpacellsegmentatorraman/HPA-Cell-Segmentation/\"","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt \nimport os\nimport cv2\nimport skimage.io as io\nimport skimage.segmentation\nimport selectivesearch\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatches\nimport scipy.ndimage as ndi\nfrom skimage import filters, measure, transform, util\nfrom skimage.morphology import (binary_erosion, closing, disk,\n                                remove_small_holes, remove_small_objects)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"**Defining the problem**  \n\nThis is a weakly supervised classification problem. The training labels are given at the image level, and we are tasked with predicting labels at the cell level. We will therefore have to identify each cell in each image, segment them, and classify them individually. \n\nThe labels themselves refer to the 'subcellur protein localization patterns of single cells'. We have 17 classes of these patterns, and a 'negative' class.\n\nThe challenge here will be segment the images correctly, and then find a way to use the imprecise labels to classify our segmented cells accurately. The problem can therefore be framed as an instance segmentation problem. \n\nFor each image in the test set, you must predict a list of instance segmentation masks and their associated detection score (Confidence). The submission csv file uses the following format:\n\nImageID,ImageWidth,ImageHeight,PredictionString  \nImageAID,ImageAWidth,ImageAHeight,LabelA1 ConfidenceA1 EncodedMaskA1 LabelA2 ConfidenceA2 EncodedMaskA2 ...  \nImageBID,ImageBWidth,ImageBHeight,LabelB1 ConfidenceB1 EncodedMaskB1 LabelB2 ConfidenceB2 EncodedMaskB2 …\n\n\nNote that a mask MAY have more than one class. If that is the case, predict separate detections for each class using the same mask.\n\n"},{"metadata":{},"cell_type":"markdown","source":"So we are submitting a list of masks, their respective labels, and the respective confidences with which these label predictions were made. \n\nThe submissions are evaluated using Mean Average Precision, or mAP. See the video below.   \nhttps://www.youtube.com/watch?v=FppOzcDvaDI"},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv('../input/hpa-single-cell-image-classification/train.csv')\nsub = pd.read_csv('../input/hpa-single-cell-image-classification/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sub.loc[0]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Labels are in the form shown above, will change it to a list of ints"},{"metadata":{"trusted":true},"cell_type":"code","source":"train.Label = train.Label.apply(lambda x: list(map(int,x.split('|'))))\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"a = sns.countplot([len(x) for x in train.Label])\na.set_title('number of labels per image')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Vast majority of images either have 1 2 or 3 labels, with at least one image with 5 labels. "},{"metadata":{"trusted":true},"cell_type":"code","source":"set(np.sum(train.Label))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"19 different training labels, want to look at their distributions in the training set. \n\nWhich are more common, do they tend to occur together?"},{"metadata":{"trusted":true},"cell_type":"code","source":"label_mat = np.zeros((len(train),19), int)\nfor image in range(len(train)):\n    for label in train.Label[image]:\n        label_mat[image, label] = 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig = plt.figure(figsize = (12,6))\na = sns.barplot(list(range(19)),label_mat.sum(axis = 0),)\na.set_title('Distribution of Label Occurances')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Quite imbalanced labels, 0 is by far the most common and 11, 18  are close to 0. Oversampling may be necessary. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize = (12,8))\na = sns.heatmap(pd.DataFrame(label_mat).corr())\na.set_title('Do labels Occur Together?')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.DataFrame(label_mat).corr()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Can see most labels are weakly correlated with each other, with some, for example 4 and 0, with strong negative correlation. \n\nNot many strong positive correlations. "},{"metadata":{},"cell_type":"markdown","source":"Each entry in train.csv has four associate images\n\nBlue - Nucleus  \nRed - Microtubules  \nYellow - Endoplasmic reticulum  \nGreen - Protein of interest\n\nFirst lets visualise a few of these separately, then find a way to combine them "},{"metadata":{"trusted":true},"cell_type":"code","source":"colours = ['_red.png', '_blue.png', '_yellow.png', '_green.png']\nTRAIN_PATHS = '../input/hpa-single-cell-image-classification/train'\ntrain_paths = [[os.path.join(TRAIN_PATHS, train.iloc[idx,0])+ colour for colour in colours] for idx in range(len(train))]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"titles = ['microtubules','nucleus', 'endoplasmic reticulum', 'protein of interest']\nfig, axs = plt.subplots(3, 4, figsize =(16,8))\nfor entry in range(3):\n    for channel in range(4):\n        img = plt.imread(train_paths[entry][channel])\n        axs[entry, channel].imshow(img)        \n        if entry == 0:\n            axs[0, channel].set_title(titles[channel])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def to_rgb(idx,paths = train_paths, gbr = False, blue_only = False):\n    red = cv2.imread(paths[idx][0],0)\n    yellow = cv2.imread(paths[idx][2],0)\n    blue = cv2.imread(paths[idx][1],0)\n    \n    if blue_only: \n        return cv2.resize(blue, (512,512))\n    else:\n        return np.dstack((cv2.resize(red,(512,512)),cv2.resize(yellow,(512,512)) ,cv2.resize(blue,(512, 512))))\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#to visualise classes, will take images with single class labels. \n\nsingle_label = []\nfor label in range(19):\n    for idx in range(len(train)):\n        if train.loc[idx, 'Label'] == [label]:\n            single_label.append(idx)\n            break\ntitles_2 = [\n'0-Nucleoplasm',\n'1-Nuclear membrane',\n'2-Nucleoli',\n'3-Nucleoli fibrillar center',\n'4-Nuclear speckles',\n'5-Nuclear bodies',\n'6-Endoplasmic reticulum',\n'7-Golgi apparatus',\n'8-Intermediate filaments',\n'9-Actin filaments',\n'10-Microtubules',\n'11-Mitotic spindle',\n'12-Centrosome',\n'13-Plasma membrane',\n'14-Mitochondria',\n'15-Aggresome',\n'16-Cytosol',\n'17-Vesicles and punctate cytosolic patterns',\n'18-Negative'\n]\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(19,1, figsize = (10, 100))\n\nfor label, idx in enumerate(single_label):\n    axs[label].imshow(to_rgb(idx))\n    axs[label].set_title(titles_2[label])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A note about 'Negative' labels from the organiser\n\n> 1) Image-level label uncertainty The image-level labels are what we refer to as weak or noisy. During annotation, the image-level labels are set per sample (i.e per a group of up to 6 images from the same sample). This means that the labels present in the majority of the images will be annotated. For negative samples it is not uncommon that lets say 4 images show no staining, while the remaining 2 show some unspecific staining or some granular pattern. If you compare the image-level label with the precise pattern observed in any given cell from this group of images, the label will be correct for the vast majority of cells, but perhaps not for all of them (as in your example 1 and 3 above).\n\n> 2) Single cell label accuracy in test set The test set consists of images where each single cell has been annotated independently. Hence the accuracy of these labels is much better, and will be correct for each cell in every image. The statement below made by Trang Le, @lnhtrang, in our notebook explaining the patterns is correct for how the test set was annotated.\n\n\nSo in other words, in some of the 'Negative' samples, one or two of the cells in the image may in fact have a protein pattern, but most will not. "},{"metadata":{},"cell_type":"markdown","source":"How to best carry out segmentation on this dataset? \n\nPopular models being used for this competition include\n\nhttps://github.com/CellProfiling/HPA-Cell-Segmentation  \n\nhttps://www.cellpose.org/\n\nThese are both pre-trained models based on variations of convolutional neural networks called U-Nets. Before going near either of those, as a newcomer to the field of computer vision and object identification, I want to explore some more basic segmentation techniques in order to appreciate the need for the tools listed above. \n\n\nThe task we have is to perform single cell classification, therefore we are going to have to produce masks over each cell we identify and classify them individually. How best to produce these masks? Why do we need a pre-trained CNN at all?\n"},{"metadata":{},"cell_type":"markdown","source":"Thresholding pretty clearly isn't going to be the answer given that the task is to segment the individual cells in the image, and the background is already completely distinct."},{"metadata":{},"cell_type":"markdown","source":"The first segmentation technique that I will implement is a graph-based technique first proposed by Felsenzwalb. It involves representing our image as an undirected graph G = (V,E) of vertices and edges, with each vertex being a single pixel and each edge connecting a pair of two vertices. An edge has a 'weight' which is given by the distance between the two vertices (pixels) in the edge. The distance metric is the difference in colour, intensity, and location. \n\nA segmentation solution is a partition of V into multiple connected components. Starting with each pixel in its own component, the segmentation is arrived at through a 'bottom-up' technique. \n\nThe basic idea is that a if the distance metric between two components is small compared to the internal difference of both those components, then merge the two components, otherwise do nothing. \n\nThe higher the scale factor 'k', the higher the threshold for a decision to label two components to be separate. So higher k will give us larger components/segments.\n\nFor more information on how the algorithm works see http://cs.brown.edu/people/pfelzens/papers/seg-ijcv.pdf\n"},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(4,4, figsize = (30, 30))\nims = [to_rgb(i) for i in range(4)]\nk = [0,500, 1000, 2000]\nfor i, image in enumerate(ims):\n    for j in range(4):\n        if j == 0:\n            axs[i,0].imshow(image)\n        else:\n            axs[i,j].imshow(skimage.segmentation.felzenszwalb(image, scale = k[j]))\n            axs[i,j].set_title('Image {} k='.format(i) + str(k[j]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, this approach runs into problems where cells overlap, as is the case with much of our dataset. "},{"metadata":{},"cell_type":"markdown","source":"Selective Search is an algorithm used to propose regions which contain objects. It is built on top of a segmentation algorithm such as Felzenszwalb. \n\nIt is a hierarchical grouping algorithm, iteratively grouping together the most similar regions (which were obtained by Felzenszwalb), until a stopping point determined by the min_size parameter. \n\nsigma is a Gaussian blur parameter to smooth the image.  \n\n\nNote, while this algorithm returns bounding boxes, not the image masks we will need for submission,it will be interesting to see if it can reliably identify individual cells, especially when their boundaries overlap. "},{"metadata":{"trusted":true},"cell_type":"code","source":"\ndef selective_search_regions(idx, scale = 500, sigma = 0.8, min_size = 100, min_region = 1000):\n    img = to_rgb(idx)\n    img_lbl,regions = selectivesearch.selective_search(img, scale=scale, sigma=sigma, min_size=min_size)\n    candidates = set()\n    for r in regions:\n        # excluding same rectangle (with different segments)\n        if r['rect'] in candidates:\n            continue\n        # excluding regions smaller than 2000 pixels\n        if r['size'] < min_region:\n            continue\n        # distorted rects\n        x, y, w, h = r['rect']\n        if w / h > 1.2 or h / w > 1.2:\n            continue\n        candidates.add(r['rect'])\n    return img_lbl, candidates\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(4,1,figsize = (10,40))\nfor idx in range(4):\n    img = to_rgb(idx)\n    img_lbl, candidates = selective_search_regions(idx)\n    # draw rectangles on the original image\n    axs[idx].imshow(img)\n    for x, y, w, h in candidates:\n        rect = mpatches.Rectangle(\n            (x, y), w, h, fill=False, edgecolor='red', linewidth=1)\n        axs[idx].add_patch(rect)\n\n    \n    \n    ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can see, the results are pretty unsatisfactory. Even after playing around with the parameters, the set of parameters which correctly identify the different cells in one image will result in far too many regions (or too few) in the next. We need a method more robust than this for our dataset of thousands of images of cells. As well, we need a model that predicts the masks as well as just identifying the objects. "},{"metadata":{},"cell_type":"markdown","source":"A U-Net is a form of CNN designed for semantic segmentation. Because our task is closer to instance segmentation, there is quite a lot of post-processing to be done on the output of a U-Net in order to produce masks for the individual cells in our images. Luckily, a lot of that work has been done for us by the creators of https://github.com/CellProfiling/HPA-Cell-Segmentation. \n\nAnother approach would be to use a purpose-built instance segmentation model such as Mask R-CNN, which outputs individual object masks, bounding boxes, and class predictions for each objects.  "},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_images(df = train, paths = train_paths, n = 4):\n    blue = []\n    rgb = []\n    for i in range(4):\n        rgb_img = to_rgb(i)/255 \n        blue_img = to_rgb(i, blue_only = True)/255\n        rgb.append(rgb_img)\n        blue.append(blue_img)\n    return rgb, blue\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"import hpacellseg.cellsegmentator as cellsegmentator\nfrom hpacellseg.utils import label_cell, label_nuclei\nNUC_MODEL = \"../input/hpacellsegmentatormodelweights/dpn_unet_nuclei_v1.pth\"\nCELL_MODEL = \"../input/hpacellsegmentatormodelweights/dpn_unet_cell_3ch_v1.pth\"\nsegmentator = cellsegmentator.CellSegmentator(\n    NUC_MODEL,\n    CELL_MODEL,\n    device=\"cuda\",\n    multi_channel_model=True,\n)\n\nN = 4\nrgb, blue = load_images(n = N)\nmasks = []\n\n\nnuc_segmentations = segmentator.pred_nuclei(blue)\ncell_segmentations = segmentator.pred_cells(rgb, precombined = True)\n\nfor nuc, cell in zip(nuc_segmentations, cell_segmentations):\n    masks.append(label_cell(nuc, cell))\n\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(cell_segmentations[0])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The above shows the output from the semantic segementation performed by U-Net. Below will be the masks once post-processing is complete. "},{"metadata":{"trusted":true},"cell_type":"code","source":"fig, axs = plt.subplots(4,3, figsize = (20,30))\n\nfor i in range(N):\n    axs[i,0].imshow(to_rgb(i))\n    axs[i,1].imshow(masks[i][1])\n    axs[i,2].imshow(to_rgb(i))\n    axs[i,2].imshow(masks[i][1], alpha = 0.5)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"TBC when I have the time. Appreciate any comments + advice"}],"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}