{"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\nimport seaborn as sns\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.\ndef get_blended_image(images): \n    \n    # blend rgby images into single array\n    blended_array = np.stack(images, 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":"# Input: image arrays\n# Output: normalized image arrays\ndef normalization(images):\n    \n    norm_images = []\n    for array in images:\n        norm_images.append(array/255)\n        \n    return images","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Introduce list of image filters as .png\n# Returns a processed image ready for the CNN and an encoded label as tensor\ndef image_prep(paths, label):\n    \n    # Transformation from .png to arrays\n    img = image_to_arrays(paths)\n    size = np.shape(img[0])[0]\n    # Image normalization\n    img = normalization(img)\n    # Process the arrays to obtain tensor of desired shapes and size\n    img = tf.image.convert_image_dtype(img, dtype=tf.float32)\n    img = tf.reshape(img, (size, size, 3))\n    img = tf.image.resize(img, IMG_SIZE)\n    \n    # Process the labels to obtain a tensor of the shapes required by the CNN model\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    \n    return img, label","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def apply_augmentation(image, label):\n    \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":"# Function for plotting the history of the CNN training\ndef 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        label_counts.append(int(num))\ncounts = pd.value_counts(label_counts)\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.title('Label counts')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The following graphs come from the notebook [HPA Single Cell Classification EDA【中文】](https://www.kaggle.com/wptouxx/hpa-single-cell-classification-eda)"},{"metadata":{"trusted":true},"cell_type":"code","source":"train['num_classes'] = train['Label'].apply(lambda r: len(r.split('|')))\ntrain['num_classes'].value_counts().plot.bar(title='Examples with multiple labels', xlabel='number of labels per example', ylabel='# train examples')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"labels = [str(i) for i in range(19)]\nfor x in labels: train[x] = train['Label'].apply(lambda r: int(x in r.split('|')))\n    \nunique_counts = {}\nfor lbl in labels:\n    unique_counts[lbl] = len(train[train.Label == lbl])\n\nfull_counts = {}\nfor lbl in labels:\n    count = 0\n    for row_label in train['Label']:\n        if lbl in row_label.split('|'): count += 1\n    full_counts[lbl] = count\n    \ncounts = list(zip(full_counts.keys(), full_counts.values(), unique_counts.values()))\ncounts = np.array(sorted(counts, key=lambda x:-x[1]))\ncounts = pd.DataFrame(counts, columns=['label', 'full_count', 'unique_count'])\n\ndf_labels = train[labels]\ncoocc = df_labels.T.dot(df_labels)\nfor i in range(19): coocc.iloc[i,i] = int(counts.unique_count[counts.label == str(i)].values[0])\n\nfig, ax = plt.subplots(figsize=(16, 13))\nsns.heatmap(coocc, cmap=\"Blues\", linewidth=0.3, cbar_kws={\"shrink\": .8})\ntitle = 'How often do individual classes cooccur in the train set?\\n(Single class occurrences are on the diagonal)\\n'\nplt.title(title, loc='center', fontsize=18)\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        plt.axis('off')\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.axis('off')\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.axis('off')\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":"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 tensorflow as tf\nfrom tensorflow.keras.applications import EfficientNetB0\nfrom tensorflow.keras.models import Sequential\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import Model\nimport wandb\nfrom tqdm import tqdm","execution_count":null,"outputs":[]},{"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":{},"cell_type":"markdown","source":"First, I want to create a separate file with the stacked filters of every image. This will accelerate the process"},{"metadata":{"trusted":true},"cell_type":"code","source":"DATA_DIR = \"/kaggle/input/hpa-single-cell-image-classification\"\ntrain = pd.read_csv(os.path.join(DATA_DIR,'train.csv'))\n\ncolours = ['_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))]\n\nif not os.path.exists('../train_blend'):\n    os.makedirs('../train_blend')\n    \n    for image in paths:\n        im = get_blended_image(image_to_arrays(image))\n        ID = image[0].split('_')[0].split('/')[-1]\n        im.save(\"../train_blend/{ID}.png\")\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"run = wandb.init(project='hpa', job_type='consume_split')\nartifact = run.use_artifact('glopezzz/hpa/split:v0', type='dataset')\nSPLIT_CSV_PATH = artifact.download()\nrun.finish()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_train_split = pd.read_csv(SPLIT_CSV_PATH+'/train_split.csv')\ndf_val_split = pd.read_csv(SPLIT_CSV_PATH+'/val_split.csv')\n\nprint(df_train_split.head(), '\\n', df_val_split.head())\nprint('\\nTrain split: ',len(df_train_split))\nprint('Val split: ',len(df_val_split))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Parameter setting"},{"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', '_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":"# **Training data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Processing the data for training:\ntraining_length = 5000\ntraining_data = []\nfor i,path in tqdm(enumerate(paths[:training_length])):\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))]))\n\ntrain_ds = train_ds.batch(batch_size=BATCH_SIZE, drop_remainder=True)\ntrain_ds = train_ds.prefetch(tf.data.experimental.AUTOTUNE)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# **Validation data**"},{"metadata":{"trusted":true},"cell_type":"code","source":"val_num = 1000\nval_data = []\nfor i,path in tqdm(enumerate(paths[training_length:training_length+val_num])):\n    img, label = image_prep(path, train['Label'][i+training_length])\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))]))\n\nval_ds = val_ds.batch(batch_size=BATCH_SIZE, drop_remainder=True)\nval_ds = val_ds.prefetch(tf.data.experimental.AUTOTUNE)","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\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":{},"cell_type":"markdown","source":"Now it's time to save the model for later use."},{"metadata":{"trusted":true},"cell_type":"code","source":"model.save('effnet_multilabel.h5')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The continuation of the project is implemented in my notebook: [Test classification and submission](https://www.kaggle.com/glopezzz/test-classification-and-submission)\n\nThere, I'll use the trained model for cell labeling using the test data.\nI'll also prepare the output for submission, with the appropiate mask encoding.\n\nGive it a look. Don't hessitate in asking anything :)"}],"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}