{"cells":[{"metadata":{"_cell_guid":"a490d5fb-3b73-4b2d-b2a0-f9cb12111f19","_uuid":"187f98c82996a15047b2562a32208ce8c02106e4"},"cell_type":"markdown","source":"# Statement of the problem\nIn this notebook I briefly describe a naive method that attempts to separate merged(overlapping) nuclei masks using OpenCV convexity analysis. The separating method can be applied to a segmentation output after applying skimage `label` or methods of sort to separated connected clusters. The analysis bellow makes the assumption that nuclei are convex-shaped (This is not 100% correct) .\n\n1. Show some examples of overlapping nuclei\n2.  Describe the convexity based separation code\n3. Apply the correction to one of the benchmar submissions (It actually improved LB 0.274 -> LB 0.300)"},{"metadata":{"_cell_guid":"459c91fc-a7aa-4d62-8ea9-0669b24390b3","_uuid":"57f768e80c7854b68e1acd1f9373e3a7d8d674c1"},"cell_type":"markdown","source":"## 1. Show some examples\nLet us read some train images and their corresponding masks."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","collapsed":true,"trusted":true},"cell_type":"code","source":"import os\nimport cv2\nimport numpy as np\nfrom matplotlib import pyplot as plt\n\ntrain_dirs = os.listdir(\"../input/data-science-bowl-2018/stage1_train/\")\nM = {}\nI = {}\nfor file_id in train_dirs[:10]:\n    train_filename = \"../input/data-science-bowl-2018/stage1_train/\"+file_id+\"/images/\"+file_id+\".png\"\n    I[file_id] = cv2.imread(train_filename) \n    train_mask_dirs = \"../input/data-science-bowl-2018/stage1_train/\"+file_id+\"/masks/\"\n    train_mask_files = os.listdir(train_mask_dirs)\n    M[file_id] = []\n    for train_mask in train_mask_files:\n        mask = cv2.imread(train_mask_dirs + train_mask, cv2.IMREAD_GRAYSCALE)\n        M[file_id].append(mask)","execution_count":13,"outputs":[]},{"metadata":{"_cell_guid":"f1ede3f9-0d59-4789-80d3-9554ca9ae84d","_uuid":"0b0466c5c65f03eeee3458826d9ee6387ff54130"},"cell_type":"markdown","source":"Plot second image (`7f34dfccd1bc2e2466ee3d6f74ff05821a0e5404e9cf2c9568da26b59f7afda5`) and the `np.amax` of all corresponding masks. Ac"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"id = 1\nimage = I[train_dirs[id]]\nmask = np.amax(M[train_dirs[id]], axis=0)\nplt.imshow(image)\nplt.show()\nplt.imshow(mask)\nplt.show()\n","execution_count":14,"outputs":[]},{"metadata":{"_cell_guid":"78b51606-0f56-49ff-81ea-67895b8abad4","_uuid":"8d294fb07b4888c06e9a8dc70a849a503b90b392"},"cell_type":"markdown","source":"Use the skimage `label` to separate the masks into connected blobs."},{"metadata":{"_cell_guid":"79530e86-13bc-4af2-a8f7-cb7ba1f9130f","_uuid":"71bf46c2b53f81d5c95cdd3540472691860d960e","trusted":true},"cell_type":"code","source":"from skimage.morphology import label\nlabeled_mask = label(mask)\nprint ('Found ', len(np.unique(labeled_mask)), ' connected masks')\nB = []\nfor i in np.unique(labeled_mask):\n    if i == 0: # id = 0 is for background\n        continue\n    mask_i = (labeled_mask==i).astype(np.uint8)\n    B.append(mask_i)\nn_blobs = len(B)","execution_count":15,"outputs":[]},{"metadata":{"_cell_guid":"5a8f00c6-c79c-4169-9eb4-d96502b1ca84","_uuid":"ffbe2ff2eb1ff6fd1cf0f73e5e5920efcc48fad0","trusted":true},"cell_type":"code","source":"f, axarr = plt.subplots(6, 7)\nf.set_figwidth(20)\nf.set_figheight(20)\n\nfor i in range(6):\n    for j in range(7):\n        k = i*7 + j\n        if k < n_blobs:\n            axarr[i][j].set_title('mask:' + str(k))\n            axarr[i][j].imshow(B[k])     \nplt.show()","execution_count":16,"outputs":[]},{"metadata":{"_cell_guid":"d1367711-0ba5-447b-b2a8-8aa21f517ddf","_uuid":"64137cddb76832f06fdb5324c4cf918dbff02373"},"cell_type":"markdown","source":"I guess by now you have noticed cases like mask 1, mask 5, mask 12 and mask 13 that deviate from nice convex shapes. "},{"metadata":{"_uuid":"4ec09267dec662fb8033bd51debd398bc5f7dfd6"},"cell_type":"markdown","source":"## 2. Describe the separation idea\nThe idea is to operate on a blob-by-blob fashion and detect the points that deviate from the convex hull that is around the blob. More details [here](http://https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_imgproc/py_contours/py_contours_more_functions/py_contours_more_functions.html). The basic idea is described in the steps below:\n1. Given a segmentation mask\n2. Apply `skimage.morphology.label` or `cv2.connectedComponents` to split it into disjoint masks\n3. For each one disjoint mask\n    \n    If the disjoint mask is not convex enough\n    \n     a.  Identify convexity defects points\n     \n      b. Draw a line of background (0) color to connect every two nearest defect points"},{"metadata":{"trusted":true,"_uuid":"dd106ac44078fe484c13fa65ea0b0ba650790baf"},"cell_type":"code","source":"from IPython.display import Image\nImage(\"../input/separation-case-image/mask_separation.png\")","execution_count":17,"outputs":[]},{"metadata":{"_cell_guid":"fbd76be0-5d25-4eab-ab50-36380c93dd10","_uuid":"a145b515bd10cede611cff8a4d2609b838d50a58"},"cell_type":"markdown","source":"\n      \n### 2.1 Convex enough\nWe use `skimage.measure.regionprops` method to characterize each individual mask. More specifically we check the ratio of `prop.convex_area/prop.filled_area`. If it is bigger than a predefined threshold then the shape of the mask is non-convex. \n\n```\nfrom skimage.measure import regionprops\n props = regionprops(mask_i, cache=False)\n  prop = props[0]\n  if prop.convex_area/prop.filled_area > 1.1:\n      print ('no convex')\n```\n \n Let us apply the above filter to the 42 masks of the example above:"},{"metadata":{"_cell_guid":"3ff18e4e-ac56-4bbe-bb14-db617d152ce9","_uuid":"25fb4ac0afa8729d55e911f99b9803664e9e4aa2","trusted":true},"cell_type":"code","source":"from skimage.measure import regionprops\nf, axarr = plt.subplots(6, 7)\nf.set_figwidth(20)\nf.set_figheight(20)\n\nfor i in range(6):\n    for j in range(7):\n        k = i*7 + j\n        if k < n_blobs:\n            props = regionprops(B[k], cache=False)\n            prop = props[0]\n            if prop.convex_area/prop.filled_area > 1.1:\n                axarr[i][j].set_title('non convex mask:' + str(k))\n            else:\n                axarr[i][j].set_title('convex mask:' + str(k))\n                \n            axarr[i][j].imshow(B[k])     \nplt.show()","execution_count":39,"outputs":[]},{"metadata":{"_cell_guid":"decc693f-6c23-425d-b516-c321b9c4a292","_uuid":"3076bd00fd7acf8aac00e77f75e47a63057e9ca0"},"cell_type":"markdown","source":"I guess by now we have a pretty good indication of convexity. We can move to correction:\n\n### 2.2 Correcting con-conex masks\nThe idea is based on the example [here](http://https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_imgproc/py_contours/py_contours_more_functions/py_contours_more_functions.html). Our aim is to locate convexity defect points and connect them with a straight line (background color) hence separating them. The defect points are returned as as list of : ` [ start point, end point, farthest point, approximate distance to farthest point ]`. We can use the approximate distance to the farthest point as indication of \"badness\". Essentially we gather all defect points and filter them using this approximate distance. In the code below we filter the defect points according to:\n![](http://)\n\n1.  `if dd[i] > 1.0` :  We need then normalized  distance to be greater than 1\n2. `dd[i]/np.max(dd) > 0.2`:  We need the point to be at least 80% far compared to the farthest point.\n\nAfter gathering all convexity defect points we begin connecting them by simply connecting the two nearest each time. This is naively implemented using two nested loops:\n```\n for i in range(len(points)):\n    f1 = points[i]\n    p1 = tuple(contour[f1][0])\n    nearest = None\n    min_dist = np.inf\n    for j in range(len(points)):\n        ...\n         dist = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) \n         if dist < min_dist:\n             ...\n             \n    # Connect point p1 to its nearest one.\n    cv2.line(thresh,p1, nearest, [0, 0, 0], 2)\n```"},{"metadata":{"_cell_guid":"4e883013-be2a-437b-bd15-02bc1d2164fc","_uuid":"faa3c251e63f1fa4f47f3d5f78b716cc74a14d16","collapsed":true,"trusted":true},"cell_type":"code","source":"def split_mask_v1(mask):\n    thresh = mask.copy().astype(np.uint8)\n    im2, contours, hierarchy = cv2.findContours(thresh, 2, 1)\n    i = 0 \n    for contour in contours:\n        if  cv2.contourArea(contour) > 20:\n            hull = cv2.convexHull(contour, returnPoints = False)\n            defects = cv2.convexityDefects(contour, hull)\n            if defects is None:\n                continue\n            points = []\n            dd = []\n\n            #\n            # In this loop we gather all defect points \n            # so that they can be filtered later on.\n            for i in range(defects.shape[0]):\n                s,e,f,d = defects[i,0]\n                start = tuple(contour[s][0])\n                end = tuple(contour[e][0])\n                far = tuple(contour[f][0])\n                d = d / 256\n                dd.append(d)\n\n            for i in range(len(dd)):\n                s,e,f,d = defects[i,0]\n                start = tuple(contour[s][0])\n                end = tuple(contour[e][0])\n                far = tuple(contour[f][0])\n                if dd[i] > 1.0 and dd[i]/np.max(dd) > 0.2:\n                    points.append(f)\n\n            i = i + 1\n            if len(points) >= 2:\n                for i in range(len(points)):\n                    f1 = points[i]\n                    p1 = tuple(contour[f1][0])\n                    nearest = None\n                    min_dist = np.inf\n                    for j in range(len(points)):\n                        if i != j:\n                            f2 = points[j]                   \n                            p2 = tuple(contour[f2][0])\n                            dist = (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) \n                            if dist < min_dist:\n                                min_dist = dist\n                                nearest = p2\n\n                    cv2.line(thresh,p1, nearest, [0, 0, 0], 2)\n    return thresh     ","execution_count":40,"outputs":[]},{"metadata":{"_cell_guid":"947276a3-3400-4ef8-8f7c-b556e1cdb776","_uuid":"7521087561d5de030804300f3773704d746f5f10"},"cell_type":"markdown","source":"Let us see the separation in action."},{"metadata":{"_cell_guid":"75315bab-480d-435b-988f-133c71549025","_uuid":"d7e2f151c0c2de738aeb15acac260c1552abede5","trusted":true},"cell_type":"code","source":"b_split = split_mask_v1(B[1])\nf, axarr = plt.subplots(1, 2)\nf.set_figheight(8)\nf.set_figwidth(16)\n\naxarr[0].set_title('Overlapping masks')\naxarr[0].imshow(B[1], cmap='gray')\naxarr[1].set_title('Separated masks')\naxarr[1].imshow(b_split, cmap='gray')\nplt.show()","execution_count":41,"outputs":[]},{"metadata":{"_cell_guid":"5bafe567-ba0d-4ec4-84d7-93f155d0f13d","_uuid":"f7ef0dbe1f90acdea3faf42a2e3a3f8a79d2e8d4"},"cell_type":"markdown","source":"## 3. Modify public kernel [Pure image processing LB 0.274](https://www.kaggle.com/ahassaine/pure-image-processing-lb-0-274)\nMost of the following code is borrowed from the Ali Hassaïne' s public kernel. "},{"metadata":{"_cell_guid":"6b9a3450-04ca-4a95-b173-28ba67623bcd","_uuid":"67f58072e191837a170fde45f65cef9c209ff73e","collapsed":true,"trusted":true},"cell_type":"code","source":"test_dirs = os.listdir(\"../input/data-science-bowl-2018/stage1_test/\")\ntest_filenames=[\"../input/data-science-bowl-2018/stage1_test/\"+file_id+\"/images/\"+file_id+\".png\" for file_id in test_dirs]\ntest_images=[cv2.imread(imagefile) for imagefile in test_filenames]","execution_count":42,"outputs":[]},{"metadata":{"_cell_guid":"d0c93c82-cdcb-445c-8590-3cefa817a713","_uuid":"203321d20c043c44cf0b9bce90268bf161e0b52c","collapsed":true,"trusted":true},"cell_type":"code","source":"from skimage.measure import regionprops\ndef process(img_rgb):\n    #green channel happends to produce slightly better results\n    #than the grayscale image and other channels\n    img_gray=img_rgb[:,:,1]#cv2.cvtColor(img_rgb, cv2.COLOR_BGR2GRAY)\n    #morphological opening (size tuned on training data)\n    circle7=cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7))\n    img_open=cv2.morphologyEx(img_gray, cv2.MORPH_OPEN, circle7)\n    #Otsu thresholding\n    img_th=cv2.threshold(img_open,0,255,cv2.THRESH_OTSU)[1]\n    #Invert the image in case the objects of interest are in the dark side\n    if(np.sum(img_th==255)>np.sum(img_th==0)):\n        img_th=cv2.bitwise_not(img_th)\n    #second morphological opening (on binary image this time)\n    bin_open=cv2.morphologyEx(img_th, cv2.MORPH_OPEN, circle7) \n    #connected components\n    cc=cv2.connectedComponents(bin_open)[1]\n    #cc=segment_on_dt(bin_open,20)\n    return cc\n\ndef rle_encoding(cc):\n    values=list(np.unique(cc))\n    values.remove(0)\n    RLEs=[]\n    for v in values:\n        dots = np.where(cc.T.flatten() == v)[0]\n        run_lengths = []\n        prev = -2\n        for b in dots:\n            if (b>prev+1):\n                run_lengths.extend((b + 1, 0))\n            run_lengths[-1] += 1\n            prev = b\n        RLEs.append(run_lengths)\n    return RLEs\n\n","execution_count":43,"outputs":[]},{"metadata":{"_cell_guid":"a27f578f-473b-4a50-9d47-e0bea20f80c5","_uuid":"559ddcb4670a6b2a7a6ce2e52dc87f1a7452d73c"},"cell_type":"markdown","source":"Except this function which does apply the separation and filtering:\n\n1. Takes as input as mask created by `cv2.connectedComponents`. That means that background is set to 0 and mask pixels have values from 1, 2, ...\n2. Breaks the masks by filtering value: `` `mask_i = (mask==i).astype(np.uint8)`\n3. Calculate the `regionprops` properties: ` props = regionprops(mask_i, cache=False)` \n4. If `prop.convex_area/prop.filled_area > 1.1` proceeds on splitting the mask\n5. Appends the splitted mask into a list\n\nFinally all masks are grouped together once more into as single one  ` np.amax(masks, axis=0)` and we call `label` to separate them once more."},{"metadata":{"_cell_guid":"8c93e781-7001-4e8d-bced-087445b9fe31","_uuid":"c4e362095c35819969916d3fd8763b00e31d3209","collapsed":true,"trusted":true},"cell_type":"code","source":"def split_and_relabel(mask):\n    masks = []\n    for i in np.unique(mask):\n        if i == 0: # id = 0 is for background\n            continue\n        mask_i = (mask==i).astype(np.uint8)\n        props = regionprops(mask_i, cache=False)\n        if len(props) > 0:\n            prop = props[0]\n            if prop.convex_area/prop.filled_area > 1.1:\n                mask_i = split_mask_v1(mask_i)\n        masks.append(mask_i)\n        \n    masks = np.array(masks)\n    masks_combined = np.amax(masks, axis=0)\n    labels = label(masks_combined)\n    return labels","execution_count":44,"outputs":[]},{"metadata":{"_cell_guid":"432a4a0d-7f69-49a3-87ff-7d939d4d851b","_uuid":"117cd9f370400b903d1b2d1d0a20a47d5c1e39ba"},"cell_type":"markdown","source":"Process each test image and create initial masks."},{"metadata":{"_cell_guid":"a95c6ea6-c909-4676-8b1a-0f640a5a65fa","_uuid":"0685f21b0a53f43307837d03a82914831d86088b","collapsed":true,"trusted":true},"cell_type":"code","source":"test_connected_components=[process(img)  for img in test_images]","execution_count":45,"outputs":[]},{"metadata":{"_cell_guid":"a663a92b-6286-4ac4-8909-05476ca17ae9","_uuid":"59ba8095edb22a4fc37d39db5b8e9628fbe0e5c5"},"cell_type":"markdown","source":"Use the splitting methology to each case:"},{"metadata":{"_cell_guid":"dbea7241-09d4-4c57-b786-d5d1db8d84b0","_uuid":"780716d24705b3c8272a1c8f878cbd4d3c5389dd","collapsed":true,"trusted":true},"cell_type":"code","source":"test_connected_components_split=[split_and_relabel(img)  for img in test_connected_components]","execution_count":46,"outputs":[]},{"metadata":{"_cell_guid":"d2019eb5-2825-457b-99a0-ab2b077d78d7","_uuid":"cd9be5545cb1fd0b66b37e08d2542f7146edb1c1"},"cell_type":"markdown","source":"Let us see what happened to three cases:"},{"metadata":{"_cell_guid":"24250fa1-3393-4a41-99ec-4daa11466eca","_uuid":"990cc24c4fc1b251bc0c1a45d6929ecc617a22ed","trusted":true},"cell_type":"code","source":"f, axarr = plt.subplots(3, 2)\nf.set_figwidth(16)\nf.set_figheight(16)\n\nfor i in range(3):\n    axarr[i, 0].set_title('Original ')\n    axarr[i, 0].imshow(test_connected_components[i])\n    axarr[i, 1].set_title('Splitted masks ')\n    axarr[i, 1].imshow(test_connected_components_split[i])\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"ed22df3a-b7b5-4255-bcd0-957b7f192009","_uuid":"d4ffedcdc6f91a7287df35a50988551092578ae0"},"cell_type":"markdown","source":"By visually inspecting we can see some good cases and some over-splitted cases. There is no free lunch after all. However I am working on some extentions:\n\n* Connecting only nearest convexity defect points is obviously a poor way\n* For masks near the margins we may not find symmetric convexity defect point"},{"metadata":{"_cell_guid":"4dfceb26-d623-408e-a0cc-fbd55366f078","_uuid":"2c3de07acfa48e318676279b771be620549a2998","collapsed":true,"trusted":true},"cell_type":"code","source":"test_RLEs=[rle_encoding(cc) for cc in test_connected_components_split]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"9aafd55d-2593-4e78-b10b-b81f22fa1c5f","_uuid":"2117631a58d77d40481a976ce11b1c349b1a498c","collapsed":true,"trusted":true},"cell_type":"code","source":"with open(\"predictions.csv\", \"w\") as myfile:\n    myfile.write(\"ImageId,EncodedPixels\\n\")\n    for i,RLEs in enumerate(test_RLEs):\n        for RLE in RLEs:\n            myfile.write(test_dirs[i]+\",\"+\" \".join([str(i) for i in RLE])+\"\\n\")","execution_count":49,"outputs":[]},{"metadata":{"collapsed":true,"trusted":true,"_uuid":"b56c5a76fe2557aa0a9906bd8b30a8c286d7c5e2"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}