{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Handlabels prize\n\nAt some point of the competition we decided to handlabel the public set because we discovered new type of glomeruli not available in train set.\nSuch glomeruli could be partially healthy (Fibrous Crescent). The contreversy outlier image is **d488c759a** which contains a lot a these glomeruli.\nAny model trained on trainset only scored bad on d488c759a. A few questions raised:\n- Is it possible that d488c759a was badly annotated?\n- If not then why trainset does not contain such type of glomeruli?\n- Is private set consistent with public set?\n\nAccording to this comment in forum:\nhttps://www.kaggle.com/c/hubmap-kidney-segmentation/discussion/234083#1281244\n\n> All images, public and private alike, went through the same annotation pipeline using the same criteria. Therefore, the annotations should be consistent. Since HuBMAP's focus is on creating a reference atlas of healthy human tissue, globally sclerotic glomeruli were excluded. However, level of injury is on a continuous scale, so some annotated glomeruli are not perfectly healthy. I believe this includes the fibrous crescent glomeruli you referred to above.\n \nand this one:\nhttps://www.kaggle.com/c/hubmap-kidney-segmentation/discussion/207884\n> We, the hackathon organizers, are taking the necessary step of completely updating the private test data. We have also added new quality control measures to the glomeruli annotation masks.\n\nWe decided that it was worth to handlabel as data should be consistent. \n\nHere are our pseudo + handlabels for public set with the following results on public LB:\n\n- 3589adb90 = 0.191+\n- 2ec3f1bb9 = 0.191+\n- aa05346ff = 0.189++\n- 57512b7f1 = 0.190+\n- d488c759a = 0.187++\n\nRelated virtual public LB from 0.948 to 0.953.\n\nAnd some examples of different glomeruli type from d488c759a.","metadata":{}},{"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 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')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Python        : ' + sys.version.split('\\n')[0])\nprint('Numpy         : ' + np.__version__)\nprint('Pandas        : ' + pd.__version__)\nprint('Rasterio      : ' + rasterio.__version__)\nprint('OpenCV        : ' + cv2.__version__)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"HOME =  \"./\"\nDATA_HOME = \"../input/hubmap-public-pseudohandlabels/\"\nTEST_HOME = \"../input/hubmap-kidney-segmentation/test/\"\n\nIMAGE_ID = \"image_id\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def read_masks(image_id, file, drop_fp=True, drop_ph=False):\n    mask_polys = {}\n    with open(file) as jsonfile:\n        data = json.load(jsonfile)\n        mask_polys[image_id] = []\n        for index in range(data.__len__()):\n            if (drop_fp is True) and (data[index]['properties'].get('name') is not None) and ('FP' in data[index]['properties'].get('name')):\n                pass # print(\"Dropping FP mask\")\n            elif (drop_ph is True) and (data[index]['properties'].get('name') is not None) and ('PH' in data[index]['properties'].get('name')):\n                pass # print(\"Dropping PH mask\")                \n            else:                \n                geom = np.array(data[index]['geometry']['coordinates'])\n                mask_polys[image_id].append(geom)\n    return mask_polys\n\n# Convert polygons to mask image\ndef polygons_to_mask(mask_polys_, h, w):    \n    generated_mask = np.zeros((h, w), dtype=np.uint8)\n    for mask_poly in mask_polys_:\n        generated_mask = cv2.fillPoly(generated_mask, pts=mask_poly.astype(np.int32), color=(255,255,255))\n    return generated_mask","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"anatomical_files = [os.path.basename(f) for f in glob.glob(DATA_HOME + \"*.json\") if \"anatomical\" in f]\nmasks_files = [os.path.basename(f) for f in glob.glob(DATA_HOME + \"*.json\") if \"anatomical\" not in f]\nmasks_files","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mask_polys = {}\nfor file in masks_files:      \n    image_id = file.replace(\".json\", \"\")\n    mask_polys[image_id] = read_masks(image_id, DATA_HOME + file)[image_id]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Compute statistics on masks: Area and points","metadata":{}},{"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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pseudo/Handlabels plots","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":{}},{"cell_type":"code","source":"for image_id, row in poly_pd.iterrows():\n    anatomical_file = TEST_HOME + image_id + \"-anatomical-structure.json\"\n    cortex_polys, medulla_polys = read_cortex_medulla(anatomical_file)\n        \n    with rasterio.open(TEST_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)\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","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some examples from d488c759a","metadata":{}},{"cell_type":"code","source":"image_id = \"d488c759a\"\n\nmask_polys_ph = {}\nmask_polys_normal = {}\nwith open(DATA_HOME + image_id + \".json\") as jsonfile:\n    data = json.load(jsonfile)\n    mask_polys_normal[image_id] = []\n    mask_polys_ph[image_id] = []\n    for index in range(data.__len__()):\n        if (data[index]['properties'].get('name') is not None) and ('FP' in data[index]['properties'].get('name')):\n            pass # print(\"Dropping FP mask\")\n        elif (data[index]['properties'].get('name') is not None) and ('PH' in data[index]['properties'].get('name')):\n            geom = np.array(data[index]['geometry']['coordinates'])\n            mask_polys_ph[image_id].append(geom)              \n        else:                \n            geom = np.array(data[index]['geometry']['coordinates'])\n            mask_polys_normal[image_id].append(geom)     ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with rasterio.open(TEST_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)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Some outliers (partially healthy or higher level of injury)","metadata":{}},{"cell_type":"code","source":"ROWS = 8\nCOLS = 6\nfig, ax = plt.subplots(ROWS, COLS, figsize=(20, 25))\ni = 0\nPAD = 40\nfor mask_poly in mask_polys_ph[image_id]:\n    x, y, w, h = cv2.boundingRect(mask_poly.astype(np.int32))\n    image = cv2.polylines(image, mask_poly.astype(np.int32), True, (255,255,255), thickness=2)    \n    img = image[y-PAD:y+h+PAD, x-PAD:x+w+PAD, :]\n    r = i%ROWS\n    c = i//ROWS\n    d = ax[r, c].imshow(img)\n    d = ax[r, c].grid(None)\n    d = ax[r, c].axis('off')\n    i = i + 1\n    if i >= ROWS*COLS: break\nplt.subplots_adjust(wspace=0.1, hspace=0.1)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Normal (healthy or low level of injury)","metadata":{}},{"cell_type":"code","source":"ROWS = 6\nCOLS = 6\nfig, ax = plt.subplots(ROWS, COLS, figsize=(20, 20))\ni = 0\nPAD = 40\nfor mask_poly in mask_polys_normal[image_id]:\n    x, y, w, h = cv2.boundingRect(mask_poly.astype(np.int32))\n    image = cv2.polylines(image, mask_poly.astype(np.int32), True, (255,255,255), thickness=2)    \n    img = image[y-PAD:y+h+PAD, x-PAD:x+w+PAD, :]\n    r = i%ROWS\n    c = i//ROWS\n    d = ax[r, c].imshow(img)\n    d = ax[r, c].grid(None)\n    d = ax[r, c].axis('off')\n    i = i + 1\n    if i >= ROWS*COLS: break\nplt.subplots_adjust(wspace=0.1, hspace=0.1)\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}