{"cells":[{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"79536998-ee20-d6db-c1f1-e531d25397b1"},"outputs":[],"source":"from collections import defaultdict\nimport csv\nimport sys\n\nimport cv2\nfrom shapely.geometry import MultiPolygon, Polygon\nimport shapely.wkt\nimport shapely.affinity\nimport numpy as np\nimport tifffile as tiff\n\ncsv.field_size_limit(sys.maxsize);"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e15a2ded-fe50-5bd5-4a39-82b426c01f66"},"outputs":[],"source":"IM_ID = '6120_2_2'\nPOLY_TYPE = '1'  # buildings\n\n# Load grid size\nx_max = y_min = None\nfor _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\ntrain_polygons = None\ncats = list()\nfor _im_id, _poly_type, _poly in csv.reader(open('../input/train_wkt_v4.csv')):\n    print(_im_id, _poly_type, _poly)\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\nim_rgb = tiff.imread('../input/three_band/{}.tif'.format(IM_ID)).transpose([1, 2, 0])\nim_size = im_rgb.shape[:2]"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"03fc92f5-663b-4bee-b168-9e607d68f442"},"outputs":[],"source":"cats"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b2df7aa2-8130-474e-389e-0d1cf4a2bf4f"},"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\nx_scaler, y_scaler = get_scalers()\n\ntrain_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":"b379f90d-a53c-d136-4f5a-d712e18e3c66"},"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\ntrain_mask = mask_for_polygons(train_polygons_scaled)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"78e0fc18-1d3e-caa8-83f6-589db46445c8"},"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":"code","execution_count":null,"metadata":{"_cell_guid":"6a717485-1d6f-febb-1736-c11a80023d88"},"outputs":[],"source":"tiff.imshow(255 * scale_percentile(im_rgb[2900:3200,2000:2300]));"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"5bd6061f-b35a-c660-8278-13a33e7d043a"},"outputs":[],"source":"def show_mask(m):\n    # hack for nice display\n    tiff.imshow(255 * np.stack([m, m, m]));\nshow_mask(train_mask[2900:3200,2000:2300])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"86244b79-d122-3914-ab10-f3ee1dfced41"},"outputs":[],"source":"from sklearn.linear_model import SGDClassifier\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.metrics import average_precision_score\n\nxs = im_rgb.reshape(-1, 3).astype(np.float32)\nys = train_mask.reshape(-1)\npipeline = make_pipeline(StandardScaler(), SGDClassifier(loss='log'))\n\nprint('training...')\n# do not care about overfitting here\npipeline.fit(xs, ys)\npred_ys = pipeline.predict_proba(xs)[:, 1]\nprint('average precision', average_precision_score(ys, pred_ys))\npred_mask = pred_ys.reshape(train_mask.shape)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c4e243e8-6255-e323-2635-0870965611de"},"outputs":[],"source":"show_mask(pred_mask[2900:3200,2000:2300])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"6a74ad55-5332-5a49-9c07-da4618aa5262"},"outputs":[],"source":"threshold = 0.3\npred_binary_mask = pred_mask >= threshold\nshow_mask(pred_binary_mask[2900:3200,2000:2300])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"89971b65-d979-2051-2e86-3dee82116eb8"},"outputs":[],"source":"# check jaccard on the pixel level\ntp, fp, fn = (( pred_binary_mask &  train_mask).sum(),\n              ( pred_binary_mask & ~train_mask).sum(),\n              (~pred_binary_mask &  train_mask).sum())\nprint('Pixel jaccard', tp / (tp + fp + fn))"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"caafeb11-b5df-b8a1-705d-ca3884972026"},"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":"code","execution_count":null,"metadata":{"_cell_guid":"9d40570d-9281-3d84-6b64-e14cab4089e5"},"outputs":[],"source":"pred_polygons = mask_to_polygons(pred_binary_mask)\npred_poly_mask = mask_for_polygons(pred_polygons)\nshow_mask(pred_poly_mask[2900:3200,2000:2300])"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"ef1f7e8d-0ff0-77b5-9f9f-43af9ff20a26"},"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":"code","execution_count":null,"metadata":{"_cell_guid":"a9ff14d2-a821-45f3-afdc-c00c12aa6292"},"outputs":[],"source":"dumped_prediction = shapely.wkt.dumps(scaled_pred_polygons)\nprint('Prediction size: {:,} bytes'.format(len(dumped_prediction)))\nfinal_polygons = shapely.wkt.loads(dumped_prediction)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"d4377f88-e3e7-790a-480e-3f004f792ec8"},"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}