{"cells":[{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"a2cb36f0-e8e2-8cea-30f8-ebff9dbbe6e0"},"outputs":[],"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."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"aa055b2f-c4a9-1e90-003d-e4c2cf3ac636"},"outputs":[],"source":"import matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport math\nfrom sklearn import mixture\nfrom sklearn.utils import shuffle\nfrom skimage import measure\nfrom glob import glob\nimport os\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\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\ntype_1_ids = type_1_ids[:30]\n\ndef 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)\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\ndef get_image_data(image_id, image_type):\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    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"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"a49ac65d-895a-e946-4cc1-98a2b2d5b9f4"},"outputs":[],"source":"def maxHist(hist):\n    maxArea = (0, 0, 0)\n    height = []\n    position = []\n    for i in range(len(hist)):\n        if (len(height) == 0):\n            if (hist[i] > 0):\n                height.append(hist[i])\n                position.append(i)\n        else: \n            if (hist[i] > height[-1]):\n                height.append(hist[i])\n                position.append(i)\n            elif (hist[i] < height[-1]):\n                while (height[-1] > hist[i]):\n                    maxHeight = height.pop()\n                    area = maxHeight * (i-position[-1])\n                    if (area > maxArea[0]):\n                        maxArea = (area, position[-1], i)\n                    last_position = position.pop()\n                    if (len(height) == 0):\n                        break\n                position.append(last_position)\n                if (len(height) == 0):\n                    height.append(hist[i])\n                elif(height[-1] < hist[i]):\n                    height.append(hist[i])\n                else:\n                    position.pop()    \n    while (len(height) > 0):\n        maxHeight = height.pop()\n        last_position = position.pop()\n        area =  maxHeight * (len(hist) - last_position)\n        if (area > maxArea[0]):\n            maxArea = (area, len(hist), last_position)\n    return maxArea\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"417c0a8c-7f2a-5bc2-30cb-53f895371520"},"outputs":[],"source":"def maxRect(img):\n    maxArea = (0, 0, 0)\n    addMat = np.zeros(img.shape)\n    for r in range(img.shape[0]):\n        if r == 0:\n            addMat[r] = img[r]\n            area = maxHist(addMat[r])\n            if area[0] > maxArea[0]:\n                maxArea = area + (r,)\n        else:\n            addMat[r] = img[r] + addMat[r-1]\n            addMat[r][img[r] == 0] *= 0\n            area = maxHist(addMat[r])\n            if area[0] > maxArea[0]:\n                maxArea = area + (r,)\n    return (int(maxArea[3]+1-maxArea[0]/abs(maxArea[1]-maxArea[2])), maxArea[2], maxArea[3], maxArea[1], maxArea[0])\n\ndef cropCircle(img):\n    if(img.shape[0] > img.shape[1]):\n        tile_size = (int(img.shape[1]*256/img.shape[0]),256)\n    else:\n        tile_size = (256, int(img.shape[0]*256/img.shape[1]))\n\n    img = cv2.resize(img, dsize=tile_size)\n            \n    gray = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY);\n    _, thresh = cv2.threshold(gray, 10, 255, cv2.THRESH_BINARY)\n\n    _, contours, _ = cv2.findContours(thresh.copy(),cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)\n\n    main_contour = sorted(contours, key = cv2.contourArea, reverse = True)[0]\n            \n    ff = np.zeros((gray.shape[0],gray.shape[1]), 'uint8') \n    cv2.drawContours(ff, main_contour, -1, 1, 15)\n    ff_mask = np.zeros((gray.shape[0]+2,gray.shape[1]+2), 'uint8')\n    cv2.floodFill(ff, ff_mask, (int(gray.shape[1]/2), int(gray.shape[0]/2)), 1)\n    #cv2.circle(ff, (int(gray.shape[1]/2), int(gray.shape[0]/2)), 3, 3, -1)\n    \n    rect = maxRect(ff)\n    img_crop = img[min(rect[0],rect[2]):max(rect[0],rect[2]), min(rect[1],rect[3]):max(rect[1],rect[3])]\n    cv2.rectangle(ff,(min(rect[1],rect[3]),min(rect[0],rect[2])),(max(rect[1],rect[3]),max(rect[0],rect[2])),3,2)\n\n    #plt.subplot(121)\n    #plt.imshow(img)\n    #plt.subplot(122)\n    #plt.imshow(ff)\n    #plt.show()\n    \n    return img_crop\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"7615b690-77fb-29a0-a4a0-5a8957c3105a"},"outputs":[],"source":"def Ra_space(img, Ra_ratio, a_threshold):\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            Ra[i*h+j, 1] = min(imgLab[i][j][1], a_threshold)\n            \n    Ra[:,0] /= max(Ra[:,0])\n    Ra[:,0] *= Ra_ratio\n    Ra[:,1] /= max(Ra[:,1])\n\n    return Ra"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"4bb8cfaa-40ec-0e3b-d672-057d3428d2d1"},"outputs":[],"source":"for k, type_ids in enumerate([type_1_ids]):\n    m = len(type_ids)\n    train_ids = sorted(type_ids)\n    counter = 0\n    \n    for i in range(m):                \n        image_id = train_ids[counter] \n        counter += 1\n\n        img = get_image_data(image_id, 'Type_%i' % (k+1))\n\n        img = cropCircle(img)\n        w = img.shape[0]\n        h = img.shape[1]\n                        \n        imgLab = cv2.cvtColor(img, cv2.COLOR_RGB2LAB);\n        \n        # Saturating the a-channel at 150 helps avoiding wrong segmentation\n        # in the case of close-up cervix pictures where the bloody os is falsly segemented as the cervix.\n        Ra = Ra_space(img, 1.0, 150) \n        a_channel = np.reshape(Ra[:,1], (w,h))\n        plt.subplot(121)\n        plt.imshow(a_channel) \n\n        g = mixture.GaussianMixture(n_components = 2, covariance_type = 'diag', random_state = 0, init_params = 'kmeans')\n        image_array_sample = shuffle(Ra, random_state=0)[:1000]\n        g.fit(image_array_sample)\n        labels = g.predict(Ra)\n        labels += 1 # Add 1 to avoid labeling as 0 since regionprops ignores the 0-label.\n    \n        # The cluster that has the highest a-mean is selected.\n        labels_2D = np.reshape(labels, (w,h))\n        gg_labels_regions = measure.regionprops(labels_2D, intensity_image = a_channel)\n        gg_intensity = [prop.mean_intensity for prop in gg_labels_regions]\n        cervix_cluster = gg_intensity.index(max(gg_intensity)) + 1\n\n        mask = np.zeros((w * h,1),'uint8')\n        mask[labels==cervix_cluster] = 255\n        mask_2D = np.reshape(mask, (w,h))\n\n        cc_labels = measure.label(mask_2D, background=0)\n        regions = measure.regionprops(cc_labels)\n        areas = [prop.area for prop in regions]\n\n        regions_label = [prop.label for prop in regions]\n        largestCC_label = regions_label[areas.index(max(areas))]\n        mask_largestCC = np.zeros((w,h),'uint8')\n        mask_largestCC[cc_labels==largestCC_label] = 255\n\n        img_masked = img.copy()\n        img_masked[mask_largestCC==0] = (0,0,0)\n        img_masked_gray = cv2.cvtColor(img_masked, cv2.COLOR_RGB2GRAY);\n            \n        _,thresh_mask = cv2.threshold(img_masked_gray,0,255,0)\n            \n        kernel = np.ones((11,11), np.uint8)\n        thresh_mask = cv2.dilate(thresh_mask, kernel, iterations = 1)\n        thresh_mask = cv2.erode(thresh_mask, kernel, iterations = 2)\n        _, contours_mask, _ = cv2.findContours(thresh_mask.copy(),cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)\n\n        main_contour = sorted(contours_mask, key = cv2.contourArea, reverse = True)[0]\n                    \n        x,y,w,h = cv2.boundingRect(main_contour)\n        cv2.rectangle(img,(x,y),(x+w,y+h),255,2)\n                        \n        plt.subplot(122)\n        plt.imshow(img)\n        plt.show()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"6cc5602b-5398-f931-f47c-d7e4cc5c711f"},"outputs":[],"source":"def maxHist(hist):\n    maxArea = (0, 0, 0)\n    height = []\n    position = []\n    for i in range(len(hist)):\n        if (len(height) == 0):\n            if (hist[i] > 0):\n                height.append(hist[i])\n                position.append(i)\n        else: \n            if (hist[i] > height[-1]):\n                height.append(hist[i])\n                position.append(i)\n            elif (hist[i] < height[-1]):\n                while (height[-1] > hist[i]):\n                    maxHeight = height.pop()\n                    area = maxHeight * (i-position[-1])\n                    if (area > maxArea[0]):\n                        maxArea = (area, position[-1], i)\n                    last_position = position.pop()\n                    if (len(height) == 0):\n                        break\n                position.append(last_position)\n                if (len(height) == 0):\n                    height.append(hist[i])\n                elif(height[-1] < hist[i]):\n                    height.append(hist[i])\n                else:\n                    position.pop()    \n    while (len(height) > 0):\n        maxHeight = height.pop()\n        last_position = position.pop()\n        area =  maxHeight * (len(hist) - last_position)\n        if (area > maxArea[0]):\n            maxArea = (area, len(hist), last_position)\n    return maxArea\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1cd4d1bc-0fa7-2bb9-85dd-652c3d045760"},"outputs":[],"source":"fname = get_filename(0,'Type_1')\nimg = cv2.imread(fname)\nimg = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\naddMat = np.zeros(img.shape)\n\naddMat[0] = img[0]\narea = maxHist(addMat[0])"}],"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":0}