{"cells":[{"metadata":{"_cell_guid":"55f2b2d6-46f1-cd55-c297-ce67a2b9afbb"},"cell_type":"markdown","source":"## Hi everyone, this kernel is about applying different segmentation methods to the images of the cervix\n\nMethods used:\n- a channel saturaion\n- Watershed\n- Edge detection\n- K-means\n\nReferences can be found at the beginning of each method"},{"metadata":{"_cell_guid":"aadb4e22-aa3a-356b-6683-fbb7796e2331","trusted":false},"cell_type":"code","source":"%matplotlib inline\n\n# 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)\nimport cv2\n\n#additional imports\nimport matplotlib.pylab as plt\nimport math\nfrom sklearn import mixture\nfrom sklearn import svm\nimport pickle\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":null,"outputs":[]},{"metadata":{"_cell_guid":"26587b4d-16a1-9b0a-fdec-56996d10f0df","trusted":false},"cell_type":"code","source":"import os\nfrom glob import glob\nTRAIN_DATA = \"../input/train\"\ntype_1_files = glob(os.path.join(TRAIN_DATA, \"Type_1\", \"*.jpg\"))\ntype_1_ids = np.array([s[len(os.path.join(TRAIN_DATA, \"Type_1\"))+1:-4] for s in type_1_files])\ntype_2_files = glob(os.path.join(TRAIN_DATA, \"Type_2\", \"*.jpg\"))\ntype_2_ids = np.array([s[len(os.path.join(TRAIN_DATA, \"Type_2\"))+1:-4] for s in type_2_files])\ntype_3_files = glob(os.path.join(TRAIN_DATA, \"Type_3\", \"*.jpg\"))\ntype_3_ids = np.array([s[len(os.path.join(TRAIN_DATA, \"Type_3\"))+1:-4] for s in type_3_files])\n\nprint(len(type_1_files), len(type_2_files), len(type_3_files))\nprint(\"Type 1\", type_1_ids[:10])\nprint(\"Type 2\", type_2_ids[:10])\nprint(\"Type 3\", type_3_ids[:10])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"26ab8ec9-347b-5670-7e39-d1f742078657","trusted":false},"cell_type":"code","source":"TEST_DATA = \"../input/test\"\ntest_files = glob(os.path.join(TEST_DATA, \"*.jpg\"))\ntest_ids = np.array([s[len(TEST_DATA)+1:-4] for s in test_files])\nprint(len(test_ids))\nprint(test_ids[:10])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"6ccb7c61-fe0c-8ac9-014a-5977be0e416b","trusted":false},"cell_type":"code","source":"ADDITIONAL_DATA = \"../input/additional\"\nadditional_type_1_files = glob(os.path.join(ADDITIONAL_DATA, \"Type_1\", \"*.jpg\"))\nadditional_type_1_ids = np.array([s[len(os.path.join(ADDITIONAL_DATA, \"Type_1\"))+1:-4] for s in additional_type_1_files])\nadditional_type_2_files = glob(os.path.join(ADDITIONAL_DATA, \"Type_2\", \"*.jpg\"))\nadditional_type_2_ids = np.array([s[len(os.path.join(ADDITIONAL_DATA, \"Type_2\"))+1:-4] for s in additional_type_2_files])\nadditional_type_3_files = glob(os.path.join(ADDITIONAL_DATA, \"Type_3\", \"*.jpg\"))\nadditional_type_3_ids = np.array([s[len(os.path.join(ADDITIONAL_DATA, \"Type_3\"))+1:-4] for s in additional_type_3_files])\n\n\nprint(len(additional_type_1_files), len(additional_type_2_files), len(additional_type_3_files))\nprint(\"Type 1\", additional_type_1_ids[:10])\nprint(\"Type 2\", additional_type_2_ids[:10])\nprint(\"Type 3\", additional_type_3_ids[:10])","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"10e5836f-4144-efc3-79ce-4cc06deb0cf4","trusted":false},"cell_type":"code","source":"def get_filename(image_id, image_type):\n    \"\"\"\n    Method to get image file path from its id and type   \n    \"\"\"\n    if image_type == \"Type_1\" or \\\n        image_type == \"Type_2\" or \\\n        image_type == \"Type_3\":\n        data_path = os.path.join(TRAIN_DATA, image_type)\n    elif image_type == \"Test\":\n        data_path = TEST_DATA\n    elif image_type == \"AType_1\" or \\\n          image_type == \"AType_2\" or \\\n          image_type == \"AType_3\":\n        data_path = os.path.join(ADDITIONAL_DATA, image_type[1:])\n    else:\n        raise Exception(\"Image type '%s' is not recognized\" % image_type)\n\n    ext = 'jpg'\n    return os.path.join(data_path, \"{}.{}\".format(image_id, ext))\n\n\ndef get_image_data(image_id, image_type, rsz_ratio=1):\n    \"\"\"\n    Method to get image data as np.array specifying image id and type\n    \"\"\"\n    fname = get_filename(image_id, image_type)\n    img = cv2.imread(fname)\n    if rsz_ratio != 1:\n        img = cv2.resize(img, dsize=(int(img.shape[1] * rsz_ratio), int(img.shape[0] * rsz_ratio)))\n    assert img is not None, \"Failed to read image : %s, %s\" % (image_id, image_type)\n    img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n    return img","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"58ec7ac9-1ccc-a31e-ccb4-6f0ccfe135ba"},"cell_type":"markdown","source":"# Definitions"},{"metadata":{"_cell_guid":"c4abe511-8bef-570c-68fa-756d6f2669db","trusted":false},"cell_type":"code","source":"# defining the color of pixels outside of the roi in RGB\nmask_color = [0, 0, 0]\n\n# a channel saturation threshold\nLOWER_A_SAT = 0\nUPPER_A_SAT = 300\n\n# resize ratio for computational speed\nresize_ratio = 0.1\n\n# figure size for plotting\nfigure_size = 7\n\n# sample images\nnormal_ids = [['1414', 'Type_1'], ['60', 'Type_2'], ['491', 'Type_3']]\ndifficult_ids = [['212', 'Type_3'], ['1093', 'Type_1'], ['1473', 'Type_1'], ['267', 'Type_1'], ['446', 'Type_1']]\nids = [['1414', 'Type_1'], ['60', 'Type_2'], ['267', 'Type_1']]","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"aeb90e7f-df31-e23b-6041-741f16668c54"},"cell_type":"markdown","source":"# a_channel and distance from center gaussian mixture\n\nTaken from chattob's kernel: https://www.kaggle.com/chattob/cervix-segmentation-gmm\n\nIspired by:  https://www.researchgate.net/publication/24041301_Automatic_Detection_of_Anatomical_Landmarks_in_Uterine_Cervix_Images"},{"metadata":{"_cell_guid":"5dc7450f-5a15-4a7a-4bbc-186d26d72a65","trusted":false},"cell_type":"code","source":"def Ra_space(img, Ra_ratio=1, a_upper_threshold=UPPER_A_SAT, a_lower_threshold=LOWER_A_SAT):\n    imgLab = cv2.cvtColor(img, cv2.COLOR_RGB2LAB);\n    w = img.shape[0]\n    h = img.shape[1]\n    Ra = np.zeros((w*h, 2))\n    for i in range(w):\n        for j in range(h):\n            R = math.sqrt((w/2-i)*(w/2-i) + (h/2-j)*(h/2-j))\n            Ra[i*h+j, 0] = R\n            a = min(imgLab[i][j][1], a_upper_threshold)\n            a = max(imgLab[i][j][1], a_lower_threshold)\n            Ra[i*h+j, 1] = a\n            \n    if Ra_ratio != 1:\n        Ra[:,0] /= max(Ra[:,0])\n        Ra[:,0] *= Ra_ratio\n        Ra[:,1] /= max(Ra[:,1])\n\n    return Ra","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1d751f4f-3c1c-69a3-e786-24b0dd6e79dd","trusted":false},"cell_type":"code","source":"def crop_roi(image, display_image=False):\n    \n    # creating the R-a feature for the image\n    Ra_array = Ra_space(image)\n    \n    # k-means gaussian mixture model\n    g = mixture.GaussianMixture(n_components = 2, covariance_type = 'diag', random_state = 0, init_params = 'kmeans')\n    g.fit(Ra_array)\n    labels = g.predict(Ra_array)\n    \n    # creating the mask array and assign the correct cluster label\n    boolean_image_mask = np.array(labels).reshape(image.shape[0], image.shape[1])\n    \n    if display_image==True:\n        outer_cluster_label = boolean_image_mask[0,0]\n    \n        new_image = image.copy()\n    \n        for i in range(boolean_image_mask.shape[0]):\n            for j in range(boolean_image_mask.shape[1]):\n                if boolean_image_mask[i, j] == outer_cluster_label:\n                    new_image[i, j] = mask_color\n    \n        plt.figure(figsize=(figure_size,figure_size))\n    \n        plt.subplot(221)\n        plt.title(\"Original image\")    \n        plt.imshow(image), plt.xticks([]), plt.yticks([])\n    \n        plt.subplot(222)\n        plt.title(\"Region of interest\")\n        plt.imshow(new_image), plt.xticks([]), plt.yticks([])\n    \n        a_channel = np.reshape(Ra_array[:,1], (image.shape[0], image.shape[1]))\n        plt.subplot(223)\n        plt.title(\"a channel\")\n        plt.imshow(a_channel, cmap='gist_heat'), plt.xticks([]), plt.yticks([])\n  \n        plt.subplot(224)\n        plt.title(\"Gaussiam mixture scatter plot\")    \n        plt.scatter(Ra_array[:,0], Ra_array[:,1], c=boolean_image_mask)\n        plt.show()\n    \n    return boolean_image_mask","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"d0fa5dd3-5220-cab3-bb78-dd15f91fd5a7","trusted":false},"cell_type":"code","source":"#ids = difficult_ids\n#ids = normal_ids\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    # watershed algorithm\n    crop_roi(image, True)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"730c06ee-c470-bce0-35db-7b6292a38e78"},"cell_type":"markdown","source":"# Watershed algorithm\n\nhttp://docs.opencv.org/3.1.0/d3/db4/tutorial_py_watershed.html"},{"metadata":{"_cell_guid":"90f0f769-6b9c-2dfd-e2e4-d33099621a3b","trusted":false},"cell_type":"code","source":"def watershed(img):\n    gray = cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)\n    ret, thresh = cv2.threshold(gray,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n    \n    # noise removal\n    kernel = np.ones((3,3),np.uint8)\n    opening = cv2.morphologyEx(thresh,cv2.MORPH_OPEN,kernel, iterations = 2)\n    \n    # sure background area\n    sure_bg = cv2.dilate(opening,kernel,iterations=3)\n     \n    # Finding sure foreground area\n    dist_transform = cv2.distanceTransform(opening,cv2.DIST_L2,5)\n    ret, sure_fg = cv2.threshold(dist_transform,0.7*dist_transform.max(),255,0)\n    \n    # Finding unknown region\n    sure_fg = np.uint8(sure_fg)\n    unknown = cv2.subtract(sure_bg,sure_fg)\n    \n    plt.figure(figsize=(figure_size,figure_size))\n    \n    plt.subplot(221)\n    plt.title(\"Original image\")    \n    plt.imshow(img), plt.xticks([]), plt.yticks([])\n    \n    plt.subplot(222)\n    plt.title(\"thres\")    \n    plt.imshow(thresh, cmap=\"Greys_r\"), plt.xticks([]), plt.yticks([])\n    \n    plt.subplot(223)\n    plt.title(\"sure_fg\")    \n    plt.imshow(sure_fg, cmap=\"Greys_r\"), plt.xticks([]), plt.yticks([])\n    \n    plt.subplot(224)\n    plt.title(\"sure_bg\")    \n    plt.imshow(sure_bg, cmap=\"Greys_r\"), plt.xticks([]), plt.yticks([])\n    ","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"fdb6a077-bb15-4e59-f9b4-b95c93fcce5e","trusted":false},"cell_type":"code","source":"#ids = difficult_ids\n#ids = normal_ids\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    # watershed algorithm\n    watershed(image)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2de13909-974d-10e7-c304-01765c76116f"},"cell_type":"markdown","source":"# Edge detection\n\nhttp://www.bogotobogo.com/python/OpenCV_Python/python_opencv3_Image_Gradient_Sobel_Laplacian_Derivatives_Edge_Detection.php"},{"metadata":{"_cell_guid":"98db0378-797f-619d-d63d-395e97e5fae4","trusted":false},"cell_type":"code","source":"def edge_detection(img0):\n    # converting to gray scale\n    gray = cv2.cvtColor(img0, cv2.COLOR_BGR2GRAY)\n\n    # remove noise\n    img = cv2.GaussianBlur(gray,(3,3),0)\n\n    # convolute with proper kernels\n    laplacian = cv2.Laplacian(img,cv2.CV_64F)\n    sobelx = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=5)  # x\n    sobely = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=5)  # y\n\n    plt.figure(figsize=(figure_size,figure_size))\n    \n    plt.subplot(2,2,1),plt.imshow(img0)\n    plt.title('Original'), plt.xticks([]), plt.yticks([])\n    plt.subplot(2,2,2),plt.imshow(laplacian,cmap = 'gray')\n    plt.title('Laplacian'), plt.xticks([]), plt.yticks([])\n    plt.subplot(2,2,3),plt.imshow(sobelx,cmap = 'gray')\n    plt.title('Sobel X'), plt.xticks([]), plt.yticks([])\n    plt.subplot(2,2,4),plt.imshow(sobely,cmap = 'gray')\n    plt.title('Sobel Y'), plt.xticks([]), plt.yticks([])\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"1dd97fdc-7640-22d6-5d97-06923ce27205","trusted":false},"cell_type":"code","source":"#ids = difficult_ids\n#ids = normal_ids[:]\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    edge_detection(image)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"15a740a8-589b-25a2-9591-3eed41f17b7c"},"cell_type":"markdown","source":"# Canny Edge Detection\nhttp://www.bogotobogo.com/python/OpenCV_Python/python_opencv3_Image_Canny_Edge_Detection.php"},{"metadata":{"_cell_guid":"b479021f-10f5-4b6f-1c5e-d759dab91b2e","trusted":false},"cell_type":"code","source":"def canny_edge(img):\n    edges = cv2.Canny(img,100,200)\n\n    plt.figure(figsize=(figure_size,figure_size))\n    \n    plt.subplot(121),plt.imshow(img,cmap = 'gray')\n    plt.title('Original Image'), plt.xticks([]), plt.yticks([])\n    plt.subplot(122),plt.imshow(edges,cmap = 'gray')\n    plt.title('Edge Image'), plt.xticks([]), plt.yticks([])\n\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"f65a4d4a-038d-c4f8-67f1-60735208c521","trusted":false},"cell_type":"code","source":"#ids = difficult_ids\n#ids = normal_ids[:]\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    canny_edge(image)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"70237adb-62db-fcdc-a502-2ca88dd00c32"},"cell_type":"markdown","source":"# K-Means color clustering\n\nhttp://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_ml/py_kmeans/py_kmeans_opencv/py_kmeans_opencv.html#kmeans-opencv"},{"metadata":{"_cell_guid":"d2b588c8-5b15-7939-91e3-d657731479d2","trusted":false},"cell_type":"code","source":"def kmeans_color(img, K=8):\n    Z = img.reshape((-1,3))\n\n    # convert to np.float32\n    Z = np.float32(Z)\n    \n    # define criteria, number of clusters(K) and apply kmeans()\n    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)\n    ret,label,center=cv2.kmeans(Z,K,None,criteria,10,cv2.KMEANS_RANDOM_CENTERS)\n\n    # Now convert back into uint8, and make original image\n    center = np.uint8(center)\n    res = center[label.flatten()]\n    res2 = res.reshape((img.shape))\n\n    \n    plt.figure(figsize=(figure_size,figure_size))\n    \n    plt.subplot(1,2,1),plt.imshow(img)\n    plt.title('Original'), plt.xticks([]), plt.yticks([])\n    \n    plt.subplot(1,2,2),plt.imshow(res2)\n    plt.title('K = %i' % K), plt.xticks([]), plt.yticks([])\n        \n    plt.show()\n    \n    return res2","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"2c03e0d1-b68a-cb7e-ccae-5957b58a4daf","trusted":false},"cell_type":"code","source":"#ids = difficult_ids[:]\n#ids = normal_ids[:]\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    kmeans_color(image, K=4)\n    kmeans_color(image, K=8)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"89f9973a-212e-b231-050b-c1825317028d"},"cell_type":"markdown","source":"# K-Means with canny edge"},{"metadata":{"_cell_guid":"c88e95dd-8962-03ce-e3bb-d6783f082d0a","trusted":false},"cell_type":"code","source":"#ids = difficult_ids[:]\n#ids = normal_ids[:]\n\nfor i in range(len(ids[:])):\n    print('Loading image %i out of %i' % (i+1, len(ids[:])))\n    image_id = ids[i]\n    image = get_image_data(image_id[0], image_id[1], resize_ratio)\n    \n    image = kmeans_color(image)\n    canny_edge(image)","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"e717661d-f3bf-cfd7-7b41-1df1eda4cba3"},"cell_type":"markdown","source":"To be continued"}],"metadata":{"_change_revision":0,"_is_fork":false,"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.0"}},"nbformat":4,"nbformat_minor":1}