{"metadata":{"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,"cells":[{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0f3bc2a8-2c2f-4d19-fb4f-ca957e651fab","_active":false},"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.\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport cv2\nimport math\nimport scipy.ndimage as ndi\nfrom sklearn import mixture\nfrom sklearn.utils import shuffle\nfrom skimage import measure\nfrom skimage.color import rgb2gray\nfrom glob import glob\nimport os\n\nfrom multiprocessing import Pool, cpu_count\nfrom sklearn.feature_extraction import image\nfrom functools import partial\nfrom sklearn import datasets, svm, metrics\nfrom sklearn.ensemble import (RandomForestClassifier, ExtraTreesClassifier,\n                              AdaBoostClassifier)\nfrom sklearn.datasets import load_iris\nfrom numpy import genfromtxt\n\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n\nTRAIN_DATA = \"../input/train\"\n\ntypes = ['Type_1','Type_2','Type_3']\ntype_ids = []\n\nfor 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])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"7d7d0f43-16cc-a9b1-8ad1-fac350cc6141","_active":false},"outputs":[],"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)\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))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"69898004-f9e1-58cd-1231-5341a8e450e1","_active":false},"outputs":[],"source":"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":"code","execution_count":null,"metadata":{"_cell_guid":"8bd2472b-e028-b0f8-cd64-c2da9708bd4d","_active":false},"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"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"12aac910-2f1a-64bd-8185-0df94fd96462","_active":false},"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])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"31c7931d-5e82-1cb5-17f3-9c85e88b9a87","_active":false},"outputs":[],"source":"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":"code","execution_count":null,"metadata":{"_cell_guid":"5cf3534a-82ca-b79a-709d-57fe8a477917","_active":true},"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":"47f504ab-1ebf-b254-e830-b0e42719f6f6","_active":false},"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":"8eca2eac-4328-7a80-5f22-849a2c34e670","_active":false},"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    ret = []\n    \n    plt_counter = 1\n    fig = plt.figure(figsize=(50, 50))\n    \n    for type in enumerate(types):\n        partial_get_and_crop = partial(get_and_crop_image, image_type = type[1])   \n\n        for image_id in image_ids[type[0]]:\n            ret.append(partial_get_and_crop(image_id))\n        \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 = (192, 256)\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            crop_img = img[ret[i][2][1]:ret[i][2][1]+ret[i][2][3],ret[i][2][0]:ret[i][2][0]+ret[i][2][2]]\n            crop_img = cv2.resize(crop_img, dsize=(192, 256))\n            \n            mask = np.zeros(img.shape,np.uint8)\n            mask[ret[i][2][1]:ret[i][2][1]+ret[i][2][3],ret[i][2][0]:ret[i][2][0]+ret[i][2][2]] = img[ret[i][2][1]:ret[i][2][1]+ret[i][2][3],ret[i][2][0]:ret[i][2][0]+ret[i][2][2]]\n            \n            ax = fig.add_subplot(all_samples, 10, plt_counter)\n            ax.imshow(img)\n\n            ax = fig.add_subplot(all_samples, 10, plt_counter+1)\n            ax.imshow(crop_img)\n\n            ax = fig.add_subplot(all_samples, 10, plt_counter+2)\n            ax.imshow(mask)\n\n            plt_counter += 3\n        \n            if i > train_samples:\n                test_data.append(rgb2gray(crop_img).flatten())\n                test_target.append(type[1])\n            else:\n                train_data.append(rgb2gray(crop_img).flatten())\n                train_target.append(type[1])\n        ret = []\n    out.close()\n    \n    plt.show()\n    \n    return"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d6eba809-9019-94f5-3d4b-15b87feda713","_active":false},"outputs":[],"source":"def model_random_forest(train_features, train_target, test_features, test_target):\n    random_forest = RandomForestClassifier(n_estimators=30)\n    random_forest.fit(train_features, train_target)\n    \n\n    random_forest_predicted = random_forest.predict(test_features)\n    random_forest_probability = random_forest.predict_proba(test_features)\n\n    print(metrics.classification_report(test_target, random_forest_predicted))\n    print(metrics.confusion_matrix(test_target, random_forest_predicted))\n    print(test_target)\n    print(random_forest_predicted)\n    print(random_forest_probability)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"a42cf981-03af-d237-ece4-5ccc09a0f303","_active":false},"outputs":[],"source":"all_samples = 10\ntrain_samples = 7\n\ntype_ids = []\n\nfor 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[:all_samples])\n    \ntrain_data = []\ntrain_target = []\ntest_data = []\ntest_target = []\n\nparallelize_image_cropping(type_ids)\n\nprint(len(train_data))\nprint(len(train_target))\nprint(len(test_data))\nprint(len(test_target))\n\nmodel_random_forest(train_data,train_target,test_data,test_target )"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"30eaf985-fe82-ba7e-2076-a3a9f44ab822","_active":false},"outputs":[],"source":"def model_random_forest(train_features, train_target, test_features, test_target):\n    random_forest = RandomForestClassifier(n_estimators=30)\n    random_forest.fit(train_features, train_target)\n    \n\n    random_forest_predicted = random_forest.predict(test_features)\n    random_forest_probability = random_forest.predict_proba(test_features)\n\n    print(metrics.classification_report(test_target, random_forest_predicted))\n    print(metrics.confusion_matrix(test_target, random_forest_predicted))\n    print(test_target)\n    print(random_forest_predicted)\n    print(random_forest_probability)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"befbbc8e-b1c8-20aa-67a9-3d34a8919331","_active":false},"outputs":[],"source":"df = pd.read_csv('rectangles.csv', sep=',',header=0)\nnumeric_headers = list(df.columns.values)\nprint(numeric_headers)\nnumpy_array = df.as_matrix()\n\nX = numpy_array[:,1][:, None]\ny = numpy_array[:,2:10]\n\nmsk = np.random.rand(len(numpy_array)) < 0.8\ntrain = numpy_array[msk]\ntest = numpy_array[~msk]"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"5a4793f4-8b7d-4055-5585-d8d7cf8046c3","_active":false},"outputs":[],"source":"model_random_forest(train[:,2:10],train[:,1][:, None],test[:,2:10],test[:,1][:, None])"}]}