{"cells":[{"metadata":{"_uuid":"8312fe142cec2e688788d4206a72c3d5ce157998","_cell_guid":"19ca1813-7f70-42f4-bce7-2a32f83887f9"},"cell_type":"markdown","source":"**Overlapping objects separation method**\n\nHi Kagglers 🍺\n\nI wan't to share my solution for separation of overlapping objects, by using intuition and simple linear algebra I have construted this method. It's not perfect but it can help you with instance separation in some cases of overlapped round objects.\n\nHow does it work:\n1.  round objects doesn't have inset corners, if corners are present you can asume that this is overlapping scenario.\n2. extract inset corners by substract convex hull from object.\n3. find all combinations of corners for separation candidate (some optimizations like crop, boundaries are used just for speed up).\n4. separation line is valid if parallel line with offset 2 is shorter than offset 5 and if offset 5 is wider than separation line. This is checked in both directions.\n\nMy solution is classic / minimal UNet with this post processing, currently at **0.382 LB**.\n\nDisclaimer: This is not a AI!"},{"metadata":{"_uuid":"5af1ff5abdf1f642a4c64b0f697b77a813f484bf","_cell_guid":"2a28ac5f-e066-4e28-9ec4-d76d5747d4a1","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":1,"outputs":[]},{"metadata":{"_uuid":"a4101399089a2d9c2354dea4e15acf17b4e07c38","_cell_guid":"38eb5559-09b7-4c97-8540-bed34d225bc8","trusted":true},"cell_type":"code","source":"import skimage.io\nimport numpy as np\nfrom skimage.morphology import label\nfrom skimage.color import label2rgb\nfrom skimage.io import imshow\nimport matplotlib.pyplot as plt\n\n# we dont need and dont like warnings, this is dangerous world we are living in\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimage_id = '1b6044e4858a9b7cee9b0028d8e54fbc8fb72e6c4424ab5b9f3859bfc72b33c5'\nimg = skimage.io.imread('../input/stage1_train/{}/images/{}.png'.format(image_id, image_id))[:, :, :3]\n\n# combine masks as one\nmasks = skimage.io.imread_collection('../input/stage1_train/{}/masks/*.png'.format(image_id)).concatenate()\nmask = np.zeros((img.shape[0], img.shape[1]), dtype=np.uint8)\nfor maskItem_ in masks:\n    mask = np.maximum(mask, (maskItem_ == 255))\n    \nprint(\"Source image\")\nimshow(img)\nplt.show()\n\nprint(\"Combined mask image with overlapping objects\")\nlabels = label(mask, neighbors=4)\nimshow(label2rgb(labels, bg_label=0))\nplt.show()","execution_count":2,"outputs":[]},{"metadata":{"_uuid":"0f87d1fa61e226e0e7b0d663d04532ef2d673f87","_cell_guid":"3c28069c-a919-402a-9ed6-488b3a53de43"},"cell_type":"markdown","source":"**There are some overlapping nuclei, let's try to separate them.**\n\nBunch of \"IF's\" and simple algebra ⏬"},{"metadata":{"_uuid":"a15dc7d9f627d507f4723d49f7b49135e779e81f","collapsed":true,"_cell_guid":"51193f3f-3aa4-46ae-998c-ee87ca545762","trusted":true},"cell_type":"code","source":"from skimage.measure import regionprops\nfrom skimage.morphology import convex_hull_image\nfrom skimage.segmentation import find_boundaries\nimport math\nfrom itertools import combinations\n\ndef extendLineToMask(y1,x1,y2,x2,mask):\n    if (y1 < 0) or (y2 < 0) or (x1 < 0) or (x2 < 0): # input validation\n        return 0,0,0,0\n    \n    tc = (np.array([y1, x1])+np.array([y2, x2]))/2.0 # get center point\n    \n    # extend line to image bound\n    if (x2 - x1 == 0): \n        # is vertical line\n        extendedLineY1 = 0\n        extendedLineX1 = x1\n        extendedLineY2 = mask.shape[0]-1\n        extendedLineX2 = x2\n    else: \n        # not vertical\n        # calculate slope\n        # use skimage line for better slope calculation\n        rrP1, ccP1 = skimage.draw.line(int(y1), int(x1), int(y2), int(x2))      \n        m = ((ccP1*rrP1).mean() - ccP1.mean()*rrP1.mean()) / ((ccP1**2).mean() - (ccP1.mean())**2)\n        if (m == 0.0):\n            return 0,0,0,0\n        # calculate b\n        b = y1-(m*x1)\n        \n        extendedLineX1 = 0\n        extendedLineY1 = m * extendedLineX1 + b # y = mx+b\n        #out of bound handling\n        if (extendedLineY1 < 0) or (extendedLineY1 > mask.shape[0]):\n            # get min/max Y\n            extendedLineY1 = max(min(mask.shape[0]-1, extendedLineY1), 0.0)\n            # recalculate X\n            extendedLineX1 = (extendedLineY1-b) / m\n                    \n        extendedLineX2 = mask.shape[1]-1\n        extendedLineY2 = m * extendedLineX2 + b\n        #out of bound handling\n        if (extendedLineY2 < 0) or (extendedLineY2 > mask.shape[0]):\n            # get min/max Y\n            extendedLineY2 = max(min(mask.shape[0]-1, extendedLineY2), 0.0)\n            # recalculate X\n            extendedLineX2 = (extendedLineY2-b) / m\n     # check infinity\n    if (math.isinf(extendedLineX1) or math.isinf(extendedLineX2) or math.isinf(extendedLineY1) or math.isinf(extendedLineY2)):\n        return 0,0,0,0\n\n    # get extended line\n    rrP1, ccP1 = skimage.draw.line(int(extendedLineY1), int(extendedLineX1), int(extendedLineY2), int(extendedLineX2))\n\n    # get index of center point\n    linecenterir = np.nonzero(rrP1 == int(tc[0]))[0]\n    linecenteric = np.nonzero(ccP1 == int(tc[1]))[0]\n    if (len(linecenteric) < len(linecenterir)):\n        linecenteri = linecenteric\n    else:\n        linecenteri = linecenterir\n    if (len(linecenteri)) > 0:\n        # trim line to mask\n        linecenterindex = linecenteri[0]\n        try: # there are some situations where issue occurs (line outside image) but I haven't had the time to add validation\n            # get left and right part of the line (this is not actually direction of line but indexation of array)\n            lineontheleft = mask[rrP1[0:linecenterindex],ccP1[0:linecenterindex]]\n            lineontheright = mask[rrP1[linecenterindex:],ccP1[linecenterindex:]]\n            # trim to mask\n            lineonthelefttrimi = np.nonzero(lineontheleft == False)[0]\n            lineontherighttrimi = np.nonzero(lineontheright == False)[0]\n            # out of image bounds > better way will be to pad image with 1px False border, but this is faster\n            if (len(lineonthelefttrimi) < 1) and (mask[rrP1[0],ccP1[0]] == True):\n                lineonthelefttrimi = [-1]\n            if (len(lineontherighttrimi) < 1) and (mask[rrP1[-1],ccP1[-1]] == True):\n                lineontherighttrimi = [len(rrP1)]\n            if (len(lineonthelefttrimi) > 0) and (len(lineontherighttrimi) > 0):\n                lineonthelefttrim = lineonthelefttrimi[-1]+1\n                lineontherighttrim = linecenterindex+lineontherighttrimi[0]\n                if (lineontherighttrim > lineonthelefttrim):\n                    extendedLineY1 = rrP1[lineonthelefttrim:lineontherighttrim][0]\n                    extendedLineX1 = ccP1[lineonthelefttrim:lineontherighttrim][0]\n                    extendedLineY2 = rrP1[lineonthelefttrim:lineontherighttrim][-1]\n                    extendedLineX2 = ccP1[lineonthelefttrim:lineontherighttrim][-1]\n                else:\n                    extendedLineY1 = 0\n                    extendedLineX1 = 0\n                    extendedLineY2 = 0\n                    extendedLineX2 = 0\n            else:\n                # invalid line\n                extendedLineY1 = extendedLineX1 = extendedLineY2 = extendedLineX2 = 0\n        except:\n            # invalid line\n            extendedLineY1 = extendedLineX1 = extendedLineY2 = extendedLineX2 = 0\n    else:\n        # invalid line\n        extendedLineY1 = extendedLineX1 = extendedLineY2 = extendedLineX2 = 0\n\n    # return line coordinates\n    return int(extendedLineY1), int(extendedLineX1), int(extendedLineY2), int(extendedLineX2)\n\n# parallel line helper function\ndef parallelLine(px1, px2, offsetPixels, length=0.0):\n    if (length == 0.0):\n        length = math.sqrt((px1[1]-px2[1])*(px1[1]-px2[1])+(px1[0]-px2[0])*(px1[0]-px2[0]))\n    x1p = px1[1] + offsetPixels * (px2[0]-px1[0]) / length\n    x2p = px2[1] + offsetPixels * (px2[0]-px1[0]) / length\n    y1p = px1[0] + offsetPixels * (px1[1]-px2[1]) / length\n    y2p = px2[0] + offsetPixels * (px1[1]-px2[1]) / length\n    return [y1p,x1p], [y2p,x2p]\n\ndef splitValidation(px1, px2, img):\n    # calculate line distance\n    delta = (px1[1]-px2[1])*(px1[1]-px2[1])+(px1[0]-px2[0])*(px1[0]-px2[0])\n    if (delta > 0):\n        L = math.sqrt(delta)\n    else:\n        L = 0\n    # get top parallel line\n    plpx1a, plpx2a = parallelLine(px1, px2, -4.0, length=L)\n    # extend line to bound\n    y1pea, x1pea, y2pea, x2pea = extendLineToMask(plpx1a[0], plpx1a[1], plpx2a[0], plpx2a[1], img)\n    # get line length\n    delta = (x1pea-x2pea)*(x1pea-x2pea)+(y1pea-y2pea)*(y1pea-y2pea)\n    if (delta > 0):\n        La = math.sqrt(delta)\n    else:\n        La = 0\n    \n    # get bottom parallel line\n    plpx1b, plpx2b = parallelLine(px1, px2, 4.0, length=L)\n    # extend line to bound\n    y1peb, x1peb, y2peb, x2peb = extendLineToMask(plpx1b[0], plpx1b[1], plpx2b[0], plpx2b[1], img)\n    # if top and bottom line are longer than split line\n    delta = (x1peb-x2peb)*(x1peb-x2peb)+(y1peb-y2peb)*(y1peb-y2peb)\n    if (delta > 0):\n        Lb = math.sqrt(delta)\n    else:\n        Lb = 0\n\n    plpx1a2, plpx2a2 = parallelLine(px1, px2, -2.0, length=L)\n    # extend line to bound\n    y1pea2, x1pea2, y2pea2, x2pea2 = extendLineToMask(plpx1a2[0], plpx1a2[1], plpx2a2[0], plpx2a2[1], img)\n    # get line length\n    delta = (x1pea2-x2pea2)*(x1pea2-x2pea2)+(y1pea2-y2pea2)*(y1pea2-y2pea2)\n    if (delta > 0):\n        La2 = math.sqrt(delta)\n    else:\n        La2 = 0\n    \n    # get bottom parallel line\n    plpx1b2, plpx2b2 = parallelLine(px1, px2, 2.0, length=L)\n    # extend line to bound\n    y1peb2, x1peb2, y2peb2, x2peb2 = extendLineToMask(plpx1b2[0], plpx1b2[1], plpx2b2[0], plpx2b2[1], img)\n    # if top and bottom line are longer than split line\n    delta = (x1peb2-x2peb2)*(x1peb2-x2peb2)+(y1peb2-y2peb2)*(y1peb-y2peb2)\n    if (delta > 0):\n        Lb2 = math.sqrt(delta)\n    else:\n        Lb2 = 0\n\n    #if 5.0 line is bigger than split line and 2.0 line is bigger than 5.0 line\n    #if -5.0 line is bigger than split line and -2.0 line is bigger than -5.0 line\n    return (La >= L-2) and (La >= La2) and (Lb >= L-2) and (Lb >= Lb2)\n\ndef getCropMaskDimensions(img):\n    cropmatrix = np.transpose(img.nonzero())\n    cropcoordmin = np.amin(cropmatrix, axis=0)\n    cropcoordmax = np.amax(cropmatrix, axis=0)\n    cropy = cropcoordmin[0]\n    cropx = cropcoordmin[1]\n    cropy2 = cropcoordmax[0]+1\n    cropx2 = cropcoordmax[1]+1\n    return cropy, cropx, cropy2, cropx2\n    \ndef separation(img):\n    # make a copy\n    inputimg_ = np.copy(img)\n    \n    # crop image to object bound > this helps speed up the processing\n    cropy, cropx, cropy2, cropx2 = getCropMaskDimensions((inputimg_ == True))\n    img_ = inputimg_[cropy:cropy2,cropx:cropx2]\n\n    # define minimal convex area\n    MIN_CONVEX_AREA = 5\n    \n    # calculate convex hull of object\n    convexhull = convex_hull_image(img_)\n    # invert convex hull\n    convexhulldiff = img_^convexhull\n    # get boundaries > will help speed up processing\n    boundaries = find_boundaries(convexhulldiff, connectivity=1, mode='inner')\n    # split inset border objects\n    label_img = label(boundaries, neighbors=8)\n    if (label_img.max() > 1):\n        # calculate region props\n        regions = regionprops(label_img)\n        # get candidates\n        objectcandidates = []\n        for props in regions:\n            # limit to min convex area size\n            if props.convex_area > MIN_CONVEX_AREA:\n                y0, x0 = props.centroid # this is not used\n                orientation = props.orientation # this is not used\n                #convert image to points > will speed up processing\n                objimg = label_img == props.label\n                points_ = []\n                for row_ in range(0,objimg.shape[0]):\n                    for col_ in range(0,objimg.shape[1]):\n                        if (objimg[row_][col_]):\n                            points_.append([row_, col_])\n\n                objectcandidates.append([orientation, points_])\n\n        #get potential pairs\n        for c in combinations(range(0,len(objectcandidates)), 2):\n            # object A\n            iobjA = c[0]\n            # object B\n            iobjB = c[1]\n            # x,y border points for objects\n            xy1 = np.array(objectcandidates[iobjA][1])\n            xy2 = np.array(objectcandidates[iobjB][1])\n            # numpy optimized way to find nearest points between two objects\n            P = np.add.outer(np.sum(xy1**2, axis=1), np.sum(xy2**2, axis=1))\n            N = np.dot(xy1, xy2.T)\n            dists = np.sqrt(P - 2*N)\n\n            # find minimal distance pixel\n            distsm_ = np.argmin(dists)\n            # convert minimal pixel to x,y grid\n            distsmp_ = divmod(distsm_, dists.shape[1])\n            # convert to absolute position\n            px1 = objectcandidates[iobjA][1][distsmp_[0]]\n            px2 = objectcandidates[iobjB][1][distsmp_[1]]\n            # separation line\n            rr, cc = skimage.draw.line(px1[0], px1[1], px2[0], px2[1])\n            line_ = img_[rr, cc]\n            # check if line is valid = first check\n            if (np.count_nonzero(line_[1:-1] == False) == 0):\n                # check if split is valid = second check\n                if (splitValidation(px1,px2,img_)):\n                    img_[rr, cc] = 0\n    # reconstruct image with cropped area\n    inputimg_[cropy:cropy2,cropx:cropx2] = img_\n    return inputimg_","execution_count":3,"outputs":[]},{"metadata":{"_uuid":"705ed5634d12d0cf77513992bf2e7a7a4aaa0971","_cell_guid":"6246acb0-a9cd-49f7-a1b3-c1bdeedd762a"},"cell_type":"markdown","source":"**Let's try it out**\n\nMaybe we can separate grapes with it 🍇"},{"metadata":{"_uuid":"5751746f05c2fe5e34e20e944d0e0b7fb6c1c353","_cell_guid":"23acde97-d173-4c94-8c21-159d7c7c6720","trusted":true},"cell_type":"code","source":"nucleiIndex = 4\n\nprint('separation of object {}'.format(nucleiIndex))\n\nimgNucleiTest = (labels == nucleiIndex)\n\n# cropping is used only to display image bigger > works the same without cropping as separation method internally crops object\ncropy, cropx, cropy2, cropx2 = getCropMaskDimensions((imgNucleiTest == True))\nimgNucleiTest = imgNucleiTest[cropy:cropy2,cropx:cropx2]\nimshow(imgNucleiTest)\nplt.show()\n\nimgNucleiTestSeparated = separation(imgNucleiTest)\nimshow(imgNucleiTestSeparated)\nplt.show()","execution_count":4,"outputs":[]},{"metadata":{"_uuid":"7455b503774a6d40148d5d53c1651ae8510c26b9","_cell_guid":"66f55e49-fd18-4e08-85b0-42c1b0232435"},"cell_type":"markdown","source":"**Looks good. Now let's try on whole image**\n\nHold my beer :)"},{"metadata":{"_uuid":"6f7d9493f97050d891e9b2fb6233490c9df46eab","_cell_guid":"a641d40d-c5a1-4308-b51f-0ec67a251497","trusted":true},"cell_type":"code","source":"reconstructedMask = np.zeros(mask.shape, dtype=np.bool)\nfor i in range(1, labels.max() + 1):\n    # separate objects\n    img_ = separation(labels == i)\n    # copy to reconstructed mask\n    reconstructedMask = reconstructedMask + img_\n\nlabelsFinal = label(reconstructedMask, neighbors=4)\n      \nprint(\"Source mask image\")\nimshow(label2rgb(labels, bg_label=0))\nplt.show()\n\nprint(\"Processed mask image\")\nimshow(label2rgb(labelsFinal, bg_label=0))\nplt.show()","execution_count":5,"outputs":[]},{"metadata":{"_uuid":"0b7fb883b17c92b769afb93e43bfc21862ac8c35","_cell_guid":"618f1efd-4412-40aa-92ff-03cfdb82d0b7"},"cell_type":"markdown","source":"**Amazing**\n\nI would like to thank to all \"brains\" that contributes to this challenge. With solution to this challenge we'll have a better chance to test and find cures faster for ill childrens.\n\nKind regards"}],"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}