{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Split adjoining cell into subsets of non-touching ones\n\nThis is the problem of graph coloring:\n\"it is a way of coloring the vertices of a graph such that no two adjacent vertices are of the same color [1]\"\n\nThe code is modified from n-color algorithm in [2] and cellpose [3].\n\nIf the cells are not touching each other, then instance segmentation is simply semantic segmentation and can be solved using UNet, etc. Hence a solution to the \"Sartorius Cell Instance Segmentation\" problem is:\n- detect cell potential candidates (e.g. centeroid of bounding box or proposal, seeds for watershed, etc)\n\n  image --> [instance detection]  --> seed\n  \n- split instances into non-overlaping subsets\n\n  seed --> [this code]  --> group1, group2 ...\n  \n- perform semantic segmentation on each subset. Use connected components labelling CCL (e.g. skimage.measure.label function) to recover instance segmentation\n\n  group-n + image --> [semantic segmentation with marker]  --> segmentation-n  --> [CCL] --> instance segmentation\n\n\nAs a final note:\n1. If the cells are not really clustered, random spilting into N subsets can also ensure they are not touching within each subset, if N is large enough. This may be useful in training.\n\n2. It is also possible to learn a network to do the subset spliting. \n\n3. This framework can be use as an add-on refinement stage of other methods like mask-rcnn. You already have a good estimate of the instance segmentation and would like to refine their masks at the full resolution. You can just split predictions of mask rcnn into non-touching subset and perform semantic segmentation with marker.\n\n[1] https://en.wikipedia.org/wiki/Graph_coloring\n\n[2] https://forum.image.sc/t/relabel-with-4-colors-like-map/33564\n\n[3] https://github.com/mouseland/cellpose\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport cv2\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom numba import njit\nfrom scipy.ndimage import generate_binary_structure\nfrom skimage.segmentation import expand_labels","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-12-03T01:08:10.708972Z","iopub.execute_input":"2021-12-03T01:08:10.709262Z","iopub.status.idle":"2021-12-03T01:08:10.714395Z","shell.execute_reply.started":"2021-12-03T01:08:10.709223Z","shell.execute_reply":"2021-12-03T01:08:10.713716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\ndef make_neighbour_idx(w, h):\n    n = generate_binary_structure(rank=2, connectivity=2)\n    n[1,1] = 0\n\n    idx = np.where(n>0)\n    idx = np.array(idx, dtype=np.uint8).T #neighbour (i,j)\n    idx = np.array(idx-[1,1])\n    idx = idx*[[w,1]]\n    idx = idx[:,0]-idx[:,1]\n    return idx\n\n@njit(fastmath=True)\ndef search_for_neighbour(image, idx):\n    line = image.ravel()\n    L = len(line)\n\n    neighbour = np.zeros((L,2), image.dtype)\n    s = 0\n    for i in range(L):\n        if line[i]==0:continue\n        for d in idx:\n            if line[i+d]==0: continue\n            if line[i]==line[i+d]: continue\n\n            neighbour[s,0] = line[i]\n            neighbour[s,1] = line[i+d]\n            s += 1\n    return neighbour[:s]\n\n\ndef find_connect(image):\n    pad = np.pad(image, 1, 'constant')\n    h,w = pad.shape\n\n    idx = make_neighbour_idx(w, h)\n    #print(idx)\n    neighbour = search_for_neighbour(pad, idx)\n    #print(neighbour)\n\n    if len(neighbour)<2:\n        return neighbour\n\n\n    neighbour.sort(axis=1)\n    i = (neighbour[:,0]<<16) + neighbour[:,1]\n\n    argsort = np.argsort(i)\n\n    i[:] = i[argsort]\n    diff = i[:-1]!=i[1:]\n\n    idx = np.where(diff)[0]+1\n    idx = np.hstack(([0], idx))\n    connect = neighbour[argsort][idx]\n    return connect\n\n\ndef make_graph(connect):\n    graph = {}\n    for i in np.unique(connect):\n        graph[i] = []\n\n    for i,j in connect:\n        graph[i].append(j)\n        graph[j].append(i)\n\n    return graph\n\n\ndef assign_color( graph, num_color=4, rand=12, depth=0, max_depth=8):\n    threshold = 1e4\n\n    if depth<max_depth:\n        node = list(graph.keys())\n        num_node = len(node)\n\n        np.random.seed(depth+1)\n        np.random.shuffle(node)\n\n        assign  = dict(zip(node, [0]*num_node))\n        counter = dict(zip(node, [0]*num_node))\n        t = 0\n        while len(node)>0 and t<threshold:\n            t+=1\n\n            k = node.pop(0)\n            counter[k] += 1\n            hist = [1e4] + [0] * num_color #history of used color\n            for p in graph[k]:\n                hist[assign[p]] += 1\n\n            if min(hist)==0:\n                assign[k] = hist.index(min(hist))\n                counter[k] = 0\n\n\n            hist[assign[k]] = 1e4\n            minc = hist.index(min(hist))\n            if counter[k]==rand:\n                counter[k] = 0\n                np.random.seed(count)\n                minc = np.random.randint(1,num_color+1)\n\n            assign[k] = minc\n            for p in graph[k]:\n                if assign[p] == minc:\n                    node.append(p)\n\n        if t==threshold:\n            print('n-color algorithm failed, trying again with num_color=%d, depth=%d'%(num_color+1,depth+1))\n            assign = assign_color(graph, num_color+1, rand, depth+1, max_depth)\n\n        return assign\n\n    else:\n        print('n-color algorithm exceeded max depth of',max_depth)\n        return None\n\n\n#---------------------------------------------\n\n#4-color algorthm based on https://forum.image.sc/t/relabel-with-4-colors-like-map/33564\ndef do_color_label(label, num_color=4):\n    # label format is 0,1,2,3... N\n\n    connect = find_connect(label)\n    graph   = make_graph(connect)\n    assign  = assign_color(graph, num_color)\n\n    lut = np.ones(label.max() + 1, dtype=np.uint8)\n    for i in assign:\n        lut[i] = assign[i]\n    lut[0] = 0\n\n    unique, inverse, count = np.unique(lut, return_inverse=True, return_counts=True)\n    lut = inverse.reshape(lut.shape)\n\n    color = lut[label]\n    return color\n\n\ndef relabel(label):\n    _, inv, area = np.unique(label, return_inverse=True, return_counts=True)\n    label = inv.reshape(label.shape)  # relabel to 0,1,2,3... N\n    return label, area\n\n\n\ndef draw_label_to_overlay(label, color=None):\n    if color is None:\n        color = np.array([\n            [  0,  0,  0],\n            [ 77,159,255],\n            [  0,255,  0],\n            [255,  0,  0],\n            [  0,255,255],\n            [255,255,  0],\n            [255,150,255],\n            #[234,178,200],\n            [  0,  0,255],\n        ])\n\n    h,w = label.shape\n    overlay = color[label]\n    return overlay\n\n","metadata":{"execution":{"iopub.status.busy":"2021-12-03T01:08:10.716885Z","iopub.execute_input":"2021-12-03T01:08:10.717345Z","iopub.status.idle":"2021-12-03T01:08:10.74886Z","shell.execute_reply.started":"2021-12-03T01:08:10.717298Z","shell.execute_reply":"2021-12-03T01:08:10.747471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"image_width  = 704\nimage_height = 520\n\ndef image_show_norm(image, mode='gray'):\n    image = (image-image.min())/(image.max()-image.min()+0.0001)\n    if mode=='gray':\n        plt.imshow(image,'gray')\n        \n    if mode=='rgb':\n        plt.imshow(image[...,::-1])\n        \n        \ndef rle_decode(rle, width=image_width, height=image_height, fill=1, dtype=np.float32):\n    s = rle.split()\n    start  = np.asarray(s[0::2], dtype=int)-1\n    length = np.asarray(s[1::2], dtype=int)\n    end = start + length\n    image = np.zeros(height * width, dtype=dtype)\n    for s, e in zip(start, end):\n        image[s:e] = fill\n    image = image.reshape(height, width) #.T\n    return image\n\n\n","metadata":{"execution":{"iopub.status.busy":"2021-12-03T01:14:05.306762Z","iopub.execute_input":"2021-12-03T01:14:05.307092Z","iopub.status.idle":"2021-12-03T01:14:05.31524Z","shell.execute_reply.started":"2021-12-03T01:14:05.307054Z","shell.execute_reply":"2021-12-03T01:14:05.314179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's try an input image \"1c10ee85de67\"","metadata":{}},{"cell_type":"code","source":"train_df = pd.read_csv('../input/sartorius-cell-instance-segmentation/train.csv')\ndf = train_df[train_df['id']=='1c10ee85de67'].reset_index(drop=True)\nlabel = np.zeros((image_height,image_width), dtype=np.int32)\nfor i,d in df.iterrows():\n    m = rle_decode(d.annotation, fill=True, dtype=np.bool)\n    label[m]=i+1\n\nimage_show_norm(label,'gray')","metadata":{"execution":{"iopub.status.busy":"2021-12-03T01:15:39.750852Z","iopub.execute_input":"2021-12-03T01:15:39.75113Z","iopub.status.idle":"2021-12-03T01:15:40.431428Z","shell.execute_reply.started":"2021-12-03T01:15:39.751099Z","shell.execute_reply":"2021-12-03T01:15:40.430558Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"After apply graph coloring,","metadata":{}},{"cell_type":"code","source":"assert len(np.unique(label)) == label.max()+1 #make sure label format: 0,1,2,3...N, is correct\n\nlabel5 = do_color_label(label, num_color=5) \nprint('input label :', np.unique(label5))\n\nexpand = expand_labels(label, distance=20)\nexpand5 = do_color_label(expand, num_color=5)\nexpand5 = expand5*(label>0)\nexpand5_overlay = draw_label_to_overlay(expand5)\n\nimage_show_norm(expand5_overlay,'rgb')","metadata":{"execution":{"iopub.status.busy":"2021-12-03T01:16:01.625531Z","iopub.execute_input":"2021-12-03T01:16:01.625817Z","iopub.status.idle":"2021-12-03T01:16:02.130159Z","shell.execute_reply.started":"2021-12-03T01:16:01.625771Z","shell.execute_reply":"2021-12-03T01:16:02.129289Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('assigned color :', np.unique(expand5))\nfor u in np.unique(expand5)[1:]:\n    m = expand5_overlay.copy()\n    m[expand5!=u]=(0,0,0)\n    plt.figure(), image_show_norm(m,'rgb')\n","metadata":{"execution":{"iopub.status.busy":"2021-12-03T01:20:02.089712Z","iopub.execute_input":"2021-12-03T01:20:02.090163Z","iopub.status.idle":"2021-12-03T01:20:03.912599Z","shell.execute_reply.started":"2021-12-03T01:20:02.090127Z","shell.execute_reply":"2021-12-03T01:20:03.911712Z"},"trusted":true},"execution_count":null,"outputs":[]}]}