{"cells": [{"source": ["# Goal\n", "So the excellent [original kernel](https://www.kaggle.com/gaborvecsei/basic-pure-computer-vision-segmentation-lb-0-229) put together by [Gabor](https://www.kaggle.com/gaborvecsei) gets 0.229 \"without even using the training data\" which is a great result, but what if we use the training data.  \n", "## Overview\n", "The idea is to take the parameters found in the original kernel and try to improve them using the IOU score as the ground criteria.\n", "1. Get a cross-validation setup working that gives us a similar value to the 0.229\n", "1. Rewrite the threshold and label methods to take all of their parameters\n", "1. Use scipy.optimize (probably sk-optimize would be better, but we'll keep it simple here) to improve the values\n", "1. Predict and submit"], "metadata": {"_cell_guid": "86e8e6f7-b9af-4927-9a60-d9ef4d5cf059", "_uuid": "c337bb4f41c90a4c45decf5d66a2aea2a41a9d7d"}, "cell_type": "markdown"}, {"source": ["from os.path import join\n", "import cv2\n", "import numpy as np # linear algebra\n", "import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n", "from glob import glob\n", "import os\n", "from skimage.io import imread\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "%matplotlib inline\n", "dsb_data_dir = os.path.join('..', 'input')\n", "stage_label = 'stage1'"], "execution_count": null, "metadata": {"_cell_guid": "e598f333-b831-4c7b-a6c7-8a1beb532e88", "_uuid": "b14671a1b5787c6aaa46cdcf51499ef45754ca42", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["all_images = glob(os.path.join(dsb_data_dir, 'stage1_*', '*', '*', '*'))\n", "img_df = pd.DataFrame({'path': all_images})\n", "img_id = lambda in_path: in_path.split('/')[-3]\n", "img_type = lambda in_path: in_path.split('/')[-2]\n", "img_group = lambda in_path: in_path.split('/')[-4].split('_')[1]\n", "img_stage = lambda in_path: in_path.split('/')[-4].split('_')[0]\n", "img_df['ImageId'] = img_df['path'].map(img_id)\n", "img_df['ImageType'] = img_df['path'].map(img_type)\n", "img_df['TrainingSplit'] = img_df['path'].map(img_group)\n", "img_df['Stage'] = img_df['path'].map(img_stage)\n", "img_df.sample(2)"], "execution_count": null, "metadata": {"_cell_guid": "427ce69b-1f31-444e-b4d3-1c512c141be9", "_uuid": "136dd5ba7ca3b388cd5b75f7418dbddb6f1cb411"}, "cell_type": "code", "outputs": []}, {"source": ["# Process and Import Training Data\n", "Here we load in the training data images and labels. We load the label images into a single index colored integer image."], "metadata": {"_cell_guid": "d900aae7-a8c7-472a-a7c1-c405c2a90d69", "_uuid": "281c2cd1244d19dd5cca2e680eb1079de7f20a47"}, "cell_type": "markdown"}, {"source": ["%%time\n", "train_df = img_df.query('TrainingSplit==\"train\"')\n", "train_rows = []\n", "group_cols = ['Stage', 'ImageId']\n", "for n_group, n_rows in train_df.groupby(group_cols):\n", "    c_row = {col_name: col_value for col_name, col_value in zip(group_cols, n_group)}\n", "    c_row['masks'] = n_rows.query('ImageType == \"masks\"')['path'].values.tolist()\n", "    c_row['images'] = n_rows.query('ImageType == \"images\"')['path'].values.tolist()\n", "    train_rows += [c_row]\n", "train_img_df = pd.DataFrame(train_rows)    \n", "IMG_CHANNELS = 3\n", "def read_and_stack(in_img_list):\n", "    return np.sum(np.stack([i*(imread(c_img)>0) for i, c_img in enumerate(in_img_list, 1)], 0), 0)\n", "\n", "def read_hist_bw(in_img_list):\n", "    return cv2.imread(in_img_list[0], cv2.IMREAD_GRAYSCALE)\n", "train_img_df['images'] = train_img_df['images'].map(read_hist_bw)\n", "train_img_df['masks'] = train_img_df['masks'].map(read_and_stack).map(lambda x: x.astype(int))\n", "train_img_df.sample(1)"], "execution_count": null, "metadata": {"_cell_guid": "ce197e18-5ba8-4994-bb01-3d2c97ecac9c", "_uuid": "8e498fa24a1836ba59b23ede176b1fc7999eef0f"}, "cell_type": "code", "outputs": []}, {"source": ["for _, c_row in train_img_df.sample(1).iterrows():\n", "    fig, (ax1, ax2) = plt.subplots(1,2, figsize = (8, 4))\n", "    ax1.imshow(c_row['images'], cmap = 'bone')\n", "    ax2.imshow(c_row['masks'], cmap = 'nipy_spectral')"], "execution_count": null, "metadata": {"_cell_guid": "aa5aef97-cb00-4569-a33d-bc58add3fb53", "_uuid": "5c2dfc5d23bd3adf7b5504678a5311f5eecccd83"}, "cell_type": "code", "outputs": []}, {"source": ["from sklearn.model_selection import train_test_split\n", "train_split_df, valid_split_df = train_test_split(train_img_df, \n", "                                                  test_size = 0.4, \n", "                                                  random_state = 2018,\n", "                                                  # ensures both splits have the different sized images\n", "                                                  stratify = train_img_df['images'].map(lambda x: '{}'.format(np.shape))\n", "                                                 )\n", "print('train', train_split_df.shape, 'valid', valid_split_df.shape)\n"], "execution_count": null, "metadata": {"_cell_guid": "9724c70f-2fd1-42aa-82ce-4e7615283c16", "_uuid": "acf47e06fd60e3c9481161f95cd3a50327224d61"}, "cell_type": "code", "outputs": []}, {"source": ["Here are the two functions from the original kernel"], "metadata": {"_cell_guid": "848e1b83-be79-4530-b49b-96846d579c94", "_uuid": "a074470349efa3bc33b15cb7e7ce4b7398ca9be8"}, "cell_type": "markdown"}, {"source": ["def gabor_threshold(image_gray):\n", "    image_gray = cv2.GaussianBlur(image_gray, (7, 7), 1)\n", "    ret, thresh = cv2.threshold(image_gray, 0, 255, cv2.THRESH_OTSU)\n", "    \n", "    _, cnts, _ = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n", "    cnts = sorted(cnts, key=cv2.contourArea, reverse=True)\n", "    max_cnt_area = cv2.contourArea(cnts[0])\n", "    \n", "    if max_cnt_area > 50000:\n", "        ret, thresh = cv2.threshold(image_gray, 0, 255, cv2.THRESH_OTSU | cv2.THRESH_BINARY_INV)\n", "    \n", "    return thresh\n", "\n", "def gabor_apply_morphology(thresh):\n", "    mask = cv2.dilate(thresh, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5)))\n", "    mask = cv2.erode(mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5)))\n", "    return mask\n", "\n", "def gabor_pipeline(in_img):\n", "    thresh = gabor_threshold(in_img)\n", "    return gabor_apply_morphology(thresh)"], "execution_count": null, "metadata": {"_cell_guid": "6c81c361-2459-4105-be42-95164307058b", "_uuid": "7357411528d58d217118632f14abd32bfd3dc7dc", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["> # IOU Metric\n", "[This kernel](https://www.kaggle.com/aglotero/another-iou-metric) has a nice IOU implementation that we just copy here"], "metadata": {"_cell_guid": "4b11d692-1548-486d-9890-5afd65253fc0", "_uuid": "e3731bdee2baf64800b6c2f2a136654ec5738e6d"}, "cell_type": "markdown"}, {"source": ["from skimage.morphology import label\n", "def iou_metric(y_true_in, y_pred_in, print_table=False):\n", "    labels = y_true_in\n", "    y_pred = label(y_pred_in > 0.5)\n", "    \n", "    true_objects = len(np.unique(labels))\n", "    pred_objects = len(np.unique(y_pred))\n", "\n", "    intersection = np.histogram2d(labels.flatten(), y_pred.flatten(), bins=(true_objects, pred_objects))[0]\n", "\n", "    # Compute areas (needed for finding the union between all objects)\n", "    area_true = np.histogram(labels, bins = true_objects)[0]\n", "    area_pred = np.histogram(y_pred, bins = pred_objects)[0]\n", "    area_true = np.expand_dims(area_true, -1)\n", "    area_pred = np.expand_dims(area_pred, 0)\n", "\n", "    # Compute union\n", "    union = area_true + area_pred - intersection\n", "\n", "    # Exclude background from the analysis\n", "    intersection = intersection[1:,1:]\n", "    union = union[1:,1:]\n", "    union[union == 0] = 1e-9\n", "\n", "    # Compute the intersection over union\n", "    iou = intersection / union\n", "\n", "    # Precision helper function\n", "    def precision_at(threshold, iou):\n", "        matches = iou > threshold\n", "        true_positives = np.sum(matches, axis=1) == 1   # Correct objects\n", "        false_positives = np.sum(matches, axis=0) == 0  # Missed objects\n", "        false_negatives = np.sum(matches, axis=1) == 0  # Extra objects\n", "        tp, fp, fn = np.sum(true_positives), np.sum(false_positives), np.sum(false_negatives)\n", "        return tp, fp, fn\n", "\n", "    # Loop over IoU thresholds\n", "    prec = []\n", "    if print_table:\n", "        print(\"Thresh\\tTP\\tFP\\tFN\\tPrec.\")\n", "    for t in np.arange(0.5, 1.0, 0.05):\n", "        tp, fp, fn = precision_at(t, iou)\n", "        if (tp + fp + fn) > 0:\n", "            p = tp / (tp + fp + fn)\n", "        else:\n", "            p = 0\n", "        if print_table:\n", "            print(\"{:1.3f}\\t{}\\t{}\\t{}\\t{:1.3f}\".format(t, tp, fp, fn, p))\n", "        prec.append(p)\n", "    \n", "    if print_table:\n", "        print(\"AP\\t-\\t-\\t-\\t{:1.3f}\".format(np.mean(prec)))\n", "    return np.mean(prec)"], "execution_count": null, "metadata": {"_cell_guid": "8525ffc1-a762-41c5-a1d5-a79fdb1b9e98", "_uuid": "67e1da6e4bd005ef030a94e758cf4779e7c72427", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["# Wrap Everything Up\n", "We have a single function to evaluate the IOU of a model on a dataset"], "metadata": {"_cell_guid": "f2fac00f-1cf9-45fa-b0f1-292b042befb2", "_uuid": "9c754c574083b67a3c2aa8176d9c5d34c0c0611e"}, "cell_type": "markdown"}, {"source": ["def calculate_iou(in_df, thresh_func):\n", "    pred_masks = valid_split_df['images'].map(thresh_func).values\n", "    gt_masks = valid_split_df['masks'].values\n", "    all_ious = [iou_metric(cur_gt, cur_pred, print_table=False) for cur_gt, cur_pred in \n", "            zip(gt_masks, pred_masks)]\n", "    return np.mean(all_ious)"], "execution_count": null, "metadata": {"_cell_guid": "8f91efc0-10ad-4de5-b172-50096d06f405", "_uuid": "60725f803240b42fe16df91261da02de683af8de", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["# Estimate our IOU with the Gabor Model\n", "Here we see the result is 0.35 which is quite a bit higher than the actual 0.229, so we have to be aware the value isn't super reliable."], "metadata": {"_cell_guid": "ba11ee11-c2fd-45f3-be81-d3aca4df9ebf", "_uuid": "5eb9013a56d1fbd141a6cc42cfb4f265f9eaa744"}, "cell_type": "markdown"}, {"source": ["%%time\n", "print('IOU', calculate_iou(valid_split_df, gabor_pipeline))"], "execution_count": null, "metadata": {"_cell_guid": "3ea493f6-b123-4f4a-9506-78bc2698b1af", "_uuid": "717fdc16fc9a233de6850d9c92407cce93a36755"}, "cell_type": "code", "outputs": []}, {"source": ["# Define a Parametric Model\n", "Here we take the basic Gabor pipeline and allow the important parameters to be adjusted so they can consequently be optimized"], "metadata": {}, "cell_type": "markdown"}, {"source": ["def parametric_pipeline(image_gray, \n", "                       blur_sigma = 1, \n", "                        dilate_iters = 1,\n", "                       dilate_size = 5, \n", "                       erode_size = 5):\n", "    # some functions don't like floats\n", "    blur_sigma = np.clip(blur_sigma, 0.01, 100)\n", "    blur_size = int(2*round(3*blur_sigma)+1)\n", "    dilate_size = int(dilate_size)\n", "    erode_size = int(erode_size)\n", "    dilate_iters = int(dilate_iters)\n", "    if blur_size>0:\n", "        image_gray = cv2.GaussianBlur(image_gray, (blur_size, blur_size), blur_sigma)\n", "    \n", "    ret, thresh = cv2.threshold(image_gray, 0, 255, cv2.THRESH_OTSU)\n", "    _, cnts, _ = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n", "    cnts = sorted(cnts, key=cv2.contourArea, reverse=True)\n", "    max_cnt_area = cv2.contourArea(cnts[0])\n", "    \n", "    if max_cnt_area > 50000:\n", "        ret, thresh = cv2.threshold(image_gray, 0, 255, cv2.THRESH_OTSU | cv2.THRESH_BINARY_INV)\n", "    mask = thresh\n", "    for i in range(dilate_iters):\n", "        if dilate_size>0:\n", "            mask = cv2.dilate(mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (dilate_size, dilate_size)))\n", "        if erode_size>0:\n", "            mask = cv2.erode(mask, cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (erode_size, erode_size)))\n", "    return mask"], "execution_count": null, "metadata": {"_cell_guid": "ee8c0810-1e18-42db-8c38-293b7f2a4568", "_uuid": "3c4ea0f0051b03ad7ce590b51719295b550a4ce2", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["# Optimization\n", "A very simple optimization routine with no knowledge about morphology, integer steps, iterations or anything else. It just serves as an example of how such a pipeline can be optimized. At the very least some random search would probably improve the results"], "metadata": {"_cell_guid": "94ef8ade-8e3b-4ef9-a7b5-73237665a56e", "_uuid": "628cad9ff87f2f42efdfa86b5f3f0d68a030b481"}, "cell_type": "markdown"}, {"source": ["from scipy.optimize import fmin\n", "from tqdm import tqdm\n", "base_x0_min = [0, 0, 0, 0]\n", "base_x0_max = [10, 10, 15, 15]\n", "\n", "def random_search_fmin(random_restart = 5, search_steps = 5):\n", "    results = []\n", "    base_x0 = (1, 1, 5, 5) # starting point\n", "    for _ in tqdm(range(random_restart)):\n", "        def inv_iou_func(scale_x0):\n", "            x0 = [s*x/10 for s,x in zip(scale_x0, base_x0)]\n", "            # score it on the training data\n", "            try:\n", "                score = calculate_iou(train_split_df, \n", "                                      lambda x: parametric_pipeline(x, *x0))\n", "            except Exception as e:\n", "                print('Arguments:', ' '.join(['%1.1f' % xi for xi in x0]))\n", "                raise ValueError(e)\n", "            print('Arguments:', ' '.join(['%1.1f' % xi for xi in x0]), \n", "                  'IOU: %2.3f' % score)\n", "            return 1-score # since we are minimizing the result\n", "\n", "        opt_scalars = fmin(inv_iou_func, \n", "                          (10, 10, 10, 10), \n", "                           xtol = 0.1,\n", "                          maxiter = search_steps)\n", "        \n", "        opt_params = [s*x/10 for s,x in zip(opt_scalars, base_x0)]\n", "        results += [(calculate_iou(train_split_df, \n", "                                  lambda x: parametric_pipeline(x, *opt_params)), \n", "                     opt_params)]\n", "        # pick a new random spot to iterate from\n", "        base_x0 = [np.random.choice(np.linspace(x_start, x_end, 10))\n", "                   for x_start, x_end in zip(base_x0_min, base_x0_max)]\n", "    n_out = sorted(results, key = lambda x: 1-x[0])\n", "    return n_out[0][1], n_out\n"], "execution_count": null, "metadata": {"_cell_guid": "1807d9c1-0f4d-4519-bf0d-156b6b61f423", "_uuid": "b20fb54fa7efa9613608e061bf7091c873ef2442"}, "cell_type": "code", "outputs": []}, {"source": ["%%time\n", "opt_params, results = random_search_fmin(5, 8)"], "execution_count": null, "metadata": {}, "cell_type": "code", "outputs": []}, {"source": ["# Calculate the Score on Hold-Out (validation)\n", "Here we calculate the score on the validation to see if we actually improved anything"], "metadata": {"_cell_guid": "71e67370-53d0-4d9b-8757-c8f93ea11f11", "_uuid": "b72b36edb48fdc98fc1a84d6ee0bee2e11da1911"}, "cell_type": "markdown"}, {"source": ["print('IOU', calculate_iou(valid_split_df, \n", "                           lambda x: parametric_pipeline(x, *opt_params)))"], "execution_count": null, "metadata": {"_cell_guid": "9253cd87-ab6e-4e4e-866c-6565a113d231", "_uuid": "2535ac4e8832d002a3b338d23cf59c97ec8222ee"}, "cell_type": "code", "outputs": []}, {"source": ["Now we load the test images and apply the algorithm to them"], "metadata": {"_cell_guid": "8174c3da-27bc-4803-9c44-36c73317de6f", "_uuid": "95338a77695bb4c45c2ae60d89c1aa19a9985b8c"}, "cell_type": "markdown"}, {"source": ["%%time\n", "test_df = img_df.query('TrainingSplit==\"test\"')\n", "test_rows = []\n", "group_cols = ['Stage', 'ImageId']\n", "for n_group, n_rows in test_df.groupby(group_cols):\n", "    c_row = {col_name: col_value for col_name, col_value in zip(group_cols, n_group)}\n", "    c_row['images'] = n_rows.query('ImageType == \"images\"')['path'].values.tolist()\n", "    test_rows += [c_row]\n", "test_img_df = pd.DataFrame(test_rows)    \n", "\n", "test_img_df['images'] = test_img_df['images'].map(read_hist_bw)\n", "print(test_img_df.shape[0], 'images to process')\n", "test_img_df.sample(1)"], "execution_count": null, "metadata": {"_cell_guid": "f0cc48cd-4393-4b59-98e9-0d0ea2877dda", "_uuid": "63abd4d28230c8346c297168f8f298d48f6b2638"}, "cell_type": "code", "outputs": []}, {"source": ["%%time\n", "test_img_df['masks'] = test_img_df['images'].map(lambda x: \n", "                                                 parametric_pipeline(x, *opt_params))"], "execution_count": null, "metadata": {"_cell_guid": "54d69648-c9a5-4d54-b534-9b9da239fa9b", "_uuid": "cc7b6df4d24c6797b50b919c4ef370e9d3a158bc"}, "cell_type": "code", "outputs": []}, {"source": ["n_img = 3\n", "fig, m_axs = plt.subplots(2, n_img, figsize = (12, 6))\n", "for (_, d_row), (c_im, c_lab) in zip(test_img_df.sample(n_img).iterrows(), \n", "                                     m_axs.T):\n", "    c_im.imshow(d_row['images'])\n", "    c_im.axis('off')\n", "    c_im.set_title('Microscope')\n", "    \n", "    c_lab.imshow(d_row['masks'])\n", "    c_lab.axis('off')\n", "    c_lab.set_title('Predicted')"], "execution_count": null, "metadata": {"_cell_guid": "3d5299a3-be04-433e-b1e8-23f807b7078e", "_uuid": "d3ad229d4ce913825bdcdd8bdab772a73b192001"}, "cell_type": "code", "outputs": []}, {"source": ["def rle_encoding(x):\n", "    '''\n", "    x: numpy array of shape (height, width), 1 - mask, 0 - background\n", "    Returns run length as list\n", "    '''\n", "    dots = np.where(x.T.flatten()==1)[0] # .T sets Fortran order down-then-right\n", "    run_lengths = []\n", "    prev = -2\n", "    for b in dots:\n", "        if (b>prev+1): run_lengths.extend((b+1, 0))\n", "        run_lengths[-1] += 1\n", "        prev = b\n", "    return run_lengths\n", "\n", "def prob_to_rles(x, cut_off = 0.5):\n", "    lab_img = label(x>cut_off)\n", "    if lab_img.max()<1:\n", "        lab_img[0,0] = 1 # ensure at least one prediction per image\n", "    for i in range(1, lab_img.max()+1):\n", "        yield rle_encoding(lab_img==i)"], "execution_count": null, "metadata": {"_cell_guid": "02a4a2d2-921e-433c-9234-0043a8e3b258", "_uuid": "d45c3c014cf81f07edfb0e11942d8ff20f9662e7", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["test_img_df['rles'] = test_img_df['masks'].map(lambda x: list(prob_to_rles(x)))"], "execution_count": null, "metadata": {"_cell_guid": "006863bf-de97-4a48-8ea9-c751b879afcc", "_uuid": "2d27f71380ee26639297bcbbc4cbf6822d126462", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": ["out_pred_list = []\n", "for _, c_row in test_img_df.iterrows():\n", "    for c_rle in c_row['rles']:\n", "        out_pred_list+=[dict(ImageId=c_row['ImageId'], \n", "                             EncodedPixels = ' '.join(np.array(c_rle).astype(str)))]\n", "out_pred_df = pd.DataFrame(out_pred_list)\n", "print(out_pred_df.shape[0], 'regions found for', test_img_df.shape[0], 'images')\n", "out_pred_df.sample(3)"], "execution_count": null, "metadata": {"_cell_guid": "f927cd7b-5cac-46e0-ba83-1287ceb75813", "_uuid": "2c18af608c436d9ef8a3984888a9c2c2db5cb58d"}, "cell_type": "code", "outputs": []}, {"source": ["out_pred_df[['ImageId', 'EncodedPixels']].to_csv('predictions.csv', index = False)"], "execution_count": null, "metadata": {"_cell_guid": "1c318658-1cf3-4238-b61a-8b890238b44e", "_uuid": "032a35371b8906477f2f4d9782d16a7e63e6deec", "collapsed": true}, "cell_type": "code", "outputs": []}, {"source": [], "execution_count": null, "metadata": {"_cell_guid": "1349ec1a-cfae-48b2-b522-a65d9384f369", "_uuid": "648d79ac4e50cfd6959773cec4ee50b02f0f94f5", "collapsed": true}, "cell_type": "code", "outputs": []}], "nbformat": 4, "metadata": {"kernelspec": {"name": "python3", "language": "python", "display_name": "Python 3"}, "language_info": {"version": "3.6.4", "nbconvert_exporter": "python", "file_extension": ".py", "pygments_lexer": "ipython3", "codemirror_mode": {"version": 3, "name": "ipython"}, "name": "python", "mimetype": "text/x-python"}}, "nbformat_minor": 1}