{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "9e340962-ed4a-8e0f-6877-f504cffc25ab"
      },
      "source": [
        "This script shows the full training and prediction pipeline for a pixel-based classifier: we create a mask, train logistic regression on one-pixel patches, make prediction for all pixels, create and smooth polygons from pixels."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "e388c17e-209b-97bb-a4d1-5ee5b8a9dbdc"
      },
      "outputs": [],
      "source": [
        "from collections import defaultdict\n",
        "import csv\n",
        "import sys\n",
        "\n",
        "import cv2\n",
        "\n",
        "from shapely.geometry import MultiPolygon, Polygon\n",
        "import shapely.wkt\n",
        "import shapely.affinity\n",
        "import numpy as np\n",
        "import tifffile as tiff\n",
        "\n",
        "csv.field_size_limit(sys.maxsize);"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "f419e71d-3ae0-9cdf-b10c-88049fafc541"
      },
      "source": [
        "We'll work on buildings (class 1) from image 6120_2_2. Fist load grid sizes and polygons."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "8ba03b67-e2a0-d855-2fe6-ae2f8e4ad897"
      },
      "outputs": [],
      "source": [
        "IM_ID = '6120_2_2'\n",
        "POLY_TYPE = '1'  # buildings\n",
        "\n",
        "# Load grid size\n",
        "x_max = y_min = None\n",
        "for _im_id, _x, _y in csv.reader(open('../input/grid_sizes.csv')):\n",
        "    if _im_id == IM_ID:\n",
        "        x_max, y_min = float(_x), float(_y)\n",
        "        break\n",
        "\n",
        "# Load train poly with shapely\n",
        "train_polygons = None\n",
        "for _im_id, _poly_type, _poly in csv.reader(open('../input/train_wkt_v4.csv')):\n",
        "    if _im_id == IM_ID and _poly_type == POLY_TYPE:\n",
        "        train_polygons = shapely.wkt.loads(_poly)\n",
        "        break\n",
        "\n",
        "# Read image with tiff\n",
        "im_rgb = tiff.imread('../input/three_band/{}.tif'.format(IM_ID)).transpose([1, 2, 0])\n",
        "im_size = im_rgb.shape[:2]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "738a815a-a1ee-aed4-19e5-a6f745458bbf"
      },
      "source": [
        "Scale polygons to match image:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9352582b-4da8-52bd-d951-96ae2e157416"
      },
      "outputs": [],
      "source": [
        "def get_scalers():\n",
        "    h, w = im_size  # they are flipped so that mask_for_polygons works correctly\n",
        "    w_ = w * (w / (w + 1))\n",
        "    h_ = h * (h / (h + 1))\n",
        "    return w_ / x_max, h_ / y_min\n",
        "\n",
        "x_scaler, y_scaler = get_scalers()\n",
        "\n",
        "train_polygons_scaled = shapely.affinity.scale(\n",
        "    train_polygons, xfact=x_scaler, yfact=y_scaler, origin=(0, 0, 0))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fa7268d5-9493-2c9a-c82d-f7ad7ec76087"
      },
      "outputs": [],
      "source": [
        "xlist=[]\n",
        "ylist=[]\n",
        "for poly in train_polygons_scaled:\n",
        "    x,y=poly.exterior.xy\n",
        "    xlist.append(np.count_nonzero(x))\n",
        "    ylist.append(np.count_nonzero(y))\n",
        "print(sum(xlist),sum(ylist))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "8ce835a0-2e65-4685-6001-e3aa211c05fa"
      },
      "source": [
        "Create a mask from polygons:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "49b3b539-00da-2adb-4d23-eeedc896a7eb"
      },
      "outputs": [],
      "source": [
        "def mask_for_polygons(polygons):\n",
        "    img_mask = np.zeros(im_size, np.uint8)\n",
        "    if not polygons:\n",
        "        return img_mask\n",
        "    int_coords = lambda x: np.array(x).round().astype(np.int32)\n",
        "    exteriors = [int_coords(poly.exterior.coords) for poly in polygons]\n",
        "    interiors = [int_coords(pi.coords) for poly in polygons\n",
        "                 for pi in poly.interiors]\n",
        "    cv2.fillPoly(img_mask, exteriors, 1)\n",
        "    cv2.fillPoly(img_mask, interiors, 0)\n",
        "    return img_mask, exteriors, interiors\n",
        "\n",
        "train_mask, exteriors, interiors = mask_for_polygons(train_polygons_scaled)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "db02bae2-d24a-9805-4add-c056dffce581"
      },
      "outputs": [],
      "source": [
        "len(exteriors)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a5411722-d823-a62c-5fb4-36793d9a9596"
      },
      "outputs": [],
      "source": [
        "len(interiors)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "5e175d53-bc8a-9ea0-05f5-03a6cd0083e2"
      },
      "outputs": [],
      "source": [
        "print(interiors)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "ed5d76c8-47a3-2a65-f6b4-46ef0e53af65"
      },
      "source": [
        "A helper for nicer display"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "78357c93-f936-a843-110a-009455a25830"
      },
      "outputs": [],
      "source": [
        "def scale_percentile(matrix):\n",
        "    w, h, d = matrix.shape\n",
        "    matrix = np.reshape(matrix, [w * h, d]).astype(np.float64)\n",
        "    # Get 2nd and 98th percentile\n",
        "    mins = np.percentile(matrix, 1, axis=0)\n",
        "    maxs = np.percentile(matrix, 99, axis=0) - mins\n",
        "    matrix = (matrix - mins[None, :]) / maxs[None, :]\n",
        "    matrix = np.reshape(matrix, [w, h, d])\n",
        "    matrix = matrix.clip(0, 1)\n",
        "    return matrix"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "be829690-6fcb-8bf9-6b75-6ca4808bd6c3"
      },
      "source": [
        "Check that image and mask are aligned.\n",
        "Image:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "92545b40-b831-8165-1dc7-fe395f900741"
      },
      "outputs": [],
      "source": [
        "tiff.imshow(255 * scale_percentile(im_rgb[2900:3200,2000:2300]));"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "e8c166c3-eb24-a9a3-998e-0b60cbb7dd6b"
      },
      "source": [
        "And mask:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "74b5d9ba-fdde-671c-757e-2608cf2d356f"
      },
      "outputs": [],
      "source": [
        "def show_mask(m):\n",
        "    # hack for nice display\n",
        "    tiff.imshow(255 * np.stack([m, m, m]));\n",
        "show_mask(train_mask[2900:3200,2000:2300])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d7d5b159-0deb-6154-5e42-eceb73290c9a"
      },
      "source": [
        "Now, let's train a very simple logistic regression classifier, just to get some noisy prediction to show how output mask is processed."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "36724807-e62f-6d9e-d074-2b10bdeb5f2b"
      },
      "outputs": [],
      "source": [
        "from sklearn.linear_model import SGDClassifier\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.pipeline import make_pipeline\n",
        "from sklearn.model_selection import cross_val_predict\n",
        "from sklearn.metrics import average_precision_score\n",
        "\n",
        "xs = im_rgb.reshape(-1, 3).astype(np.float32)\n",
        "ys = train_mask.reshape(-1)\n",
        "pipeline = make_pipeline(StandardScaler(), SGDClassifier(loss='log'))\n",
        "\n",
        "print('training...')\n",
        "# do not care about overfitting here\n",
        "pipeline.fit(xs, ys)\n",
        "pred_ys = pipeline.predict_proba(xs)[:, 1]\n",
        "print('average precision', average_precision_score(ys, pred_ys))\n",
        "pred_mask = pred_ys.reshape(train_mask.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "82444eb3-6619-c5e0-8b75-e84861fc02b8"
      },
      "source": [
        "Now check predictions:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "62b059c5-3b85-ff83-844f-28b0eb27f113"
      },
      "outputs": [],
      "source": [
        "show_mask(pred_mask[2900:3200,2000:2300])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "1e45d1b3-ff2a-f922-9d9f-8685043aef26"
      },
      "source": [
        "We must choose a threshold to turn it into a binary mask:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "bd6bce17-ea6a-b5b3-5e39-6973717585cf"
      },
      "outputs": [],
      "source": [
        "threshold = 0.3\n",
        "pred_binary_mask = pred_mask >= threshold\n",
        "show_mask(pred_binary_mask[2900:3200,2000:2300])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "37a76e2e-d88e-13e1-8959-9b6abc3bca09"
      },
      "source": [
        "Now it's possible to check Jaccard on the pixel level:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "da81d13e-9555-67b3-31fe-d336c01af0d3"
      },
      "outputs": [],
      "source": [
        "# check jaccard on the pixel level\n",
        "tp, fp, fn = (( pred_binary_mask &  train_mask).sum(),\n",
        "              ( pred_binary_mask & ~train_mask).sum(),\n",
        "              (~pred_binary_mask &  train_mask).sum())\n",
        "print('Pixel jaccard', tp / (tp + fp + fn))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "7dc1b380-3e44-9169-e6da-ae65982cd071"
      },
      "source": [
        "Next is the most interesting bit, creating polygons from bit masks. Please see inline comments:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ca144939-f1bd-c3c1-83c4-56d0febf7695"
      },
      "outputs": [],
      "source": [
        "def mask_to_polygons(mask, epsilon=10., min_area=10.):\n",
        "    # first, find contours with cv2: it's much faster than shapely\n",
        "    image, contours, hierarchy = cv2.findContours(\n",
        "        ((mask == 1) * 255).astype(np.uint8),\n",
        "        cv2.RETR_CCOMP, cv2.CHAIN_APPROX_TC89_KCOS)\n",
        "    # create approximate contours to have reasonable submission size\n",
        "    approx_contours = [cv2.approxPolyDP(cnt, epsilon, True)\n",
        "                       for cnt in contours]\n",
        "    if not contours:\n",
        "        return MultiPolygon()\n",
        "    # now messy stuff to associate parent and child contours\n",
        "    cnt_children = defaultdict(list)\n",
        "    child_contours = set()\n",
        "    assert hierarchy.shape[0] == 1\n",
        "    # http://docs.opencv.org/3.1.0/d9/d8b/tutorial_py_contours_hierarchy.html\n",
        "    for idx, (_, _, _, parent_idx) in enumerate(hierarchy[0]):\n",
        "        if parent_idx != -1:\n",
        "            child_contours.add(idx)\n",
        "            cnt_children[parent_idx].append(approx_contours[idx])\n",
        "    # create actual polygons filtering by area (removes artifacts)\n",
        "    all_polygons = []\n",
        "    for idx, cnt in enumerate(approx_contours):\n",
        "        if idx not in child_contours and cv2.contourArea(cnt) >= min_area:\n",
        "            assert cnt.shape[1] == 1\n",
        "            poly = Polygon(\n",
        "                shell=cnt[:, 0, :],\n",
        "                holes=[c[:, 0, :] for c in cnt_children.get(idx, [])\n",
        "                       if cv2.contourArea(c) >= min_area])\n",
        "            all_polygons.append(poly)\n",
        "    # approximating polygons might have created invalid ones, fix them\n",
        "    all_polygons = MultiPolygon(all_polygons)\n",
        "    if not all_polygons.is_valid:\n",
        "        all_polygons = all_polygons.buffer(0)\n",
        "        # Sometimes buffer() converts a simple Multipolygon to just a Polygon,\n",
        "        # need to keep it a Multi throughout\n",
        "        if all_polygons.type == 'Polygon':\n",
        "            all_polygons = MultiPolygon([all_polygons])\n",
        "    return all_polygons"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d26dc686-a0f6-979a-4cf6-1dd00541766f"
      },
      "source": [
        "Turn our prediction to polygons, and then turn back into a mask to check what it looks like:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "f8e31a8f-b2e3-7f84-f733-b0b8602520df"
      },
      "outputs": [],
      "source": [
        "pred_polygons = mask_to_polygons(pred_binary_mask)\n",
        "pred_poly_mask = mask_for_polygons(pred_polygons)\n",
        "show_mask(pred_poly_mask[2900:3200,2000:2300])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "d18792ef-7138-4fc0-6ff9-1dc3f860f104"
      },
      "source": [
        "Now to create a submission we just scale back to original coordinates"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "6f9f1e0b-465b-2e86-f059-1b76e5505c49"
      },
      "outputs": [],
      "source": [
        "scaled_pred_polygons = shapely.affinity.scale(\n",
        "    pred_polygons, xfact=1 / x_scaler, yfact=1 / y_scaler, origin=(0, 0, 0))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "9d066236-0de6-d722-8247-129a6bf32661"
      },
      "source": [
        "Checking submission size:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "7ffd7429-1931-5a7f-acd3-8f5c314b38b6"
      },
      "outputs": [],
      "source": [
        "dumped_prediction = shapely.wkt.dumps(scaled_pred_polygons)\n",
        "print('Prediction size: {:,} bytes'.format(len(dumped_prediction)))\n",
        "final_polygons = shapely.wkt.loads(dumped_prediction)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "27cbb32c-a33c-820f-9279-0a70b54f6b99"
      },
      "source": [
        "Now the litmus test: check Jaccard compared to **original** polygons"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "fa3b414c-aa0c-29d6-50c5-a2bae13e4414"
      },
      "outputs": [],
      "source": [
        "print('Final jaccard',\n",
        "      final_polygons.intersection(train_polygons).area /\n",
        "      final_polygons.union(train_polygons).area)"
      ]
    }
  ],
  "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
}