{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Purpose\nThis notebook creates a new cleaner data set from the images provided.  For each sample 4 new images will be created for each cell segmented out of the original 4 images.\n\nFor example, if ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0 has 3 cells:\n```\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_red.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_green.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_blue.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_yellow.png\n```\nWill become:\n```\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_red_1.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_red_2.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_red_3.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_green_1.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_green_2.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_green_3.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_blue_1.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_blue_2.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_blue_3.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_yellow_1.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_yellow_2.png\n./ca035f36-bbc9-11e8-b2bc-ac1f6b6435d0_yellow_3.png\n```\n\nReferences:\n\nhttps://www.kaggle.com/christopherworley/human-protein-atlas-segmentation\n    \nhttps://github.com/CellProfiling/HPA-Cell-Segmentation\n"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"\nimport pandas as pd \nimport numpy as np\n\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\n\nfrom PIL import Image  \nfrom IPython import display \n\nimport os\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install 'git+https://github.com/haoxusci/pytorch_zoo@master#egg=pytorch_zoo'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install https://github.com/CellProfiling/HPA-Cell-Segmentation/archive/master.zip","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import hpacellseg.cellsegmentator as cellsegmentator\nfrom hpacellseg.utils import label_cell, label_nuclei","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Parameters"},{"metadata":{"trusted":true},"cell_type":"code","source":"ROOT_PATH = '/kaggle/input/hpa-single-cell-image-classification/'\n\nCHANNELS = np.array(['red', 'green', 'blue', 'yellow'])\nCHANNEL_DESCRIPTIONS = np.array(['Microtubule Channels', 'Protein of Interest', 'Nuclei Channels', 'Endoplasmic Reticulum'])\n\nCHANNEL_RED = 0    # microtubule channels\nCHANNEL_GREEN = 1  # protein of interest\nCHANNEL_BLUE = 2   # nuclei channels\nCHANNEL_YELLOW = 3 # endoplasmic reticulum \n\nCHANNEL_SIZE = len(CHANNELS)\nSAMPLE_SIZE = -1\n\n# image size for single cell\nCELL_IMAGE_SIZE = 2048\n\n# Plot or do not plot images.  Plotting will slow things down\nPLOT_OUTPUT = False\n\n# data directory\nOUTPUT_DIR = \"./data\"\nif os.path.exists(OUTPUT_DIR) == False:\n    os.mkdir(OUTPUT_DIR)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Functions"},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# Get array of all images for given sample id\n#\ndef get_images(id):\n    images = list()\n    for channel in CHANNELS:\n        path = ROOT_PATH + 'train/{}_{}.png'.format(id, channel)\n        image = Image.open(path) \n        images.append( np.asarray(image) )\n    return images","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# Get single image that blends all RGBY into RGB\n#\ndef get_blended_image(id): \n    # get rgby images for sample\n    images = get_images(id)\n\n    # blend rgby images into single array\n    blended_array = np.stack((\n            np.maximum(images[0], images[3]),\n            np.maximum(images[1], images[3]),\n            images[2]\n        ), 2)\n\n    # Create PIL Image\n    blended_image = Image.fromarray( np.uint8(blended_array) )\n    return blended_image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# Crop given image and create new image with cropped image centered\n#\ndef center_mass(image, top, right, bottom, left):\n    centered_image = np.zeros((CELL_IMAGE_SIZE, CELL_IMAGE_SIZE))\n    x1, y1 = int((centered_image.shape[0] / 2) - ((right-left)/2)), int((centered_image.shape[1] / 2) - ((bottom-top)/2))\n    x2, y2 = x1 + (right - left), y1 + (bottom - top)\n    centered_image[y1:y2, x1:x2] = image[top:bottom, left:right]\n    return np.uint8(centered_image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# Determine boundries of object in given array\n#\ndef get_cell_bounds(image_array):\n    top = -1\n    right = -1\n    bottom = -1\n    left = -1\n    \n    # find upper and lower bounds\n    for index in range(image_array.shape[0]):\n        is_empty = np.sum( np.unique( image_array[index] ) ) == 0\n        if top == -1 and is_empty == False:\n            top = index\n        elif is_empty == False:\n            bottom = index\n            \n    # find left and right bounds\n    for index in range(image_array.shape[1]):\n        is_empty = np.sum( np.unique( image_array[:,index] ) ) == 0\n        if left == -1 and is_empty == False:\n            left = index\n        elif is_empty == False:\n            right = index\n    \n    return (top, right, bottom, left)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#\n# Create mask for given images\n#\ndef create_cell_mask(images):\n        \n    # Segment nuclie\n    nuc_segmentations = segmentator.pred_nuclei([np.asarray( images[CHANNEL_BLUE] )])\n\n    # For full cells\n    cell_segmentations = cell_segmentations = segmentator.pred_cells([\n        [np.asarray( images[CHANNEL_RED] )],\n        [np.asarray( images[CHANNEL_YELLOW] )],\n        [np.asarray( images[CHANNEL_BLUE] )]\n    ])\n\n    nuclei_mask, cell_mask = label_cell(nuc_segmentations[0], cell_segmentations[0])\n    \n    return cell_mask\n    \n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Read Trainning Data"},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train = pd.read_csv(ROOT_PATH + 'train.csv')\nprint(\"Trainning data length: {}\".format(len(df_train)))\ndf_train.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# if sample size is set then reduce trainning set accordingly\nif SAMPLE_SIZE > -1:\n    df_train = df_train.sample(SAMPLE_SIZE)\n    df_train.reset_index(inplace=True);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"NUC_MODEL = \"./nuclei-model.pth\"\nCELL_MODEL = \"./cell-model.pth\"\nsegmentator = cellsegmentator.CellSegmentator(\n    NUC_MODEL,\n    CELL_MODEL,\n    scale_factor=0.25,\n    device=\"cuda\",\n    padding=False,\n    multi_channel_model=True,\n)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Main Loop Through Sampled Data\n\n1. Loop through samples\n2. Loop through all cells identified\n3. Create new image with single cell centered\n4. Write new cell image out"},{"metadata":{"trusted":true},"cell_type":"code","source":"\nCOLUMN_COUNT = CHANNEL_SIZE + 1\nROW_COUNT = 2\n\nsample_count = 1\n\n# setup figure to plot blened image and centered images for each channel\nfig = plt.figure(figsize=(40,10))\n\nfor sample_index, df_sample in df_train.iterrows():\n    images = get_images(df_sample['ID'])\n    \n    # Blended image\n    print(\"\\rGenerating blended image (count={}, current_id={})                  \".format(sample_count, df_sample['ID']), end=\"\")\n    blended_image = get_blended_image(df_sample['ID'])\n    blended_array = np.asarray(blended_image)\n\n    # Plot blended image\n    if PLOT_OUTPUT == True:\n        ax = fig.add_subplot(1, COLUMN_COUNT, 1)\n        ax.set_title(\"Blended Cell\")\n        plt.imshow(np.asarray(blended_image))\n\n        for channel_index in range(len(CHANNELS)):\n            ax = fig.add_subplot(ROW_COUNT, COLUMN_COUNT, channel_index + 2)\n            ax.set_title(\"Original {} {}\".format(sample_count, CHANNEL_DESCRIPTIONS[channel_index]))\n            plt.imshow(images[channel_index], cmap=plt.get_cmap('bone'))\n        \n    # Cell mask\n    print(\"\\rCreating single cell data set (count={}, current_id={})              \".format(sample_count, df_sample['ID']), end=\"\")\n    cell_mask = create_cell_mask(images)\n    \n    # Get unique vector of segment numbers\n    numbers = np.unique(cell_mask)\n    numbers = np.delete(numbers, [0])\n\n    # Get cell bounds \n    (top, right, bottom, left) = get_cell_bounds(blended_array)\n    \n    for number in numbers:\n        print(\"\\rIsolating single cell (count={}, cell_number={}, current_id={})   \".format(sample_count, number, df_sample['ID']), end=\"\")\n        # Isolate the cell within the cell mask\n        isolated_mask = np.where(cell_mask == number, True, False)\n        isolated_mask_multi_channel = np.stack((isolated_mask, isolated_mask, isolated_mask), axis=2)\n        \n        # Get boundries of the isolated blended cell\n        blended_isolated_image = np.where(isolated_mask_multi_channel == True, blended_array, 0)\n        (top, right, bottom, left) = get_cell_bounds(blended_isolated_image)\n        \n        # create new single cell images for each channel\n        index = 7\n        for channel_index in range(len(CHANNELS)):\n            isolated_image = np.where(isolated_mask == True, images[channel_index], 0)\n            centered_image = center_mass(isolated_image, top, right, bottom, left)\n            Image.fromarray(centered_image).save(\"./{}/{}_{}_{}.png\".format(OUTPUT_DIR, df_sample['ID'], CHANNELS[channel_index], number), \"PNG\")\n\n            if PLOT_OUTPUT == True:\n                ax = fig.add_subplot(ROW_COUNT, COLUMN_COUNT, index)\n                ax.set_title(\"Clean Cell {} {}\".format(number, CHANNEL_DESCRIPTIONS[channel_index]))\n                plt.imshow(centered_image, cmap=plt.get_cmap('bone'))\n                index = index + 1\n\n        if PLOT_OUTPUT == True:\n            display.clear_output(wait=True)\n            display.display(fig)\n    print(\"\\rSingle cell data set completed (count={}, current_id={})               \".format(sample_count, df_sample['ID']), end=\"\")\n    sample_count = sample_count + 1\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"","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}