{"nbformat": 4, "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"file_extension": ".py", "pygments_lexer": "ipython3", "codemirror_mode": {"version": 3, "name": "ipython"}, "nbconvert_exporter": "python", "version": "3.6.4", "mimetype": "text/x-python", "name": "python"}}, "nbformat_minor": 1, "cells": [{"cell_type": "markdown", "source": ["# Goal\n", "So the excellent [original kernel](https://www.kaggle.com/ahassaine/pure-image-processing-lb-0-274) put together by [Ali](https://www.kaggle.com/ahassaine) gets 0.274 \"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.274\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": {"_uuid": "c337bb4f41c90a4c45decf5d66a2aea2a41a9d7d", "_cell_guid": "86e8e6f7-b9af-4927-9a60-d9ef4d5cf059"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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'"], "metadata": {"collapsed": true, "_uuid": "b14671a1b5787c6aaa46cdcf51499ef45754ca42", "_cell_guid": "e598f333-b831-4c7b-a6c7-8a1beb532e88"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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)"], "metadata": {"_uuid": "136dd5ba7ca3b388cd5b75f7418dbddb6f1cb411", "_cell_guid": "427ce69b-1f31-444e-b4d3-1c512c141be9"}}, {"cell_type": "markdown", "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": {"_uuid": "281c2cd1244d19dd5cca2e680eb1079de7f20a47", "_cell_guid": "d900aae7-a8c7-472a-a7c1-c405c2a90d69"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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])[:,:,1]\n", "\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)"], "metadata": {"_uuid": "8e498fa24a1836ba59b23ede176b1fc7999eef0f", "_cell_guid": "ce197e18-5ba8-4994-bb01-3d2c97ecac9c"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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')"], "metadata": {"_uuid": "5c2dfc5d23bd3adf7b5504678a5311f5eecccd83", "_cell_guid": "aa5aef97-cb00-4569-a33d-bc58add3fb53"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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"], "metadata": {"_uuid": "acf47e06fd60e3c9481161f95cd3a50327224d61", "_cell_guid": "9724c70f-2fd1-42aa-82ce-4e7615283c16"}}, {"cell_type": "markdown", "source": ["Here are the two functions from the original kernel"], "metadata": {"_uuid": "a074470349efa3bc33b15cb7e7ce4b7398ca9be8", "_cell_guid": "848e1b83-be79-4530-b49b-96846d579c94"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["def ali_pipeline(img_green):\n", "    #green channel happends to produce slightly better results\n", "    #than the grayscale image and other channels\n", "    #morphological opening (size tuned on training data)\n", "    circle7=cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(7,7))\n", "    img_open=cv2.morphologyEx(img_green, cv2.MORPH_OPEN, circle7)\n", "    #Otsu thresholding\n", "    img_th=cv2.threshold(img_open,0,255,cv2.THRESH_OTSU)[1]\n", "    #Invert the image in case the objects of interest are in the dark side\n", "    if(np.sum(img_th==255)>np.sum(img_th==0)):\n", "        img_th=cv2.bitwise_not(img_th)\n", "    #second morphological opening (on binary image this time)\n", "    bin_open=cv2.morphologyEx(img_th, cv2.MORPH_OPEN, circle7) \n", "    #connected components\n", "    cc=cv2.connectedComponents(bin_open)[1]\n", "    #cc=segment_on_dt(bin_open,20)\n", "    return cc"], "metadata": {"collapsed": true, "_uuid": "7357411528d58d217118632f14abd32bfd3dc7dc", "_cell_guid": "6c81c361-2459-4105-be42-95164307058b"}}, {"cell_type": "markdown", "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": {"_uuid": "e3731bdee2baf64800b6c2f2a136654ec5738e6d", "_cell_guid": "4b11d692-1548-486d-9890-5afd65253fc0"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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 = y_pred_in\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)"], "metadata": {"collapsed": true, "_uuid": "67e1da6e4bd005ef030a94e758cf4779e7c72427", "_cell_guid": "8525ffc1-a762-41c5-a1d5-a79fdb1b9e98"}}, {"cell_type": "markdown", "source": ["# Wrap Everything Up\n", "We have a single function to evaluate the IOU of a model on a dataset"], "metadata": {"_uuid": "9c754c574083b67a3c2aa8176d9c5d34c0c0611e", "_cell_guid": "f2fac00f-1cf9-45fa-b0f1-292b042befb2"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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)"], "metadata": {"collapsed": true, "_uuid": "60725f803240b42fe16df91261da02de683af8de", "_cell_guid": "8f91efc0-10ad-4de5-b172-50096d06f405"}}, {"cell_type": "markdown", "source": ["# Estimate our IOU with the Ali Model\n", "Here we see the result is 0.42 which is quite a bit higher than the actual 0.274, so we have to be aware the value isn't super reliable."], "metadata": {"_uuid": "5eb9013a56d1fbd141a6cc42cfb4f265f9eaa744", "_cell_guid": "ba11ee11-c2fd-45f3-be81-d3aca4df9ebf"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["%%time\n", "print('IOU', calculate_iou(valid_split_df, ali_pipeline))"], "metadata": {"_uuid": "717fdc16fc9a233de6850d9c92407cce93a36755", "_cell_guid": "3ea493f6-b123-4f4a-9506-78bc2698b1af"}}, {"cell_type": "markdown", "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": {"_uuid": "8cfe7529e3e8a88ae874aa60e6049535767686c9", "_cell_guid": "749c24e0-d324-4e60-89f5-5cfc4549a2bc"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["def parametric_pipeline(img_green,\n", "                invert_thresh_pd = 10,\n", "                circle_size_x = 7,\n", "                circle_size_y = 7,\n", "                ):\n", "    circle_size_x = np.clip(int(circle_size_x), 1, 30)\n", "    circle_size_y = np.clip(int(circle_size_y), 1, 30)\n", "    \n", "    #green channel happends to produce slightly better results\n", "    #than the grayscale image and other channels\n", "    #morphological opening (size tuned on training data)\n", "    circle7=cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(circle_size_x, circle_size_y))\n", "    img_open=cv2.morphologyEx(img_green, cv2.MORPH_OPEN, circle7)\n", "    #Otsu thresholding\n", "    img_th=cv2.threshold(img_open,0,255,cv2.THRESH_OTSU)[1]\n", "    #Invert the image in case the objects of interest are in the dark side\n", "    if(np.sum(img_th==255)>((invert_thresh_pd/10.0)*np.sum(img_th==0))):\n", "        img_th=cv2.bitwise_not(img_th)\n", "    #second morphological opening (on binary image this time)\n", "    bin_open=cv2.morphologyEx(img_th, cv2.MORPH_OPEN, circle7) \n", "    #connected components\n", "    cc=cv2.connectedComponents(bin_open)[1]\n", "    #cc=segment_on_dt(bin_open,20)\n", "    return cc"], "metadata": {"collapsed": true, "_uuid": "3c4ea0f0051b03ad7ce590b51719295b550a4ce2", "_cell_guid": "ee8c0810-1e18-42db-8c38-293b7f2a4568"}}, {"cell_type": "markdown", "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. The random search is also very primitive, a better package would improve this massively."], "metadata": {"_uuid": "628cad9ff87f2f42efdfa86b5f3f0d68a030b481", "_cell_guid": "94ef8ade-8e3b-4ef9-a7b5-73237665a56e"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["from scipy.optimize import fmin_powell\n", "from tqdm import tqdm\n", "base_x0_min = [0, 1, 1]\n", "base_x0_max = [25, 30, 30]\n", "\n", "def random_search_fmin(random_restart = 5, search_steps = 5):\n", "    results = []\n", "    base_x0 = (10, 7, 7) # starting point\n", "    for _ in tqdm(range(random_restart)):\n", "        def inv_iou_func(x0):\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_params = fmin_powell(inv_iou_func, \n", "                          base_x0, \n", "                           direc = np.array([0.5, -1.5, -1.5]),\n", "                          xtol = 0.25,\n", "                          maxfun = search_steps)\n", "        \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"], "metadata": {"collapsed": true, "_uuid": "b20fb54fa7efa9613608e061bf7091c873ef2442", "_cell_guid": "1807d9c1-0f4d-4519-bf0d-156b6b61f423"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["%%time\n", "opt_params, results = random_search_fmin(5, 10)"], "metadata": {"scrolled": false, "_uuid": "3659a5ff3d3045dbb832606ebee9d3de9754e2f0", "_cell_guid": "50b04527-fbe5-4183-b76d-d4d5378de011"}}, {"cell_type": "markdown", "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": {"_uuid": "b72b36edb48fdc98fc1a84d6ee0bee2e11da1911", "_cell_guid": "71e67370-53d0-4d9b-8757-c8f93ea11f11"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["print('Opt Parameters', opt_params)\n", "print('IOU', calculate_iou(valid_split_df, \n", "                           lambda x: parametric_pipeline(x, *opt_params)))"], "metadata": {"_uuid": "2535ac4e8832d002a3b338d23cf59c97ec8222ee", "_cell_guid": "9253cd87-ab6e-4e4e-866c-6565a113d231"}}, {"cell_type": "markdown", "source": ["Now we load the test images and apply the algorithm to them"], "metadata": {"_uuid": "95338a77695bb4c45c2ae60d89c1aa19a9985b8c", "_cell_guid": "8174c3da-27bc-4803-9c44-36c73317de6f"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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)"], "metadata": {"_uuid": "63abd4d28230c8346c297168f8f298d48f6b2638", "_cell_guid": "f0cc48cd-4393-4b59-98e9-0d0ea2877dda"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["%%time\n", "test_img_df['masks'] = test_img_df['images'].map(lambda x: \n", "                                                 parametric_pipeline(x, *opt_params))"], "metadata": {"_uuid": "cc7b6df4d24c6797b50b919c4ef370e9d3a158bc", "_cell_guid": "54d69648-c9a5-4d54-b534-9b9da239fa9b"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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')"], "metadata": {"_uuid": "d3ad229d4ce913825bdcdd8bdab772a73b192001", "_cell_guid": "3d5299a3-be04-433e-b1e8-23f807b7078e"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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)"], "metadata": {"collapsed": true, "_uuid": "d45c3c014cf81f07edfb0e11942d8ff20f9662e7", "_cell_guid": "02a4a2d2-921e-433c-9234-0043a8e3b258"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["test_img_df['rles'] = test_img_df['masks'].map(lambda x: list(prob_to_rles(x)))"], "metadata": {"collapsed": true, "_uuid": "2d27f71380ee26639297bcbbc4cbf6822d126462", "_cell_guid": "006863bf-de97-4a48-8ea9-c751b879afcc"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "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)"], "metadata": {"_uuid": "2c18af608c436d9ef8a3984888a9c2c2db5cb58d", "_cell_guid": "f927cd7b-5cac-46e0-ba83-1287ceb75813"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": ["out_pred_df[['ImageId', 'EncodedPixels']].to_csv('predictions.csv', index = False)"], "metadata": {"collapsed": true, "_uuid": "032a35371b8906477f2f4d9782d16a7e63e6deec", "_cell_guid": "1c318658-1cf3-4238-b61a-8b890238b44e"}}, {"outputs": [], "cell_type": "code", "execution_count": null, "source": [], "metadata": {"collapsed": true, "_uuid": "648d79ac4e50cfd6959773cec4ee50b02f0f94f5", "_cell_guid": "1349ec1a-cfae-48b2-b522-a65d9384f369"}}]}