{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Masks - Quick EDA - New data\n\nWhat do we learn here:\n\n- From tiny masks to large masks\n- 3 groups based on masks aera\n- Masks are in Cortex\n- Some outlier masks outside tissue or in Medulla\n- Masks footprint compared 1024x1024 tile"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import os, sys, random, gc, math, glob, time, pathlib\nimport numpy as np\nimport pandas as pd\nimport io, timeit, os, gc, pickle, psutil\nimport warnings\nimport cv2\nimport gdal\nimport osgeo\nimport json\nimport rasterio\nfrom rasterio.windows import Window\n\nimport seaborn as sns\nsns.set()\nsns.set_context(\"paper\", font_scale=1.2) \n\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nimport matplotlib.colors as cols\n\nwarnings.filterwarnings('ignore')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print('Python        : ' + sys.version.split('\\n')[0])\nprint('Numpy         : ' + np.__version__)\nprint('Pandas        : ' + pd.__version__)\nprint('Rasterio      : ' + rasterio.__version__)\nprint('GDal          : ' + osgeo.gdal.__version__)\nprint('OpenCV        : ' + cv2.__version__)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"HOME =  \"./\"\nDATA_HOME = \"../input/hubmap-kidney-segmentation/\"\nTRAIN_HOME = DATA_HOME + \"train/\"\nTEST_HOME = DATA_HOME + \"test/\"\n\nIMAGE_ID = \"image_id\"","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"anatomical_files = [os.path.basename(f) for f in glob.glob(TRAIN_HOME + \"*.json\") if \"anatomical\" in f]\nmasks_files = [os.path.basename(f) for f in glob.glob(TRAIN_HOME + \"*.json\") if \"anatomical\" not in f]\nmasks_files","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"mask_polys = {}\nfor file in masks_files:\n    with open(TRAIN_HOME + file) as jsonfile:\n        data = json.load(jsonfile)\n        key = file.split(\".\")[0]\n        mask_polys[key] = []\n        for index in range(data.__len__()):\n            if (data[index]['properties']['classification']['name'] == 'glomerulus'):             \n                geom = np.array(data[index]['geometry']['coordinates'])\n                mask_polys[key].append(geom)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Compute statistics on masks: Area and points"},{"metadata":{"trusted":true},"cell_type":"code","source":"mask_areas = {}\nmask_points = {}\nstats = []\nfor item, polys in mask_polys.items():\n    areas = []\n    points = []\n    for p in polys:\n        p = p.squeeze(axis=0)\n        area = cv2.contourArea(p.astype(np.float32))\n        areas.append(area)\n        points.append(len(p))\n    mask_areas[item] = areas\n    mask_points[item] = points\n    stats.append((item, len(polys), np.min(areas), np.max(areas), np.mean(areas), np.median(areas), np.std(areas), np.min(points), np.mean(points), np.max(points)))\npoly_pd = pd.DataFrame(stats, columns=[IMAGE_ID, \"total\", \"min\", \"max\", \"mean\", \"median\", \"std\", \"min_pts\", \"mean_pts\", \"max_pts\"])\npoly_pd = poly_pd.sort_values([\"median\"]).reset_index(drop=True)\npoly_pd = poly_pd.set_index(IMAGE_ID)\npoly_pd","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fix, ax = plt.subplots(1,2, figsize=(22, 5))\nfor item, row in poly_pd.iterrows():\n    d = sns.distplot(pd.DataFrame(mask_points[item]), ax=ax[0], label=item)\n    d.set_title(\"Mask points distribution\")\n    d = sns.distplot(pd.DataFrame(mask_areas[item]), ax=ax[1], label=item)\n    d.set_title(\"Mask area distribution\")\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_cortex_medulla(file):\n    cortex_polys = []\n    medulla_polys = []\n    with open(file) as jsonfile:\n        data = json.load(jsonfile)    \n        for index in range(data.__len__()):\n            if (data[index]['properties']['classification']['name'] == 'Cortex'):\n                geom = np.array(data[index]['geometry']['coordinates'])\n                cortex_polys.append(geom)                \n            if (data[index]['properties']['classification']['name'] == 'Medulla'):\n                geom = np.array(data[index]['geometry']['coordinates'])\n                medulla_polys.append(geom)\n    return cortex_polys, medulla_polys","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Shape with height, width\ndef make_grid(shape, window=1024, min_overlap=0):\n    \"\"\"\n        Return Array of size (N,4), where N - number of tiles,\n        2nd axis represente slices: x1,x2,y1,y2 \n    \"\"\"\n    y, x = shape\n    nx = x // (window - min_overlap) + 1\n    x1 = np.linspace(0, x, num=nx, endpoint=False, dtype=np.int64)\n    x1[-1] = x - window\n    x2 = (x1 + window).clip(0, x)\n    ny = y // (window - min_overlap) + 1\n    y1 = np.linspace(0, y, num=ny, endpoint=False, dtype=np.int64)\n    y1[-1] = y - window\n    y2 = (y1 + window).clip(0, y)\n    slices = np.zeros((nx,ny, 4), dtype=np.int64)\n    \n    for i in range(nx):\n        for j in range(ny):\n            slices[i,j] = x1[i], x2[i], y1[j], y2[j]    \n    return slices.reshape(nx*ny,4)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Display cortex (blue), medulla (green) and masks (red) for each train image.\nAlso display 1024x1024 grid (white) to see mask footprint compared to tile size "},{"metadata":{"trusted":true},"cell_type":"code","source":"for image_id, row in poly_pd.iterrows():\n    anatomical_file = TRAIN_HOME + image_id + \"-anatomical-structure.json\"\n    cortex_polys, medulla_polys = read_cortex_medulla(anatomical_file)\n        \n    with rasterio.open(TRAIN_HOME + image_id + \".tiff\") as file:\n        if file.count == 3:\n            image = file.read([1,2,3]).transpose(1,2,0).copy()\n        else:\n            h, w = (file.height, file.width)\n            subdatasets = file.subdatasets\n            if len(subdatasets) > 0:\n                image = np.zeros((h, w, len(subdatasets)), dtype=np.uint8)\n                for i, subdataset in enumerate(subdatasets, 0):\n                    with rasterio.open(subdataset) as layer:\n                        image[:,:,i] = layer.read(1) # np.moveaxis(image_, 0, -1) #.squeeze(axis=2)  \n                \n    \n    # Green = Medulla\n    if len(medulla_polys) > 0:\n        for medulla_poly in medulla_polys:\n            image = cv2.polylines(image, medulla_poly.astype(np.int32), True, (0,255,0), thickness=30)\n\n    # Blue = Cortex\n    if len(cortex_polys) > 0:\n        for cortex_poly in cortex_polys:\n            if len(cortex_poly) > 1:\n                for cortex_pts in cortex_poly:\n                    image = cv2.polylines(image, np.expand_dims(np.array(cortex_pts[0]).astype(np.int32), axis=0), True, (0,0,255), thickness=30)\n            else:\n                image = cv2.polylines(image, cortex_poly.astype(np.int32), True, (0,0,255), thickness=30)\n    \n    \n    # Red = Mask\n    for mask_poly in mask_polys[image_id]:\n        image = cv2.polylines(image, mask_poly.astype(np.int32), True, (255,0,0), thickness=30)        \n    \n    fix, ax = plt.subplots(1,1, figsize=(32, 30))\n\n    # make_grid to compare mask footprint to tile size\n    size = image.shape\n    boxes = make_grid((size[0], size[1]), window=1024)\n    \n    for i, box in enumerate(boxes):\n        x1, y1 = box[0], box[2]\n        x2, y2 = box[1], box[3]\n        image = cv2.rectangle(image, (x1, y1), (x2, y2), color=(255,255,255), thickness=8) \n\n    image = cv2.resize(image, (image.shape[1]//4, image.shape[0]//4))\n    d = ax.imshow(image)\n    d = ax.set_title(\"%s, %dx%d, masks=%d, median area=%.1f, blue=Cortex, green=Medulla\" % (image_id, size[0], size[1], len(mask_polys[image_id]), poly_pd.loc[image_id][\"median\"]))\n    d = ax.grid(None)\n    d = ax.axis('off')\n    plt.show()\n    \n    del image\n    \n    # break","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}