{ "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", "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": "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\n", "\n", "train_mask = mask_for_polygons(train_polygons_scaled)" ] }, { "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", "\n", "mask = np.zeros_like(im_rgb)\n", "for i in range(3): \n", " mask[:,:,i] = train_mask.copy()\n", "\n", "masked_image = im_rgb[mask]\n", " \n", "tiff.imshow(masked_image[2900:3200,2000:2300]);\n", "#show_mask(train_mask[2900:3200,2000:2300])\n", "\n" ] }, { "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\n" ] }, { "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.5.2" } }, "nbformat": 4, "nbformat_minor": 0 }