{"cells":[{"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 os\nimport numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport tensorflow as tf\nfrom matplotlib import pyplot as plt\nimport hpacellseg.cellsegmentator as cellsegmentator\nfrom hpacellseg.utils import label_cell, label_nuclei","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **FUNCTIONS**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Input: list of image filters as png\n# Output: list of image filters as np.arrays\ndef image_to_arrays(path):\n    \n    image_arrays = list()\n    for image in path:\n        array = np.asarray(Image.open(image))\n        image_arrays.append(array)\n        \n    return image_arrays","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"> Reference: [Human Protein Atlas - Segmentation](https://www.kaggle.com/christopherworley/human-protein-atlas-segmentation#Functions)"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get single image that blends all RGBY into RGB\n# Introduce the images as arrays. Can use the function above.\n\ndef get_blended_image(images): \n    # get rgby images for sample\n\n    # blend rgby images into single array\n    blended_array = np.stack(images[:-1], 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":"# Introduce list of image filters\n# Returns a processed image ready for the CNN and an encoded label as tensor\ndef image_prep(paths, label):\n\n    img = image_to_arrays(paths)\n    size = np.shape(img[0])[0]\n    img = tf.image.convert_image_dtype(img, dtype=tf.float32)\n    img = tf.reshape(img, (1, size, size, 3))\n    img = tf.image.resize(img, IMG_SIZE)\n\n    label = tf.strings.split(label, sep='|')\n    label = tf.strings.to_number(label, out_type=tf.int32)\n    label = tf.reduce_sum(tf.one_hot(indices=label, depth=19), axis=0)\n    label = tf.reshape(label, (1, 19))\n    \n    return img, label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def apply_augmentation(image, label):\n    aug_img = tf.numpy_function(func=aug_fn, inp=[image], Tout=tf.float32)\n    aug_img.set_shape((IMG_SIZE[0], IMG_SIZE[0], 3))\n    \n    return aug_img, label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def plot_hist(hist):\n    plt.plot(hist.history[\"accuracy\"])\n    plt.plot(hist.history[\"val_accuracy\"])\n    plt.title(\"model accuracy\")\n    plt.ylabel(\"accuracy\")\n    plt.xlabel(\"epoch\")\n    plt.legend([\"train\", \"validation\"], loc=\"upper left\")\n    plt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/hpa-single-cell-image-classification\"\n\ntrain = pd.read_csv(os.path.join(DATA_DIR,'train.csv'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"colours = ['_red.png', '_blue.png', '_yellow.png', '_green.png']\nTRAIN = '../input/hpa-single-cell-image-classification/train'\npaths = [[os.path.join(TRAIN, train.iloc[idx,0])+ colour for colour in colours] for idx in range(len(train))]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# *Data Analisys...*"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's check out the label distribution frequency.\nlabel_counts = []\nfor label in train['Label']:\n    sep = label.split('|')\n    for num in sep:\n        labels.append(int(num))\ncounts = pd.value_counts(labels)\n\n# It's an ugly plot, but I'm trying to save some time here...\nplt.bar(x = counts.index,height=counts)\nplt.xticks(counts.index)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"titles = ['microtubules', 'nuclei', 'endoplasmic reticulum', 'protein of interest']\nfig, axs = plt.subplots(3, 4, figsize =(16,8))\nfor entry in range(3):\n    for channel in range(4):\n        img = plt.imread(paths[entry][channel])\n        axs[entry, channel].imshow(img)        \n        if entry == 0:\n            axs[0, channel].set_title(titles[channel])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Segmentation using [HPA-Cell-Segmentation](https://github.com/CellProfiling/HPA-Cell-Segmentation)**"},{"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\nimage = paths[4]\narrays = image_to_arrays(image)\nnuclei = arrays[1]\ncell = arrays[:-1]\n\n# Nuclei segmentation\nnuc_segmentations = segmentator.pred_nuclei([nuclei])\n\nf, ax = plt.subplots(1, 2, figsize=(16,16))\nax[0].imshow(arrays[1])\nax[0].set_title('Original Nucleis', size=20)\nax[1].imshow(nuc_segmentations[0])\nax[1].set_title('Segmented Nucleis', size=20)\nplt.show()\n\n# Cell segmentation\ninter_step = [[i] for i in image[:-1]]\ncell_segmentations = segmentator.pred_cells(inter_step)\n\nf, ax = plt.subplots(1, 2, figsize=(16,16))\nax[0].imshow(get_blended_image(arrays))\nax[0].set_title('Original Cells', size=20)\nax[1].imshow(cell_segmentations[0])\nax[1].set_title('Segmented Cells', size=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Visualizing the masks**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Nuclei mask\nnuclei_mask = label_nuclei(nuc_segmentations[0])\n# Cell masks\ncell_nuclei_mask, cell_mask = label_cell(nuc_segmentations[0], cell_segmentations[0])\n# Plotting\nf, ax = plt.subplots(1, 3, figsize=(16,16))\nax[0].imshow(nuclei_mask)\nax[0].set_title('Nuclei Mask', size=20)\nax[1].imshow(cell_nuclei_mask)\nax[1].set_title('Cell Nuclei Mask', size=20)\nax[2].imshow(cell_mask)\nax[2].set_title('Cell Mask', size=20)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's check the results of the segmentation"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Let's stack the original image and the segmentation mask, to see how the segmentation worked out\nplt.figure(figsize=(20,20))\nplt.imshow(get_blended_image(arrays))\nplt.imshow(cell_mask, alpha=0.5)\nplt.title('Segmentation results', size=40)\nplt.axis('off')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Cell separation**\n\nThe objective of this project is to label each cell in the image. Therefore each cell in the image must be separated."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Unique vector of cell_mask numbers\nnumbers = set(np.ravel(cell_mask))\nnumbers.remove(0)\n\nfig = plt.figure(figsize=(25,6*len(numbers)/4))\nindex = 1\n\nax = fig.add_subplot(len(numbers)//4+1, 4, index)\nax.set_title(\"Complete Cell Mask\", size=20)\nplt.imshow(cell_mask)\n\nindex += 1\nfor number in numbers:\n    isolated_cell = np.where(cell_mask==number, cell_mask, 0)\n    ax = fig.add_subplot(len(numbers)//4+1, 4, index)\n    ax.set_title(\"Segment {number}\", size=20)\n    plt.imshow(isolated_cell)\n    index += 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now that the segmentation is complete. We should be able to train an image classification model to identify each cell within the image.\n\nThe main problem is that the labels are given for each image, therefore we don't really know which of the cells in the image may represent such label.\nMaybe the CNN is able to understand the pattern given the same label for every cell of the image, although it can lead to high misslabeling."},{"metadata":{},"cell_type":"markdown","source":"# **TRAINING MODEL SETUP**"},{"metadata":{},"cell_type":"markdown","source":"Reference: [HPA: Multi-Label Classification with TF and W&B](https://www.kaggle.com/ayuraj/hpa-multi-label-classification-with-tf-and-w-b)"},{"metadata":{},"cell_type":"markdown","source":"Imports."},{"metadata":{"trusted":true},"cell_type":"code","source":"from tensorflow.keras.applications import EfficientNetB0\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import Model\nimport wandb","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Parameter setting"},{"metadata":{"trusted":true},"cell_type":"code","source":"LABELS= {\n0: \"Nucleoplasm\",\n1: \"Nuclear membrane\",\n2: \"Nucleoli\",\n3: \"Nucleoli fibrillar center\",\n4: \"Nuclear speckles\",\n5: \"Nuclear bodies\",\n6: \"Endoplasmic reticulum\",\n7: \"Golgi apparatus\",\n8: \"Intermediate filaments\",\n9: \"Actin filaments\",\n10: \"Microtubules\",\n11: \"Mitotic spindle\",\n12: \"Centrosome\",\n13: \"Plasma membrane\",\n14: \"Mitochondria\",\n15: \"Aggresome\",\n16: \"Cytosol\",\n17: \"Vesicles and punctate cytosolic patterns\",\n18: \"Negative\"\n}","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We'll use EfficientNetB0 model, which requires an image dimension of (224,224,3).Therefor, we can only pass a 3 filter image... \n#We'll put aside the yellow filter for now.\nIMG_SIZE = [224, 224]\nBATCH_SIZE = 64\nAUTOTUNE = tf.data.experimental.AUTOTUNE\n\ncolours = ['_red.png', '_blue.png', '_green.png']\nTRAIN = '../input/hpa-single-cell-image-classification/train'\npaths = [[os.path.join(TRAIN, train.iloc[idx,0])+ colour for colour in colours] for idx in range(len(train))]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Training data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Processing the data for training:\ntraining_data = []\nfor i,path in enumerate(paths[:500]):\n    img, label = image_prep(path, train['Label'][i])\n    training_data.append([img,label])\n\ntrain_ds = tf.data.Dataset.from_tensor_slices(([training_data[i][0] for i in range(len(training_data))], [training_data[i][1] for i in range(len(training_data))]))\nlen(train_ds)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Validation data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"val_data = []\nstart_img = 500\nval_num = 100\nfor i,path in enumerate(paths[start_img:start_img+val_num]):\n    img, label = image_prep(path, train['Label'][i+start_img])\n    val_data.append([img,label])\n\nval_ds = tf.data.Dataset.from_tensor_slices(([val_data[i][0] for i in range(len(val_data))], [val_data[i][1] for i in range(len(val_data))]))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **CNN Model**"},{"metadata":{"trusted":true},"cell_type":"code","source":"base_model = EfficientNetB0(include_top=False, weights='imagenet')\nbase_model.trainable = True\n\ninputs = layers.Input((IMG_SIZE[0], IMG_SIZE[0], 3))\n\nx = base_model(inputs, training=True)\nx = layers.GlobalAveragePooling2D()(x)\nx = layers.Dropout(0.5)(x)\noutputs = layers.Dense(len(LABELS), activation='sigmoid')(x)\n\ntf.keras.backend.clear_session()\n\nmodel = Model(inputs, outputs)\nmodel.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tf.keras.backend.clear_session()\n\nearlystopper = tf.keras.callbacks.EarlyStopping(\n    monitor='val_loss', patience=10, verbose=0, mode='min',\n    restore_best_weights=True\n)\n\nmodel.compile('adam', 'binary_crossentropy', metrics=[tf.keras.metrics.AUC(multi_label=True)])\n#model.compile('adam', 'binary_crossentropy', metrics=['accuracy'])\n#run = wandb.init(entity='ayush-thakur', project='hpa', job_type='train')\n\nhist = model.fit(train_ds, \n          epochs=50,\n          validation_data=val_ds,\n          verbose=1,\n          callbacks=[earlystopper]\n                )\n#plot_hist(hist)\n#run.finish()","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}