{"cells":[{"metadata":{},"cell_type":"markdown","source":"# 🗺️Getting our Bearings🗺️\n## Goal: Look at images of the same type and predict the labels of each individual cell within those images.\n* Task 1: Develop segmentation to cut images into images containing single-cells.\n* Task 2: Develop and train a CNN to predict labels (using green channel only) of single-cells\n* Task 3: Predict labels for all cells in \"test\" images - This is our submission."},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls ../input/hpa-single-cell-image-classification","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Imports\n\n#Much of the code used here is based on: https://keras.io/examples/keras_recipes/tfrecord/\n\nimport pandas as pd\nimport numpy as np\nimport tensorflow as tf\ntf.compat.v1.enable_eager_execution()\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport matplotlib.image as mpimg\nfrom mpl_toolkits.axes_grid1 import ImageGrid\n\nimport glob\n\nfrom PIL import Image\n\nfrom keras.models import Sequential\nfrom keras.layers import Conv2D, Conv3D, MaxPooling2D, Flatten, Dense, Activation, Dropout\nfrom keras.optimizers import Adam,SGD,Adagrad,Adadelta,RMSprop\n\nfrom sklearn.model_selection import train_test_split\nfrom functools import partial\n\nlabel_map = {'0': 'Nucleoplasm',\n             '1': 'Nuclear membrane',\n             '2': 'Nucleoli',\n             '3': 'Nucleoli fibrillar center',\n             '4': 'Nuclear speckles',\n             '5': 'Nuclear bodies',\n             '6': 'Endoplasmic reticulum',\n             '7': 'Golgi apparatus',\n             '8': 'Intermediate filaments',\n             '9': 'Actin filaments',\n             '10': 'Microtubules',\n             '11': 'Mitotic spindle',\n             '12': 'Centrosome',\n             '13': 'Plasma membrane',\n             '14': 'Mitochondria',\n             '15': 'Aggresome',\n             '16': 'Cytosol',\n             '17': 'Vesicles and punctate cytosolic patterns',\n             '18': 'Negative'\n            }\n\nsample_submit = pd.read_csv('../input/hpa-single-cell-image-classification/sample_submission.csv')\nsample_submit.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_csv = pd.read_csv('../input/hpa-single-cell-image-classification/train.csv')\ntrain_csv.head(3)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Spicy...we not only have balance issues, we have multiple cell types per image\ntrain_csv['Label'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Filtering down to a list of images with single cell types:\ntrain_csv[train_csv['Label'].isin(label_map.keys())]['Label'].value_counts()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# 🔬 Image Exploration 🔬\nThe first real challenge is going to be unsupervised segmentation so let's see what we're up against."},{"metadata":{"trusted":true},"cell_type":"code","source":"train_dir = '/kaggle/input/hpa-single-cell-image-classification/train_tfrecords/*'\ntest_dir = '/kaggle/input/hpa-single-cell-image-classification/test_tfrecords/*'\ntrain_list = glob.glob(train_dir)\ntest_list = glob.glob(test_dir)\n\nprint(\"Train Files: \", len(train_list))\nprint(\"Test Files: \", len(test_list))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# A utility script for reading features out of TFRecords:\n# I think I found it on stack exchange a while ago...\n\ndef list_record_features(tfrecords_path):\n    # Dict of extracted feature information\n    features = {}\n    # Iterate records\n    for rec in tf.data.TFRecordDataset([str(tfrecords_path)]):\n        # Get record bytes\n        example_bytes = rec.numpy()\n        # Parse example protobuf message\n        example = tf.train.Example()\n        example.ParseFromString(example_bytes)\n        # Iterate example features\n        for key, value in example.features.feature.items():\n            # Kind of data in the feature\n            kind = value.WhichOneof('kind')\n            # Size of data in the feature\n            size = len(getattr(value, kind).value)\n            # Check if feature was seen before\n            if key in features:\n                # Check if values match, use None otherwise\n                kind2, size2 = features[key]\n                if kind != kind2:\n                    kind = None\n                if size != size2:\n                    size = None\n            # Save feature data\n            features[key] = (kind, size)\n    return features\n\nrecord_features = list_record_features(glob.glob(train_dir+\"*\")[0])\nprint(record_features)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# We start with my basic data generator for unpacking tfrecords.\n# This will get tuned in for the dataset in question.\n\ndef decode_image(image):\n    image = tf.image.decode_jpeg(image)\n    image = tf.cast(image, tf.float32)\n    #image = tf.image.resize_with_crop_or_pad(image, 512, 512)\n    image = tf.image.resize(image, size=[im_size[0], im_size[1]])\n    return image\n\n\ndef read_tfrecord(example, labeled=True):\n    tfrecord_format = (\n        {\n            \"image\": tf.io.FixedLenFeature([], tf.string),\n            \"target\": tf.io.FixedLenFeature([], tf.string),\n            \"image_name\": tf.io.FixedLenFeature([], tf.string)\n        }\n        if labeled\n        else {\"image\": tf.io.FixedLenFeature([], tf.string),}\n    )\n    example = tf.io.parse_single_example(example, tfrecord_format)\n    image = decode_image(example[\"image\"])\n    if labeled:\n        #label = tf.one_hot(example[\"target\"], depth = 19)\n        label = example['target']\n        #return example\n        name = example[\"image_name\"]\n        return image, label, name\n    else:\n        #return example\n        return image\n\n\ndef load_dataset(filenames, labeled=True):\n    ignore_order = tf.data.Options()\n    ignore_order.experimental_deterministic = False  # disable order, increase speed\n    dataset = tf.data.TFRecordDataset(\n        filenames\n    )  # automatically interleaves reads from multiple files\n    dataset = dataset.with_options(\n        ignore_order\n    )  # uses data as soon as it streams in, rather than in its original order\n    dataset = dataset.map(\n        partial(read_tfrecord, labeled=labeled), num_parallel_calls=tf.data.experimental.AUTOTUNE\n    )\n    # returns a dataset of (image, label) pairs if labeled=True or just images if labeled=False\n    return dataset\n\n\ndef get_dataset(filenames, batch_size=10, labeled=True):\n    image_size = im_size\n    dataset = load_dataset(filenames, labeled)\n    # Normally I would shuffle the dataset here but we need our channels not to get jumbled (for now :p )\n    #dataset = dataset.shuffle(2048)\n    dataset = dataset.prefetch(buffer_size=tf.data.experimental.AUTOTUNE)\n    dataset = dataset.batch(batch_size)\n    return dataset","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"im_size = [200,200,4]\nbatch_size = 25\n\ntrain_dataset = get_dataset(train_list, batch_size=batch_size, labeled=True)\nvalid_dataset = get_dataset(test_list, batch_size=batch_size, labeled=True)\nimage_batch, label_batch, name_batch = next(iter(train_dataset))\n\ndef show_batch(image_batch, label_batch):\n    plt.figure(figsize=(10, 10))\n    for n in range(25):\n        ax = plt.subplot(5, 5, n + 1)\n        ax.title.set_text(str(label_batch[n].numpy()).replace(\"b\", \"\").replace(\"'\", \"\"))\n        plt.imshow(image_batch[n] / 255.0)\n        plt.axis(\"off\")\n        \nshow_batch(image_batch, name_batch)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch_size = 1000\ntrain_dataset = get_dataset(train_list, batch_size=batch_size, labeled=True)\nimage_batch, label_batch, name_batch = next(iter(train_dataset))\nname_batch\n\n# To produce composite images we run into something of a problem.\n# Our data-generator produces info out of order.\n# We would need to take the batch, match the images with their mates and provide the model with that stack for segmentation.\n# That's where I'll leave off 2/23","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Masking\nSo we found we have an R, G, B and Y channel reflecting tagged proteins in each image. Let's see if we can get the HPA-Cell-Segmentator working."},{"metadata":{"trusted":true,"collapsed":true},"cell_type":"code","source":"!pip install 'git+https://github.com/haoxusci/pytorch_zoo@master#egg=pytorch_zoo'\n!pip install 'https://github.com/CellProfiling/HPA-Cell-Segmentation/archive/master.zip'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image_batch[0].numpy()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Example useage from the HPA github:\n\nNUC_MODEL = \"./nuclei-model.pth\"\nCELL_MODEL = \"./cell-model.pth\"\n\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\n# Passing the numpy arrays of the images as they are (float32) causes: OverflowError: cannot convert float infinity to integer\n# TODO: Fix above error and replace this comment with the solution.\n\n# Testing:\n# From SkyDevour - takes composite image and returns image masks:\ndef get_masks(imgs, test=True):\n    try:\n        images = [[img[:, :, 0] for img in imgs], \n                  [img[:, :, 3] for img in imgs], \n                  [img[:, :, 2] for img in imgs],]\n    \n        nuc_segmentations = segmentator.pred_nuclei(images[2])\n        cell_segmentations = segmentator.pred_cells(images)\n        cell_masks = []\n        for i in tqdm(range(len(cell_segmentations)), desc='Labeling cells..'):\n            _, cell_mask = label_cell(nuc_segmentations[i], cell_segmentations[i])\n            cell_masks.append(cell_mask)\n        return cell_masks\n    except:\n        raise ValueError('Segmentation failed')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"composite = np.concatenate((image_batch[1].numpy(),image_batch[0].numpy(), image_batch[3].numpy()), axis = 2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(composite)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def vis_integrated_gradients_masks_test(img_idx, conf_threshold=0.01, mask_height=2048, mask_width=2048, \n                                        max_cell_level_conf_2_image_level_conf=0.005, test_ids=test_ids,\n                                        model=model_rgby, quantile_level=0.9, figsize=7):\n    image_id = test_ids[img_idx]\n    img = [cv2.resize(cv2.imread(os.path.join(TEST_IMGS_FOLDER, f'{image_id}_{color}.png'), cv2.IMREAD_GRAYSCALE),\n                      (mask_height, mask_width))\n           for color in ['red','green','blue','yellow']]\n    img = np.stack(img, axis=-1)\n    mask = get_masks([img])[0]\n    n_cells = mask.max()\n    cell_2_max_conf = dict()   \n    \n    img = cv2.resize(img, (IMG_HEIGHT, IMG_WIDTH)).astype(np.float32)/255.\n    predictions_test = model.predict(np.expand_dims(img, 0))\n    \n    for class_i, class_name in enumerate(class_names[:-1]):\n        class_conf_score = predictions_test[0][class_i]\n        if class_conf_score > conf_threshold:\n            try:\n                explanation = explainer.explain(([img], None), model, class_i, n_steps=15)\n                explanation_img = cv2.resize(explanation, (mask_height, mask_width))\n                explanation_total_level = np.quantile(explanation_img.flatten(), quantile_level)\n            except:\n                continue\n\n            plt.figure(figsize=(figsize, figsize))\n            plt.imshow(mask)\n            plt.imshow(explanation_img, alpha=0.7)\n            plt.xticks([])\n            plt.yticks([])\n            plt.title(f'{test_ids[img_idx]}\\n{class_name} ({class_conf_score:.2f}): raw Grad-CAMs', fontsize=22)\n            plt.show()\n\n            masks_all = np.zeros((mask_height, mask_width))\n            coord_2_conf = dict()\n            for cell_i in range(1, n_cells + 1):\n                cell_mask_bool = mask == cell_i\n                cell_explanation_perc = np.quantile(explanation_img[cell_mask_bool], quantile_level)\n                cell_conf = np.clip(cell_explanation_perc*class_conf_score/explanation_total_level, 0, class_conf_score)\n                if cell_conf/class_conf_score >= max_cell_level_conf_2_image_level_conf and cell_conf > 1e-3:\n                    masks_all[cell_mask_bool] = 1\n                    mask_pixels_x, mask_pixels_y = np.where(cell_mask_bool)\n                    coord_2_conf[(int(mask_pixels_y.mean()), int(mask_pixels_x.mean()))] = cell_conf\n                    if not cell_i in cell_2_max_conf:\n                        cell_2_max_conf[cell_i] = cell_conf\n                    else:\n                        cell_2_max_conf[cell_i] = max(cell_conf, cell_2_max_conf[cell_i])\n\n            plt.figure(figsize=(figsize, figsize)) \n            plt.imshow(masks_all)\n            for coords, conf in coord_2_conf.items():\n                conf_rounded = np.round(conf*100)/100\n                plt.scatter(*coords, s=700, color='red', marker=r\"$ {} $\".format(conf_rounded))\n            plt.xticks([])\n            plt.yticks([])\n            plt.title(f'{test_ids[img_idx]}\\n{class_name}: cell-level predictions', fontsize=22)\n            plt.show()\n\n    masks_all = np.zeros((mask_height, mask_width))\n    coord_2_conf = dict()\n    for cell_i in range(1, n_cells + 1):\n        if not cell_i in cell_2_max_conf:\n            cell_conf = 0.99\n        else:\n            cell_conf = 1 - cell_2_max_conf[cell_i]\n        if cell_conf >= conf_threshold:\n            cell_mask_bool = mask == cell_i\n            masks_all[cell_mask_bool] = 1\n            mask_pixels_x, mask_pixels_y = np.where(cell_mask_bool)\n            coord_2_conf[(int(mask_pixels_y.mean()), int(mask_pixels_x.mean()))] = cell_conf\n\n    plt.figure(figsize=(9, 9)) \n    plt.imshow(masks_all)\n    for coords, conf in coord_2_conf.items():\n        conf_rounded = np.round(conf*100)/100\n        plt.scatter(*coords, s=700, color='red', marker=r\"$ {} $\".format(conf_rounded))\n    plt.xticks([])\n    plt.yticks([])\n    plt.title(f'{test_ids[img_idx]}\\n{class_names[-1]}: cell-level predictions', fontsize=22)\n    plt.show()\n    \n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"type(np.asarray(images[CHANNEL_BLUE])[0][0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(segmentator.pred_nuclei([(np.asarray(images[CHANNEL_BLUE]))])[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"np.asarray(images[CHANNEL_BLUE])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"image = image_batch[0].numpy()\nimage = tf.cast(image, tf.float32)\nimage = image.numpy()\nimage\n#segmentator.pred_nuclei(image)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"nuclei_mask = label_nuclei(nuc_segmentations[0])\nprint(np.shape(nuclei_mask))\n\ncell_nuclei_mask, cell_mask = label_cell(nuc_segmentations[0], cell_segmentations[0])\nprint(np.shape(nuclei_mask))\nprint(np.shape(cell_mask))\n\nfig = plt.figure(figsize=(25,25))\n\n# Nuclei mask\nnuclei_image = Image.fromarray( np.uint8(nuclei_mask) )\n\nax = fig.add_subplot(1, 4, 1)\nax.set_title(\"Nuclei Mask\")\nplt.imshow(np.asarray(nuclei_image))\n\n# Cell nuclei mask\ncell_nuclei_image = Image.fromarray( np.uint8(cell_nuclei_mask) )\n\nax = fig.add_subplot(1, 4, 2)\nax.set_title(\"Cell Nuclei Mask\")\nplt.imshow(np.asarray(cell_nuclei_image))\n\n# Cell mask\ncell_image = Image.fromarray( np.uint8(cell_mask) )\n\nax = fig.add_subplot(1, 4, 3)\nax.set_title(\"Cell Mask\")\nplt.imshow(np.asarray(cell_image))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Get unique vector of segment numbers\nnumbers = set( np.ravel(cell_mask) )\nnumbers.remove(0)\n\nfig = plt.figure(figsize=(25,6*len(numbers)/4))\nindex = 1\n\n# plot original cell mask from above\nax = fig.add_subplot((len(numbers)//4)+1, 4, index)\nax.set_title(\"Complete Cell Mask\")\nplt.imshow(np.asarray(cell_image))\nindex = index + 1\n\nfor number in numbers:\n    # set all other 'numbers' to zero in cell mask\n    isolated = np.where(cell_mask == number, cell_mask, 0)\n\n    # plot isolated image\n    ax = fig.add_subplot((len(numbers)//4)+1, 4, index)\n    ax.set_title(\"Segment: {}\".format(number))\n\n    plt.imshow(isolated)\n    index = index + 1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Classification Model\nNow that we have a segmentation strategy we can bake it into our data-generator."},{"metadata":{"trusted":true},"cell_type":"code","source":"# So let's tweak our image decoder slightly.\n# We want to segment cells, then only pass the green channel back since that's what we're using to predict classes.\n\ndef decode_image(image):\n    image = tf.image.decode_jpeg(image)\n    image = tf.cast(image, tf.float32)\n    image = tf.image.resize(image, size=[im_size[0], im_size[1]])\n    return image","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#So now our generator is giving us:\n\nim_size = [200,200,4]\nbatch_size = 25\n\ntrain_dataset = get_dataset(train_list, batch_size=batch_size, labeled=True)\nvalid_dataset = get_dataset(test_list, batch_size=batch_size, labeled=True)\nimage_batch, label_batch = next(iter(train_dataset))\n\ndef show_batch(image_batch, label_batch):\n    plt.figure(figsize=(10, 10))\n    for n in range(25):\n        ax = plt.subplot(5, 5, n + 1)\n        ax.set_title(label_batch[n])\n        plt.imshow(image_batch[n] / 255.0)\n        plt.axis(\"off\")\n        \nshow_batch(image_batch, label_batch)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#Global variables for training\nBATCH_SIZE = 25\nim_size = [299, 299, 3]\nkeras_batch_size = BATCH_SIZE\nlearning_rate = 0.095\nepochs = 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"sample_train = next(iter(train_dataset))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"str_format = str(sample_train[1][4].numpy())\nstr_format.r","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So far we have a dataset of images paired to lists of correct outputs contained in strings with format b'x|y|z...'"},{"metadata":{"trusted":true},"cell_type":"code","source":"str(sample_train[1][0].numpy())","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = tf.keras.applications.InceptionV3(\n    include_top=True,\n    weights=None,\n    input_tensor=None,\n    input_shape=None,\n    pooling=None,\n    classes=1,\n    classifier_activation=\"softmax\",\n)\n\nmodel.compile(optimizer='adam', loss=tf.keras.losses.categorical_crossentropy,\n              metrics=tf.keras.metrics.categorical_accuracy)\n\nmodel.fit(train_dataset, validation_data=valid_dataset, batch_size=keras_batch_size, epochs=epochs)","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}