{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"dc0d5f94-4604-1484-7f66-89250caa1cba"},"source":"This notebook is not going to use polygons from the training set, we're going to use raster data only =)\n\nLet's prepare the playground first"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"08e681f5-3988-032e-5838-046eb2632301"},"outputs":[],"source":"import os\nimport numpy as np\nimport tifffile as tiff\nimport cv2\nimport matplotlib.pyplot as plt\nfrom skimage.segmentation import slic, mark_boundaries\n\ndef stretch_8bit(bands, lower_percent=2, higher_percent=98):\n    out = np.zeros_like(bands)\n    for i in range(3):\n        a = 0 #np.min(band)\n        b = 255  #np.max(band)\n        c = np.percentile(bands[:,:,i], lower_percent)\n        d = np.percentile(bands[:,:,i], higher_percent)        \n        t = a + (bands[:,:,i] - c) * (b - a) / (d - c)    \n        t[t<a] = a\n        t[t>b] = b\n        out[:,:,i] =t\n    return out.astype(np.uint8)    \n    \ndef RGB(image_id):\n    filename = os.path.join('..', 'input', 'three_band', '{}.tif'.format(image_id))\n    img = tiff.imread(filename)\n    img = np.rollaxis(img, 0, 3)    \n    return img\n    \ndef M(image_id):\n    filename = os.path.join('..', 'input', 'sixteen_band', '{}_M.tif'.format(image_id))\n    img = tiff.imread(filename)    \n    img = np.rollaxis(img, 0, 3)\n    return img"},{"cell_type":"markdown","metadata":{"_cell_guid":"16a12784-895e-8b81-9fad-502593d305b1"},"source":"Load image and adjust contrast a little bit."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"373546b9-58d2-5230-8e88-5f46d232e64e"},"outputs":[],"source":"    image_id = '6010_4_1'\n    rgb = RGB(image_id)\n    rgb1 = stretch_8bit(rgb)"},{"cell_type":"markdown","metadata":{"_cell_guid":"7275c6fd-1b65-9fb4-28cf-aa1e99a90907"},"source":"Region of interest\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0b1e2dc3-f92b-e168-69c5-ce57fe487939"},"outputs":[],"source":"    \n    y1,y2,x1,x2 = 0, 2160, 0, 3300\n    region = rgb1[y1:, x1:, :]\n    plt.figure()\n    plt.imshow(region)"},{"cell_type":"markdown","metadata":{"_cell_guid":"ff5d09d1-42c1-8b2b-eb85-939f3484232c"},"source":"Now, instead of RGB we'll use something a little bit different :-)\n"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c44674f9-163c-1341-327f-c334017e01b7"},"outputs":[],"source":"    m = M(image_id)    \n    m = cv2.resize(m, tuple(reversed(rgb.shape[:2])))\n    \n    img = np.zeros_like(rgb)\n    img[:,:,0] = m[:,:,6] #nir1\n    img[:,:,1] = m[:,:,4] #red\n    img[:,:,2] = rgb[:,:,2] #blue\n    img = stretch_8bit(img)\n    region = img[y1:, x1:, :]\n    plt.figure()\n    plt.imshow(region)"},{"cell_type":"markdown","metadata":{"_cell_guid":"878f6c16-523b-e084-a14a-a48d4071b924"},"source":"Trees are red. To start, let's get lazy and forget about yellow and purple rooftops  and focus on blue ones."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"98ccf409-c79e-08a2-d24b-b59b518913a6"},"outputs":[],"source":"    blue_mask = cv2.inRange(region, np.array([15,115,200]), np.array([80,200,255]))    \n    mask = cv2.bitwise_and(region, region, mask=blue_mask)\n    mask = cv2.cvtColor(mask, cv2.COLOR_RGB2GRAY)    \n    plt.figure()\n    plt.imshow(mask, cmap='gray')"},{"cell_type":"markdown","metadata":{"_cell_guid":"7858a7d1-f0f6-11e0-ab6f-ae38710c8da1"},"source":"The idea is to create a mask with regions we know for sure are blue rooftops. Hold on..."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"cd4d20bd-4082-4be8-48e9-8122bbdb4f7f"},"outputs":[],"source":"    segments = slic(region, n_segments=100, compactness=20.0, \n                    max_iter=10, sigma=5, spacing=None, multichannel=True, \n                    convert2lab=True, enforce_connectivity=False, \n                    min_size_factor=10, max_size_factor=3, slic_zero=False)\n    boundaries = mark_boundaries(region, segments, color=(0,255,0))\n    plt.figure()\n    plt.imshow(boundaries)"},{"cell_type":"markdown","metadata":{"_cell_guid":"292da583-9a8d-0ba6-19b0-5afcf22e5953"},"source":"Some magic from scikit-image :-) They're segmenting based on colour and distance using k-means. Of course the parameters need to be adjusted and depend on the size of region, but not bad as a starting point."},{"cell_type":"markdown","metadata":{"_cell_guid":"46dc2334-deaa-20fa-f33d-de98d743f5e4"},"source":"Of those regions, let's select the blue ones..."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"b59531c6-b92c-a76d-9e8c-a93872a66d8f"},"outputs":[],"source":"    out = np.zeros_like(mask)\n    for i in range(np.max(segments)):\n        s = segments == i\n        s_size = np.sum(s)\n        s_count = np.sum([1 for x in mask[s].ravel() if x>0])\n        #print(s_count, s_size)\n        if s_count > 0.1*s_size:\n            out[s] = 255\n        \n    plt.figure()\n    plt.imshow(out, cmap='gray')\n    "},{"cell_type":"markdown","metadata":{"_cell_guid":"37f3a75d-08cc-4632-d7d9-2c959ab3a0c1"},"source":"And the final output:"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"9f113c50-a926-90fe-9396-ca35ca4ded22"},"outputs":[],"source":"    out2 = cv2.bitwise_and(region, region, mask=out)\n    fig, ax = plt.subplots(1, 2)\n    ax[0].imshow(region)\n    ax[1].imshow(out2)"},{"cell_type":"markdown","metadata":{"_cell_guid":"4ad47bd9-d970-c4e8-d6fb-b0c4d7b8bb79"},"source":"You could add a purple mask and a yellow mask using cv2.bitwise_or(mask1, mask2), etc."}],"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}