{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"a797778f-d4b5-94c4-3629-603a80aba56d"},"source":""},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"56dfcbc9-19fd-c042-918e-978f3a24932f"},"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 \nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom PIL import Image\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\nimport matplotlib.pyplot as plt\n\nimport seaborn as sns\nfrom scipy.misc import imread\nfrom glob import glob\nimport random\n\nimport cv2\nimport math\nfrom sklearn import mixture\nfrom sklearn.utils import shuffle\nfrom skimage import measure\nfrom glob import glob\nfrom keras.preprocessing.image import ImageDataGenerator, array_to_img, img_to_array, load_img\nfrom keras.models import Sequential\nfrom keras.layers import Conv2D, MaxPooling2D\nfrom keras.layers import Activation, Dropout, Flatten, Dense\nfrom keras import backend as K\nimport os\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/train\"]).decode(\"utf8\"))\n\nsize = 256, 256\n\nim=Image.open('../input/train/Type_1/10.jpg')\nim.thumbnail(size, Image.ANTIALIAS)\n#print (im.format, im.size, im.mode)\nplt.imshow(im)\n\n\n\n\nTRAIN_DATA = \"../input/train\"\ntype_1_files = glob(os.path.join(TRAIN_DATA, \"Type_1\", \"*.jpg\"))\n\ntype_1_ids = np.array([s[len(os.path.join(TRAIN_DATA, \"Type_1\"))+1:-4] for s in type_1_files])\n\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]\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\n\n# Any results you write to the current directory are saved as output."},{"cell_type":"markdown","metadata":{"_cell_guid":"e12a0a26-cd70-8904-3c68-d78674f62254"},"source":"Number of images"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"5e44a4d9-6a2f-ea14-92fd-e49bee632204"},"outputs":[],"source":"sub_folders = check_output([\"ls\", \"../input/train/\"]).decode(\"utf8\").strip().split('\\n')\ncount_dict = {}\nfor sub_folder in sub_folders:\n    num_of_files = len(check_output([\"ls\", \"../input/train/\"+sub_folder]).decode(\"utf8\").strip().split('\\n'))\n    print(\"{0} photos of cervix type {1} \".format(num_of_files, sub_folder))\n\n    count_dict[sub_folder] = num_of_files\n    \nplt.figure(figsize=(12,4))\nsns.barplot(list(count_dict.keys()), list(count_dict.values()), alpha=0.8)\nplt.xlabel('Cervix types', fontsize=12)\nplt.ylabel('Number of Images in train', fontsize=12)\nplt.title(\"train dataset\")\n\nplt.show()"},{"cell_type":"markdown","metadata":{"_cell_guid":"c494f58c-fc13-d5ec-39a6-095448dbe719"},"source":"Look at additional"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"80780d0a-ea83-707e-b6eb-a24d184d772f"},"outputs":[],"source":"sub_folders = check_output([\"ls\", \"../input/additional/\"]).decode(\"utf8\").strip().split('\\n')\ncount_dict = {}\nfor sub_folder in sub_folders:\n    num_of_files = len(check_output([\"ls\", \"../input/additional/\"+sub_folder]).decode(\"utf8\").strip().split('\\n'))\n    print(\"{0} photos of cervix type {1} \".format(num_of_files, sub_folder))\n\n    count_dict[sub_folder] = num_of_files\n    \nplt.figure(figsize=(12,4))\nsns.barplot(list(count_dict.keys()), list(count_dict.values()), alpha=0.8)\nplt.xlabel('Cervix types', fontsize=12)\nplt.ylabel('Number of Images in additional', fontsize=12)\nplt.title(\"Additional dataset\")\nplt.show()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"4fd6089d-aa5b-052b-7986-d37ea5a89513"},"outputs":[],"source":"num_test_files = len(check_output([\"ls\", \"../input/test/\"]).decode(\"utf8\").strip().split('\\n'))\nprint(\"Number of test images present :\", num_test_files)"},{"cell_type":"markdown","metadata":{"_cell_guid":"8e1d6d65-468e-dfee-efcb-6005f46e0d75"},"source":"First, we crop the image in order to remove the circular frames that might be present. This is done by finding the largest inscribed rectangle to the thresholded image. The image is then cropped to this rectangle. (see these videos for an explanation of the algorithm: https://www.youtube.com/watch?v=g8bSdXCG-lA, https://www.youtube.com/watch?v=VNbkzsnllsU)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e4a60c3c-3644-07f4-68b4-a98cb6ba8a4f"},"outputs":[],"source":"train_path = \"../input/train/\"\nsub_folders = check_output([\"ls\", train_path]).decode(\"utf8\").strip().split('\\n')\ndifferent_file_sizes = {}\nfor sub_folder in sub_folders:\n    file_names = check_output([\"ls\", train_path+sub_folder]).decode(\"utf8\").strip().split('\\n')\n    for file_name in file_names:\n        im_array = imread(train_path+sub_folder+\"/\"+file_name)\n        size = \"_\".join(map(str,list(im_array.shape)))\n        different_file_sizes[size] = different_file_sizes.get(size,0) + 1\n\nplt.figure(figsize=(12,4))\nsns.barplot(list(different_file_sizes.values()), list(different_file_sizes.keys()), alpha=0.8)\nplt.ylabel('Image size', fontsize=12)\nplt.xlabel('Number of Images in train', fontsize=12)\nplt.title(\"Image sizes present in train dataset\")\nplt.show()"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"924d5cc2-5c65-e626-a998-6b8064bc8f5c"},"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            \n\ndef 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"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c1dccd2a-53aa-0489-4cea-9119be5b8e76"},"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":"0b2642b3-c6be-defc-337a-c3bc2f93343c"},"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()"}],"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}