{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "ccc128d6-b369-cfbc-8a83-d49a138cc912"
      },
      "source": [
        "This kernels aims at segmenting the cervix using the technique presented in this paper: https://www.researchgate.net/publication/24041301_Automatic_Detection_of_Anatomical_Landmarks_in_Uterine_Cervix_Images"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fb67719e-6367-0dac-9c81-85acc797b305"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import cv2\n",
        "import math\n",
        "from sklearn import mixture\n",
        "from sklearn.utils import shuffle\n",
        "from skimage import measure\n",
        "from glob import glob\n",
        "import os\n",
        "from multiprocessing import Pool, cpu_count\n",
        "from functools import partial\n",
        "from subprocess import check_output\n",
        "print(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n",
        "\n",
        "TRAIN_DATA = \"../input/train\"\n",
        "\n",
        "types = ['Type_1']#,'Type_2','Type_3']\n",
        "type_ids = []\n",
        "\n",
        "for type in enumerate(types):\n",
        "    type_i_files = glob(os.path.join(TRAIN_DATA, type[1], \"*.jpg\"))\n",
        "    type_i_ids = np.array([s[len(TRAIN_DATA)+8:-4] for s in type_i_files])\n",
        "    type_ids.append(type_i_ids[:5])\n",
        "\n",
        "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)\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",
        "def 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": "markdown",
      "metadata": {
        "_cell_guid": "d905eee8-7e3a-11f4-0ef6-8aa3f5e11395"
      },
      "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": "b3f0968b-7f53-1646-b624-ad2c137a0e3f"
      },
      "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",
        "\n",
        "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",
        "\n",
        "def 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",
        "    \n",
        "    rect = maxRect(ff)\n",
        "    rectangle = [min(rect[0],rect[2]), max(rect[0],rect[2]), min(rect[1],rect[3]), max(rect[1],rect[3])]\n",
        "    img_crop = img[rectangle[0]:rectangle[1], rectangle[2]:rectangle[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",
        "    return [img_crop, rectangle, tile_size]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "2b329aab-9144-3163-3033-fde55b737cc5"
      },
      "source": [
        "\u201cFor an initial delineation of the cervix, we use two features: \n",
        "\n",
        " - the *a* color channel of the source image in Lab color space (the higher the value of *a* , the \u201credder\u201d the pixel color)\n",
        " - *R*, the distance of a pixel from the image center. The *R* feature provides spatial information and supports the extraction of continuous regions within the image plane.\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "bd01c6dc-c40c-a9c4-5c0c-f93d44568553"
      },
      "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": "markdown",
      "metadata": {
        "_cell_guid": "08069a8d-3546-64f4-74a8-514f5ecaff0d"
      },
      "source": [
        "\"The image is separated next into two clusters in the 2-D (*a-R*) feature space; we use Gaussian mixture modeling, initialized by a K-means procedure.\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "295958e3-b4f9-7943-c44d-7954a242c352"
      },
      "outputs": [],
      "source": [
        "def get_and_crop_image(image_id, image_type):\n",
        "    img = get_image_data(image_id, image_type)\n",
        "    initial_shape = img.shape\n",
        "    [img, rectangle_cropCircle, tile_size] = cropCircle(img)\n",
        "    imgLab = cv2.cvtColor(img, cv2.COLOR_RGB2LAB);\n",
        "    w = img.shape[0]\n",
        "    h = img.shape[1]\n",
        "    Ra = Ra_space(imgLab, 1.0, 150)\n",
        "    a_channel = np.reshape(Ra[:,1], (w,h))\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",
        "    cv2.drawContours(img, main_contour, -1, 255, 3)\n",
        "    \n",
        "    x,y,w,h = cv2.boundingRect(main_contour)\n",
        "    \n",
        "    rectangle = [x+rectangle_cropCircle[2],\n",
        "                 y+rectangle_cropCircle[0],\n",
        "                 w,\n",
        "                 h,\n",
        "                 initial_shape[0],\n",
        "                 initial_shape[1],\n",
        "                 tile_size[0],\n",
        "                 tile_size[1]]\n",
        "\n",
        "    return [image_id, img, rectangle]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "503a6882-7e6c-ce4f-aeb5-5766c7b6b5c5"
      },
      "outputs": [],
      "source": [
        "def parallelize_image_cropping(image_ids):\n",
        "    out = open('rectangles.csv', \"w\")\n",
        "    out.write(\"image_id,type,x,y,w,h,img_shp_0_init,img_shape1_init,img_shp_0,img_shp_1\\n\")\n",
        "    imf_d = {}\n",
        "    p = Pool(cpu_count())\n",
        "    for type in enumerate(types):\n",
        "        partial_get_and_crop = partial(get_and_crop_image, image_type = type[1])    \n",
        "        ret = p.map(partial_get_and_crop, image_ids[type[0]])\n",
        "        for i in range(len(ret)):\n",
        "            out.write(image_ids[type[0]][i])\n",
        "            out.write(',' + str(type[1]))\n",
        "            out.write(',' + str(ret[i][2][0]))\n",
        "            out.write(',' + str(ret[i][2][1]))\n",
        "            out.write(',' + str(ret[i][2][2]))\n",
        "            out.write(',' + str(ret[i][2][3]))\n",
        "            out.write(',' + str(ret[i][2][4]))\n",
        "            out.write(',' + str(ret[i][2][5]))\n",
        "            out.write(',' + str(ret[i][2][6]))\n",
        "            out.write(',' + str(ret[i][2][7]))\n",
        "            out.write('\\n')\n",
        "            img = get_image_data(image_ids[type[0]][i], type[1])\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",
        "            img = cv2.resize(img, dsize=tile_size)\n",
        "            cv2.rectangle(img,\n",
        "                          (ret[i][2][0], ret[i][2][1]), \n",
        "                          (ret[i][2][0]+ret[i][2][2], ret[i][2][1]+ret[i][2][3]),\n",
        "                          255,\n",
        "                          2)\n",
        "            plt.imshow(img)\n",
        "            plt.show()\n",
        "        ret = []\n",
        "    out.close()\n",
        "\n",
        "    return"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "745db4a2-e1ab-772a-9de7-86fc2eede735"
      },
      "outputs": [],
      "source": [
        "parallelize_image_cropping(type_ids)"
      ]
    }
  ],
  "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
}