{"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":"<center>\n    <h1>HuBMAP - Quick Exploratory Data Analysis</h1>\n<center>","metadata":{}},{"cell_type":"markdown","source":"<center>\n<img src=\"https://hubmapconsortium.org/wp-content/uploads/2019/01/HuBMAP-Retina-Logo-Color.png\">\n</center>","metadata":{}},{"cell_type":"markdown","source":"## Introduction\n\nThis is a quick exploratory analysis with the to get familiar with the dataset and to identify possible hurdles the might pop up down the line.\n\n> Credit to the original [notebook](https://www.kaggle.com/code/yuriikochurovskyi/hubmap-image-eda-step-by-step-beginner-friendly) which I based this one off.","metadata":{}},{"cell_type":"code","source":"import os\nimport cv2\nimport tifffile\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-06-22T19:58:13.911642Z","iopub.execute_input":"2022-06-22T19:58:13.912029Z","iopub.status.idle":"2022-06-22T19:58:13.918498Z","shell.execute_reply.started":"2022-06-22T19:58:13.912Z","shell.execute_reply":"2022-06-22T19:58:13.917465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='section1'></a>\n    \n## Importing and processing image data\n### The Dataset\nThe dataset is comprised of TIFF files. \n- The training set is a collection of \".tiff\" files\n- The public test set has some additional \".tiff\" files. \n\nThe training set includes annotations in both RLE-encoded and unencoded (JSON) forms. The annotations denote segmentations of glomeruli. \n\nFile **train.csv** contains the unique IDs for each image, as well as an RLE-encoded representation of the mask for the objects in the image. \n\nRLE or Run Length Encoding converts a matrix into a vector and returns the position/starting point of the first pixel from where we observe an object (identified by a 1) and gives us a count of how many pixels from that pixel we see the series of 1s. For example coded Mask will look like [1 1 1 0 0 1 1], running RLE would give us 1 3 6 2, which means 3 pixels from the zeroth pixel (inclusive) and 2 pixels from the 5th pixel we see a series of 1s\n\nFor the begining, let's open end review **train.csv**, it contains all RLE-masks related to each images_IDs","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('../input/hubmap-organ-segmentation/train.csv')\n\nimage_list = ['10044', '10274', '10392']\ninput_dir = '../input/hubmap-organ-segmentation/train_images'\noutput_dir = '.'\ndf['id'] = df['id'].astype(str)\ndf","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","execution":{"iopub.status.busy":"2022-06-22T19:58:14.121183Z","iopub.execute_input":"2022-06-22T19:58:14.122713Z","iopub.status.idle":"2022-06-22T19:58:14.3006Z","shell.execute_reply.started":"2022-06-22T19:58:14.122668Z","shell.execute_reply":"2022-06-22T19:58:14.299086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"TBD","metadata":{}},{"cell_type":"code","source":"def resize_im(im_name, scale_percent):\n    image_path = os.path.join(input_dir, im_name+'.tiff')\n    im_read = tifffile.imread(image_path)\n    width = int(im_read.shape[1] * scale_percent / 100)\n    height = int(im_read.shape[0] * scale_percent / 100)\n    dim = (width, height)\n    print('File name: {}, original size: {}, resized to: {}'.format(im_name, (im_read.shape[0], im_read.shape[1]), (width, height)))\n    resized = cv2.resize(im_read, dim, interpolation=cv2.INTER_AREA)\n    image_path = os.path.join(output_dir, ('r_' + im_name+'.tiff'))\n    tifffile.imwrite(image_path, resized)    ","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:14.303448Z","iopub.execute_input":"2022-06-22T19:58:14.303913Z","iopub.status.idle":"2022-06-22T19:58:14.314198Z","shell.execute_reply.started":"2022-06-22T19:58:14.30387Z","shell.execute_reply":"2022-06-22T19:58:14.312738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Resizing results:","metadata":{}},{"cell_type":"code","source":"for im in image_list:\n    resize_im(im, 5)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:14.38435Z","iopub.execute_input":"2022-06-22T19:58:14.384793Z","iopub.status.idle":"2022-06-22T19:58:15.944507Z","shell.execute_reply.started":"2022-06-22T19:58:14.384758Z","shell.execute_reply":"2022-06-22T19:58:15.943251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's do the same with masks. I will decode relevant masks from the train.csv file and then resize and save it in separate file.\n\nThe function for RLE encoding:","metadata":{}},{"cell_type":"code","source":"def rle2mask(mask_rle, shape):\n    '''\n    mask_rle: run-length as string formated (start length)\n    shape: (width,height) of array to return\n    Returns numpy array, 1 - mask, 0 - background\n    '''\n    s = mask_rle.split()\n    starts, lengths = [np.asarray(x, dtype=int) for x in (s[0:][::2], s[1:][::2])]\n    starts -= 1\n    ends = starts + lengths\n    img = np.zeros(shape[0]*shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        img[lo:hi] = 1\n    return img.reshape(shape).T","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:15.946807Z","iopub.execute_input":"2022-06-22T19:58:15.948075Z","iopub.status.idle":"2022-06-22T19:58:15.956777Z","shell.execute_reply.started":"2022-06-22T19:58:15.948025Z","shell.execute_reply":"2022-06-22T19:58:15.955529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here is the function, which read RLE-mask from the DataFrame, resize it with some scale (in percent) and store it in the folder /output","metadata":{}},{"cell_type":"code","source":"def resize_mask(im_name, scale_percent):    \n    im_read = tifffile.imread(os.path.join(input_dir, im_name +'.tiff'))\n    mask_rle = df[df[\"id\"] == im_name][\"rle\"].values[0]\n    mask = rle2mask(df[df[\"id\"] == im_name][\"rle\"].values[0], (im_read.shape[1], im_read.shape[0]))*255\n    width = int(im_read.shape[1] * scale_percent / 100)\n    height = int(im_read.shape[0] * scale_percent / 100)\n    dim = (width, height)\n    print('File name: {}, original size: {}, resized to: {}'.format(im_name, (im_read.shape[0], im_read.shape[1]), (width, height)))\n    resized = cv2.resize(mask, dim, interpolation=cv2.INTER_AREA)\n    image_path = os.path.join(output_dir, (im_name + '.tiff'))\n    tifffile.imwrite(image_path, resized)    ","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:15.958493Z","iopub.execute_input":"2022-06-22T19:58:15.959473Z","iopub.status.idle":"2022-06-22T19:58:15.972385Z","shell.execute_reply.started":"2022-06-22T19:58:15.959427Z","shell.execute_reply":"2022-06-22T19:58:15.970973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Resizing results:","metadata":{}},{"cell_type":"code","source":"for im in image_list:\n    print(im)\n    resize_mask(im, 5)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:15.975922Z","iopub.execute_input":"2022-06-22T19:58:15.976948Z","iopub.status.idle":"2022-06-22T19:58:16.140305Z","shell.execute_reply.started":"2022-06-22T19:58:15.976896Z","shell.execute_reply":"2022-06-22T19:58:16.138905Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All resized files:","metadata":{}},{"cell_type":"code","source":"os.listdir(output_dir)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:16.142494Z","iopub.execute_input":"2022-06-22T19:58:16.143291Z","iopub.status.idle":"2022-06-22T19:58:16.152287Z","shell.execute_reply.started":"2022-06-22T19:58:16.143241Z","shell.execute_reply":"2022-06-22T19:58:16.150854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='section2'></a>\n## Plotting samples of training images\n\nAnother step is plotting these files with relevant masks. Now, when images are resized it’s not a problem.\nThe function for plotting:","metadata":{}},{"cell_type":"code","source":"def show_image(image_id):\n    fig, ax = plt.subplots(nrows=2, ncols=1, figsize=(16, 32))\n    image_path = os.path.join(output_dir, '{}.tiff'.format(image_id))\n    mask_path = os.path.join(output_dir, 'r_{}.tiff'.format(image_id))    \n    image = tifffile.imread(image_path)\n    mask = tifffile.imread(mask_path)\n    if len(mask.shape) == 2: hybr = image[:, :] + mask[:, :]/2\n    else: hybr = image[:, :] + mask[:,: , 0]/2\n    ax[0].imshow(image)\n    ax[0].axis('off')\n    ax[0].set_title('Real Image')\n    ax[1].imshow(hybr)\n    ax[1].axis('off')\n    ax[1].set_title('Masks')\n    plt.show()    ","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:16.154346Z","iopub.execute_input":"2022-06-22T19:58:16.154877Z","iopub.status.idle":"2022-06-22T19:58:16.167995Z","shell.execute_reply.started":"2022-06-22T19:58:16.154833Z","shell.execute_reply":"2022-06-22T19:58:16.16682Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\nshow_image(image_list[0])","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:16.16946Z","iopub.execute_input":"2022-06-22T19:58:16.169943Z","iopub.status.idle":"2022-06-22T19:58:16.624175Z","shell.execute_reply.started":"2022-06-22T19:58:16.169908Z","shell.execute_reply":"2022-06-22T19:58:16.622867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\nshow_image(image_list[1])","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:16.625895Z","iopub.execute_input":"2022-06-22T19:58:16.626196Z","iopub.status.idle":"2022-06-22T19:58:17.071592Z","shell.execute_reply.started":"2022-06-22T19:58:16.626169Z","shell.execute_reply":"2022-06-22T19:58:17.069936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\nshow_image(image_list[2])","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.073317Z","iopub.execute_input":"2022-06-22T19:58:17.073693Z","iopub.status.idle":"2022-06-22T19:58:17.535629Z","shell.execute_reply.started":"2022-06-22T19:58:17.07366Z","shell.execute_reply":"2022-06-22T19:58:17.534458Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='section3'></a>\n## Image tiling\n\nFor the beginning I will split **'10044.tiff'** with original size for tiles of size 1024x1024 and store all files into the folder **split**:\n-\tImages will be stored in the folder **split/images/**\n-\tMask-files will be stored in the folder **split/masks/**\n","metadata":{}},{"cell_type":"code","source":"os.makedirs('../working/split/images', exist_ok = True)\nos.makedirs('../working/split/masks', exist_ok = True)\nim_name = '10044.tiff'\nimage_path = os.path.join(input_dir, im_name)\ndf = pd.read_csv('../input/hubmap-organ-segmentation/train.csv')\ndf['id'] = df['id'].astype(str)\nsplit_size = 1024\nim = tifffile.imread(os.path.join(input_dir, im_name))\nmask_rle = df[df[\"id\"] == im_name[:-5]][\"rle\"].values[0]\nmask = rle2mask(df[df[\"id\"] == im_name[:-5]][\"rle\"].values[0], (im.shape[1], im.shape[0]))*255\nfor r in range(0, im.shape[0], split_size):\n    for c in range(0, im.shape[1], 1024):\n        im_tile = im[r: r + split_size, c: c + split_size]\n        mask_tile = mask[r: r + split_size, c: c + split_size]\n        # here I filter images with 0-mask and white borders around.\n        if (np.sum(mask_tile)==0):\n            if ((2 * split_size <= r <= (im.shape[0] - 2 * split_size)) and \\\n                (2 * split_size <= c <= (im.shape[1] - 2 * split_size))):\n                tifffile.imwrite(f\"split/images/img{r}_{c}.png\", im_tile)\n                tifffile.imwrite(f\"split/masks/img{r}_{c}.png\", mask_tile)\n        else:\n            tifffile.imwrite(f\"split/images/img{r}_{c}.png\", im_tile)\n            tifffile.imwrite(f\"split/masks/img{r}_{c}.png\", mask_tile)\n            ","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.540342Z","iopub.execute_input":"2022-06-22T19:58:17.540954Z","iopub.status.idle":"2022-06-22T19:58:17.805272Z","shell.execute_reply.started":"2022-06-22T19:58:17.540916Z","shell.execute_reply":"2022-06-22T19:58:17.803956Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As a result I received a set of images and masks (label) for a model training. Let’s count just for information: ","metadata":{}},{"cell_type":"code","source":"len(os.listdir('split/images'))","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.806998Z","iopub.execute_input":"2022-06-22T19:58:17.807453Z","iopub.status.idle":"2022-06-22T19:58:17.815574Z","shell.execute_reply.started":"2022-06-22T19:58:17.807407Z","shell.execute_reply":"2022-06-22T19:58:17.814164Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And some statistics. Let's calculate the areas of the masks at each images, but for the beginning files with 0-mask. For comfortable calculations I will use Pandas DataFrame where index is file name:","metadata":{}},{"cell_type":"code","source":"mask_list = os.listdir('split/masks')\ndf=pd.DataFrame(index=mask_list)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.817623Z","iopub.execute_input":"2022-06-22T19:58:17.81817Z","iopub.status.idle":"2022-06-22T19:58:17.827657Z","shell.execute_reply.started":"2022-06-22T19:58:17.818123Z","shell.execute_reply":"2022-06-22T19:58:17.826644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Area calculation function:","metadata":{}},{"cell_type":"code","source":"def area_calc(image_id):\n    mask_path = os.path.join(mask_dir, '{}'.format(image_id))\n    mask = cv2.imread(mask_path)\n    return int(np.count_nonzero(mask) / 3)\n","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.829489Z","iopub.execute_input":"2022-06-22T19:58:17.83067Z","iopub.status.idle":"2022-06-22T19:58:17.84016Z","shell.execute_reply.started":"2022-06-22T19:58:17.830623Z","shell.execute_reply":"2022-06-22T19:58:17.839365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Calculate and write mask areas (sum) per image to the DataFrame:","metadata":{}},{"cell_type":"code","source":"mask_dir = 'split/masks'\nmask_areas=[]\nfor msk in mask_list:\n    mask_areas.append(area_calc(msk))\ndf['area'] = mask_areas\n","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.841249Z","iopub.execute_input":"2022-06-22T19:58:17.842793Z","iopub.status.idle":"2022-06-22T19:58:17.878099Z","shell.execute_reply.started":"2022-06-22T19:58:17.84271Z","shell.execute_reply":"2022-06-22T19:58:17.877277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Total images:', len(df))\nprint('Non-zero images:', len(df[df['area']!=0]))","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.879629Z","iopub.execute_input":"2022-06-22T19:58:17.880212Z","iopub.status.idle":"2022-06-22T19:58:17.887933Z","shell.execute_reply.started":"2022-06-22T19:58:17.880155Z","shell.execute_reply":"2022-06-22T19:58:17.886487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='section4'></a>\n## Mask-area per image distribution and some statistics\n\nNon-zero Image distribution: ","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\nfig, ax = plt.subplots(1,1,figsize=(18,8))\nax.hist(df[df['area']!=0].values, bins=50, color='deeppink', edgecolor='black')  # `density=False` would make counts\nax.set_title('Non-zero Image destribution. Image File: {}     Total images: {}'.format(im_name, \n                                                                                       str(len(df[df['area']!=0]))), \n             fontsize=20)\nax.set_ylabel('Quantity', fontsize=16)\nax.set_xlabel('Area(pixels)', fontsize=16);\nax.grid()\n","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:17.890103Z","iopub.execute_input":"2022-06-22T19:58:17.890559Z","iopub.status.idle":"2022-06-22T19:58:18.274054Z","shell.execute_reply.started":"2022-06-22T19:58:17.890527Z","shell.execute_reply":"2022-06-22T19:58:18.272583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sorted = df[df['area']!=0].sort_values(by=['area'])\nsmallest_list = df_sorted.head(5)['area'].index\nlargest_list = df_sorted.tail(5)['area'].index\nzero_list = df[df['area']==0].head(5)['area'].index\nprint('Smallest:', list(smallest_list))\nprint('Largest:', list(largest_list))\nprint('Zero_list:', list(zero_list))\n","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:18.275702Z","iopub.execute_input":"2022-06-22T19:58:18.276023Z","iopub.status.idle":"2022-06-22T19:58:18.288843Z","shell.execute_reply.started":"2022-06-22T19:58:18.275993Z","shell.execute_reply":"2022-06-22T19:58:18.286995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's take a look at the pictures with the largest and the smallest areas, bur first of all I will sort non-zero mask areas:","metadata":{}},{"cell_type":"markdown","source":"<a id='section5'></a>\n## Plotting images with the smallest and the largest mask ares\n\nA bit modified function for small image plotting:","metadata":{}},{"cell_type":"code","source":"def show_image(image_name):\n    fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(32, 16))\n    image_path = os.path.join('../working/split/images', image_name)\n    mask_path = os.path.join('../working/split/masks', image_name)\n    image = tifffile.imread(image_path)\n    mask = tifffile.imread(mask_path)\n    if len(mask.shape)==2:    \n        hybr = image[:, :, 0] + mask[:, :]/2\n    else:\n        hybr = image[:, :, 0] + mask[:,: , 0]/2\n    ax[0].imshow(image)\n    ax[0].axis('off')\n    ax[0].set_title('Real Image')\n    ax[1].imshow(hybr)\n    ax[1].axis('off')\n    ax[1].set_title('Masks')\n    plt.show()\n    ","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:18.290996Z","iopub.execute_input":"2022-06-22T19:58:18.291455Z","iopub.status.idle":"2022-06-22T19:58:18.302678Z","shell.execute_reply.started":"2022-06-22T19:58:18.291418Z","shell.execute_reply":"2022-06-22T19:58:18.301359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plot 5 images with the smallest mask areas","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\nfor file in smallest_list:\n    show_image(file)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:18.304885Z","iopub.execute_input":"2022-06-22T19:58:18.305322Z","iopub.status.idle":"2022-06-22T19:58:24.589269Z","shell.execute_reply.started":"2022-06-22T19:58:18.305276Z","shell.execute_reply":"2022-06-22T19:58:24.587917Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plot 5 images with the largest mask areas","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\nfor file in largest_list:\n    show_image(file)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:24.590979Z","iopub.execute_input":"2022-06-22T19:58:24.59182Z","iopub.status.idle":"2022-06-22T19:58:31.174284Z","shell.execute_reply.started":"2022-06-22T19:58:24.591779Z","shell.execute_reply":"2022-06-22T19:58:31.173051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"And some Zero-mask images:","metadata":{}},{"cell_type":"code","source":"%matplotlib inline\nfor file in zero_list:\n    show_image(file)","metadata":{"execution":{"iopub.status.busy":"2022-06-22T19:58:31.175905Z","iopub.execute_input":"2022-06-22T19:58:31.176221Z","iopub.status.idle":"2022-06-22T19:58:31.182971Z","shell.execute_reply.started":"2022-06-22T19:58:31.176192Z","shell.execute_reply":"2022-06-22T19:58:31.181725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a id='section6'></a>\n## Conclusion\n\nHere is brief overview of the images provided by **kaggle** for this semantic segmentation competition.\nLet's see what UNet Neural Network monsters you guys are going to come up with!\n\n[Jump on top](#section0)","metadata":{}},{"cell_type":"markdown","source":"> writing code under pressure is so much adrenaline.. I won't sleep for a week!","metadata":{}}]}