{"cells":[{"metadata":{},"cell_type":"markdown","source":"# HuBMAP - Exploratory Data Analysis (EDA)\nThis notebook provides brief exploratory data analysis for the new HuBMAP data set. The full kidney images in the training dataset are visualized with the glomeruli FTUs highlighted. A brief analysis of the shape of the glomerulis follows.\n\n# References\nThe following references were used in this notebook.\n- Reading images: https://www.kaggle.com/ihelon/hubmap-exploratory-data-analysis\n- RLE encoding: https://www.kaggle.com/paulorzp/rle-functions-run-lenght-encode-decode"},{"metadata":{},"cell_type":"markdown","source":"# Imports"},{"metadata":{"trusted":true},"cell_type":"code","source":"import cv2\nimport datetime\nimport gc\nimport glob\nimport math\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport os\nimport pandas as pd\nimport skimage.morphology\nimport sys\nimport tensorflow as tf\nimport tifffile","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"base_path = '../input/hubmap-kidney-segmentation'\n\nplot_full_image = True\n\n# Number of glomeruli to display for each image\nnum_glom_display = 5\n\n# Number of glomberuli to save as tiff files.\nnum_glom_save = 5\n\nglob_scale = 0.25","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Utility Functions"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def rle_to_image(rle_mask, image_shape):\n    \"\"\"\n    Converts an rle string to an image represented as a numpy array.\n    Reference: https://www.kaggle.com/paulorzp/rle-functions-run-lenght-encode-decode\n\n    :param rle_mask: string with rle mask.\n    :param image_shape: (width, height) of array to return\n    :return: Image as a numpy array. 1 = mask, 0 = background.\n    \"\"\"\n\n    # Processing\n    s = rle_mask.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    image = np.zeros(image_shape[0] * image_shape[1], dtype=np.uint8)\n    for lo, hi in zip(starts, ends):\n        image[lo:hi] = 1\n\n    return image.reshape(image_shape).T\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# File Structure\nThe files in the root of the dataset are shown below. The dataset consists of 2 directories that contain training and test images and 3 csv-files with additional information about the images.\n## Directory Contents"},{"metadata":{"trusted":true},"cell_type":"code","source":"print('\\n'.join(os.listdir(base_path)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training Images\nThe train directory contains 15 images for training."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_files = sorted(glob.glob(os.path.join(base_path, 'train/*.tiff')))\nprint(f'Number of training images: {len(train_files)}')\nprint('\\n'.join(train_files))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Test Images\nThe test directory contains 5 images for testing."},{"metadata":{"trusted":true},"cell_type":"code","source":"test_files = sorted(glob.glob(os.path.join(base_path, 'test/*.tiff')))\nprint(f'Number of test images: {len(test_files)}')\nprint('\\n'.join(test_files))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Train.csv\nThe masks indicating a glomeruli FTUs are stored in rle format in the train.csv for each training image id."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(os.path.join(base_path, 'train.csv'))\ndisplay(df_train)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Sample_Sumbission.csv\nThe sample_submission.csv files shows the format of the submissions files consisting of the test image id and an rle encoded masks."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_submission = pd.read_csv(os.path.join(base_path,'sample_submission.csv'))\ndisplay(df_submission)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Patient Data\nHuBMAP-20-dataset_information.csv contains additional information about each image such as image size and anonymized patient data."},{"metadata":{"trusted":true},"cell_type":"code","source":"df_info = pd.read_csv(os.path.join(base_path,'HuBMAP-20-dataset_information.csv'))\ndisplay(df_info)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.options.display.float_format = '{:,.1f}'.format\ndf_info.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Training Image Analysis"},{"metadata":{},"cell_type":"markdown","source":"## Width and Height Distribution\nThe training images do not have consistent dimensions. This has to be corrected when loading the images. They have on of the following shapes:\n- [height, width, channel]\n- [channel, height, width]\n- [1, 1, channel, height, width]"},{"metadata":{"trusted":true},"cell_type":"code","source":"for f in train_files + test_files:\n    image = tifffile.imread(f)\n    print(f'Image {f} shape: {image.shape}', flush=True)\n    del image\n    gc.collect()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The size of the images varies greatly as well. "},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.scatter(df_info['width_pixels'], df_info['height_pixels'])\nplt.title('Image Height and Width')\nplt.xlabel('Width')\nplt.ylabel('Height')\nplt.xlim(0, df_info['width_pixels'].max() * 1.1)\nplt.ylim(0, df_info['height_pixels'].max() * 1.1)\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Image Utilitity Functions"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def overlay_image_mask(image, mask, mask_color=(0,255,0), alpha=1.0):\n    im_f= image.astype(np.float32)\n#     if mask.ndim == 2:\n#         mask = np.expand_dims(mask,-1)        \n    mask_col = np.expand_dims(np.array(mask_color)/255.0, axis=(0,1))\n    return (im_f + alpha * mask * (np.mean(0.8 * im_f + 0.2 * 255, axis=2, keepdims=True) * mask_col - im_f)).astype(np.uint8)\n\n\ndef overlay_image_mask_original(image, mask, mask_color=(0,255,0), alpha=1.0):\n    return  np.concatenate((image, overlay_image_mask(image, mask)), axis=1)\n\ndef get_image_id(image_file):\n    return os.path.splitext(os.path.split(image_file)[1])[0]\n\n\ndef read_image(image_file, scale=1.0):\n    image = tifffile.imread(image_file).squeeze()\n    if image.shape[0] == 3:\n        image = np.transpose(image, (1,2,0))\n    \n    orig_shape = image.shape\n    if scale != 1.0:\n        image = cv2.resize(image, (0,0), fx=scale, fy=scale)\n    return image, orig_shape\n\n\ndef read_mask(image_file, image_shape, scale=1.0):\n    image_id = get_image_id(image_file)\n    train_info = df_train.loc[df_train['id'] == image_id]\n    rle = train_info['encoding'].values[0] if len(train_info) > 0 else None\n    if rle is not None:\n        mask = rle_to_image(rle, (image_shape[1], image_shape[0]))\n        if scale != 1.0:\n            mask = cv2.resize(mask, (0,0), fx=scale, fy=scale)\n        return np.expand_dims(mask,-1)\n    else:\n        return None        \n\n    \ndef read_image_mask(image_file, scale=1.0):\n    image, image_shape = read_image(image_file, scale)\n    mask = read_mask(image_file, image_shape, scale)\n    return image, mask\n\n\ndef get_tile(image, mask, x, y, tile_size, scale=1.0):\n    x = round(x * scale)\n    y = round(y * scale)\n    size = int(round(tile_size / 2 * scale))\n    image_s = image[y-size:y+size, x-size:x+size, :] \n    mask_s = mask[y-size:y+size, x-size:x+size, :]\n    return image_s, mask_s\n\n\ndef get_particles(mask, scale=1.0):\n    num, labels, stats, centroids = cv2.connectedComponentsWithStats(mask)\n    df_particles = pd.DataFrame(dict(zip(['x','y','left','top','width','height','area'],\n                               [(centroids[1:,0]) / scale,\n                                (centroids[1:,1]) / scale,\n                                (stats[1:,cv2.CC_STAT_LEFT]) / scale,\n                                (stats[1:,cv2.CC_STAT_TOP]) / scale,\n                                (stats[1:,cv2.CC_STAT_WIDTH]) / scale,\n                                (stats[1:,cv2.CC_STAT_HEIGHT]) / scale,\n                                (stats[1:,cv2.CC_STAT_AREA]) / (scale * scale)])))\n    df_particles.sort_values(['x','y'], inplace=True, ignore_index=True)\n    df_particles['no'] = range(len(df_particles))\n    return df_particles\n\n\ndef analyze_image(image_file):\n    image_id = get_image_id(image_file)\n    image, image_shape = read_image(image_file, glob_scale)\n    mask = read_mask(image_file, image_shape, glob_scale)\n\n    mask_full = read_mask(image_file, image_shape, scale=1.0)\n    df_glom = get_particles(mask_full, scale=1.0)\n    df_glom['id'] = image_id\n    del mask_full\n    gc.collect()\n    \n    info = df_info[df_info['image_file'] == f'{image_id}.tiff']\n    print(f'Image ID:        {image_id:}')\n    print(f'Image Size:      {info[\"width_pixels\"].values[0]} x {info[\"height_pixels\"].values[0]}')\n    print(f'Patient No:      {info[\"patient_number\"].values[0]}')\n    print(f'Sex:             {info[\"sex\"].values[0]}')\n    print(f'Age:             {info[\"age\"].values[0]}')\n    print(f'Race:            {info[\"race\"].values[0]}')\n    print(f'Height:          {info[\"height_centimeters\"].values[0]} cm')\n    print(f'Weight:          {info[\"weight_kilograms\"].values[0]} kg')\n    print(f'BMI:             {info[\"bmi_kg/m^2\"].values[0]} kg/m^2')\n    print(f'Laterality:      {info[\"laterality\"].values[0]}')\n    print(f'Percent Cortex:  {info[\"percent_cortex\"].values[0]} %')\n    print(f'Percent Medulla: {info[\"percent_medulla\"].values[0]} %')\n    \n    # Plot full image\n    if plot_full_image:\n        scale = 0.1\n        image_small = cv2.resize(image, (0,0), fx=scale, fy=scale)\n        mask_small = cv2.resize(mask, (0,0), fx=scale, fy=scale)\n        mask_small = np.expand_dims(mask_small,-1) \n    \n        plt.figure(figsize=(16, 16))\n        plt.imshow(overlay_image_mask(image_small, mask_small))\n        plt.axis('off')\n\n    # Plot glomeruli images\n    fig_cols = 5\n    fig_rows = int(math.ceil(num_glom_display/fig_cols))\n    plt.figure(figsize=(4 * fig_cols, 4 * fig_rows))\n    if num_glom_save > 0 and not os.path.exists(image_id):\n        os.mkdir(image_id)\n    for i in range(min(max(num_glom_display, num_glom_save), len(df_glom))):\n        image_s, mask_s = get_tile(image,mask, df_glom['x'][i], df_glom['y'][i], 1000, scale=glob_scale)\n        ovl = overlay_image_mask(image_s, mask_s)\n        if i < num_glom_display:\n            plt.subplot(fig_rows, fig_cols, i+1)\n            plt.imshow(ovl)\n            plt.axis('off')\n        if i < num_glom_save:\n            cv2.imwrite(f'{image_id}_{i:03}.png', cv2.cvtColor(ovl, cv2.COLOR_RGB2BGR))    \n    \n    del image, mask\n    gc.collect()\n    return df_glom\n\n\ndef plot_glom(df, image_id, glom_no):\n    image, mask = read_image_mask(os.path.join(base_path, f'train/{image_id}.tiff'), scale=glob_scale)\n    glom = df.loc[(df['id'] == image_id) & (df['no'] == glom_no)]\n    im, ma = get_tile(image, mask, glom['x'].iloc[0], glom['y'].iloc[0], 1000, scale=glob_scale)\n    del image, mask\n    gc.collect()\n    plt.figure(figsize=(16,8))\n    plt.imshow(overlay_image_mask_original(im, ma))\n    plt.title(f'Image: {image_id}, Glomeruli No: {glom_no}, Area: {glom[\"area\"].iloc[0]}')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training Images With Glomerulis"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = pd.DataFrame()\ndf_glom = df_glom.append(analyze_image(train_files[0]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[1]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[2]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[3]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[4]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[5]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[6]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[7]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[8]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[9]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[10]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[11]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[12]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[13]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom = df_glom.append(analyze_image(train_files[14]), ignore_index=True)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Glomerulis\n## Basic Statistics"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom.to_csv('glomeruli.csv')\ndisplay(df_glom)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Glomerulis Per Image"},{"metadata":{"trusted":true},"cell_type":"code","source":"g = df_glom.groupby('id')\nplt.bar(g.size().index, g.size().values)\nplt.title('Number of Glomerulis in Image')\nplt.xticks(rotation=90)\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Glomeruli Width, Height and Area Distribution"},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20,5))\nplt.subplot(1,3,1)\nplt.hist(df_glom['width'], bins=40, density=True)\nplt.title('Width Distribution')\nplt.grid()\nplt.subplot(1,3,2)\nplt.hist(df_glom['height'], bins=40, density=True)\nplt.title('Height Distribution')\nplt.grid()\nplt.subplot(1,3,3)\nplt.hist(df_glom['area'], bins=40, density=True)\nplt.title('Area Distribution')\nplt.grid()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Glomerulis by Size"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_glom.sort_values('area', inplace=True)\ndf_glom","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5 Smallest Glomerulis"},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(5):\n    plot_glom(df_glom, df_glom['id'].iloc[i], df_glom['no'].iloc[i])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## 5 Largest Glomerulis"},{"metadata":{"trusted":true},"cell_type":"code","source":"for i in range(len(df_glom)-5, len(df_glom)):\n    plot_glom(df_glom, df_glom['id'].iloc[i], df_glom['no'].iloc[i])","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}