{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction"},{"metadata":{},"cell_type":"markdown","source":"In remote sensing, it is common to design and use sensors which are sensitive to different spectral wavelengths than the visible spectrum. These types of sensors can be qualified as multispectral or hyperspectral, depending on the spectral resolution. Multispectral cameras usually have the usual RGB channels, and a few other spectral bands. The additional spectral bands will be selected in function of the phenonema that we want to observe. For example, in the case of plants, it has been demonstrated that plant health can be evaluated using the [Normalized Difference Vegetation Index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index). Due to these sensors being highly specific, methods have been developed to highlight green vegetation in imagery from widely-used RGB cameras. By using such vegetation indices for RGB images and an Otsu thresholding method, the vegetation can be extracted from RGB images purely by solely using \"traditional\" image processing methods."},{"metadata":{},"cell_type":"markdown","source":"The three main objective of this notebook are:\n* Use two of the best performing vegetation indices for RGB images mentioned in [this recent paper](https://link.springer.com/article/10.1007/s11119-020-09735-1).\n* Try to train a fully-convolutional Unet by filtering the vegetation in images using the best vegetation index. Once the groundtruth is extracted, we can apply advanced data augmentation techniques on the images of healthy plants to mimic diseases. We proceed this way as diseases will lead to diverse leaf coloration types that will negatively impact our vegetation indices. Therefore, we might obtain a more robust model that also segmentates diseased leaves correctly, where vegetation indices cannot.\n* Open the conversation around the use of vegetation indices for disease classification, and whether they can actually help. As some other deep learning models have been improved thanks to preprocessing (e.g. [automatic color equilization in medical imagery problems](https://ieeexplore.ieee.org/abstract/document/9162010)), I have good hopes something can be done here too!"},{"metadata":{"trusted":true},"cell_type":"code","source":"!git clone https://github.com/jakeret/unet\n!pip install ../working/unet/","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom PIL import Image\nimport os\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm\nfrom sklearn.utils import shuffle\nfrom sklearn.utils import class_weight\nfrom sklearn.preprocessing import minmax_scale\nimport random\nimport cv2\nfrom imgaug import augmenters as iaa\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport tensorflow as tf\nfrom tensorflow.keras import layers\nfrom tensorflow.keras import Input\nfrom tensorflow.keras.models import Model\nfrom tensorflow.keras.layers import Dense, Flatten, Dropout, Activation\nfrom tensorflow.keras.layers import BatchNormalization, GlobalAveragePooling2D\nfrom tensorflow.keras.callbacks import ModelCheckpoint, ReduceLROnPlateau, EarlyStopping\n\nfrom unet import utils\nfrom unet.datasets import circles\nimport unet","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_folder = '../input/cassava-leaf-disease-classification/train_images/'\nsamples_df = pd.read_csv(\"../input/cassava-leaf-disease-classification/train.csv\")\nsamples_df[\"label\"] = samples_df[\"label\"].astype(\"str\")\nsamples_df.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Only keep images of healthy plants as the vegetation indices may not work as well on diseased plants with yellow and brown patches."},{"metadata":{"trusted":true},"cell_type":"code","source":"samples_df = samples_df.query(\"label=='4'\")","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_percentage = 0.8\ntraining_item_count = int(len(samples_df)*training_percentage)\nvalidation_item_count = len(samples_df)-int(len(samples_df)*training_percentage)\ntraining_df = samples_df[:training_item_count]\nvalidation_df = samples_df[training_item_count:]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test 2 vegetation indices for RGB images"},{"metadata":{},"cell_type":"markdown","source":"Below, we are implementing functions that return the process vegetation indices:\n* Elliptical Color Index (ECI) - [link to the paper](https://link.springer.com/article/10.1007/s11119-020-09735-1)\n* Color Index Vegetation Extraction (CIVE) - [link to the paper](https://ieeexplore.ieee.org/abstract/document/1225492)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_ECI_band(img):\n    '''\n    Return the ECI band calculated from an RGB image between 0 and 255\n    using the formula below:    \n    ECI = (red_channel-1)^2 + green_channel^2/0.16\n    '''\n    img = img/255.\n    img = cv2.GaussianBlur(img,(35,35),0)\n    ECI_band = np.power(img[:,:,0]-1,2) + np.power(img[:,:,1],2)/0.16\n    normalized_ECI_band = (ECI_band/ECI_band.max()*255).astype(np.uint8)\n    return normalized_ECI_band\n\n\ndef get_CIVE_band(img):\n    '''\n    Return the CIVE band calculated from an RGB image between 0 and 255\n    using the formula below:\n    CIVE = 0.441*red_channel - 0.881*green_channel + 0.385*blue_channel + 18.787\n    '''\n    img = cv2.GaussianBlur(img,(35,35),0)\n    CIVE_band = 0.441*img[:,:,0] - 0.881*img[:,:,1] + 0.385*img[:,:,2] + 18.787\n    normalized_CIVE_band = (((CIVE_band+abs(CIVE_band.min()))/CIVE_band.max())).astype(np.uint8)\n    return normalized_CIVE_band\n\n\ndef apply_ECI_mask(img, vegetation_index_band):\n    '''\n    Apply a binary mask on an image and return the masked image\n    '''\n    ret, otsu = cv2.threshold(vegetation_index_band,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n    masked_img = cv2.bitwise_and(img,img,mask = otsu)\n    return masked_img\n\n\ndef apply_CIVE_mask(img, vegetation_index_band):\n    '''\n    Apply a binary mask on an image and return the masked image\n    '''\n    ret, otsu = cv2.threshold(vegetation_index_band,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)\n    masked_img = cv2.bitwise_and(img,img,mask = otsu)\n    return masked_img","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's try to mask the non-vegetation using ECI+Otsu."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 20))\nitems = 6\nfor idx, image_id in enumerate(samples_df.image_id[:items]):\n    img_path = training_folder+image_id\n    img = np.array(Image.open(img_path))\n    ax = plt.subplot(items, 3, idx*3 + 1)\n    ax.set_title(\"original\")\n    plt.imshow(img)\n    \n    ECI_band =  get_ECI_band(img)\n    ax = plt.subplot(items, 3, idx*3 + 2)\n    ax.set_title(\"ECI band\")\n    plt.imshow(ECI_band)\n    \n    masked_img = apply_ECI_mask(img, ECI_band)\n    ax = plt.subplot(items, 3, idx*3 + 3)\n    ax.set_title(\"ECI+Otsu\")\n    plt.imshow(masked_img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now, let's assess the performance of CIVE+Otsu on Cassava leaves segmentation."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 20))\nitems = 6\nfor idx, image_id in enumerate(samples_df.image_id[:items]):\n    img_path = training_folder+image_id\n    img = np.array(Image.open(img_path))\n    ax = plt.subplot(items, 3, idx*3 + 1)\n    ax.set_title(\"original\")\n    plt.imshow(img)\n    \n    CIVE_band = get_CIVE_band(img)\n    ax = plt.subplot(items, 3, idx*3 + 2)\n    ax.set_title(\"CIVE band\")\n    plt.imshow(CIVE_band)\n    \n    masked_img = apply_CIVE_mask(img, CIVE_band)\n    ax = plt.subplot(items, 3, idx*3 + 3)\n    ax.set_title(\"CIVE+Otsu\")\n    plt.imshow(masked_img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At first, it appears that CIVE is more precise than ECI. When put them side by side for comparison as below, the CIVE+Otsu method does seem to be outperforming the ECI+Otsu method."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 20))\nitems = 6\nfor idx, image_id in enumerate(samples_df.image_id[10:10+items]):\n    img_path = training_folder+image_id\n    img = np.array(Image.open(img_path))\n    ax = plt.subplot(items, 3, idx*3 + 1)\n    ax.set_title(\"original\")\n    plt.imshow(img)\n    \n    ECI_band =  get_ECI_band(img)\n    CIVE_band = get_CIVE_band(img)\n\n    masked_img = apply_ECI_mask(img, ECI_band)\n    ax = plt.subplot(items, 3, idx*3 + 2)\n    ax.set_title(\"ECI+Otsu\")\n    plt.imshow(masked_img)\n    \n    masked_img = apply_CIVE_mask(img, CIVE_band)\n    ax = plt.subplot(items, 3, idx*3 + 3)\n    ax.set_title(\"CIVE+Otsu\")\n    plt.imshow(masked_img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Build a Cassava leaf segmentation model"},{"metadata":{},"cell_type":"markdown","source":"From here, I will build a custom generator to feed data into a [Unet](https://github.com/damienpontifex/tf-keras-unet). The custom generator will randomly crop a 512x512 images from the original images, before applying CIVE and the Otsu thresholding to only extract the vegetation mask. The vegetation mask will be kept as groundtruth. As CIVE only works correctly on healthy plants, I have filtered out all diseased samples but will, instead, apply data augmentation techniques mimicking diseases (or at least, trying!). "},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_input_image(img_path, crop_size=512):\n    '''\n    Randomly select a 512x512 pixel area from the full image\n    '''\n    img = Image.open(img_path)\n    img_height, img_width = img.size\n    img = np.array(img)\n    y = random.randint(0,img_height-crop_size)\n    x = random.randint(0,img_width-crop_size)\n    cropped_img = img[x:x+crop_size , y:y+crop_size,:]\n    \n    return cropped_img\n\n\ndef get_groundtruth_mask(img, crop_size=512):\n    '''\n    Generate the groundtruth mask using the CIVE band\n    '''\n    img = cv2.GaussianBlur(img,(35,35),0)\n    cive_band = 0.441*img[:,:,0] - 0.881*img[:,:,1] + 0.385*img[:,:,2] + 18.787\n    normalized_cive_band = (((cive_band+abs(cive_band.min()))/cive_band.max())).astype(np.uint8)\n    ret, otsu_mask = cv2.threshold(normalized_cive_band,0,1,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)\n    \n    veg_mask = otsu_mask.astype(np.uint8)\n    soil_mask = np.where(otsu_mask==1, 0, 1).astype(np.uint8)\n    masks = np.transpose(np.array([soil_mask, veg_mask]),(1,2,0))\n    \n    return masks","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"augmentation = iaa.Sequential([\n    iaa.Multiply((0.5, 1.5)),\n    iaa.Affine(scale={\"x\": (1, 1.2), \"y\":(1, 1.2)}),\n    iaa.Sometimes(0.4,\n                 iaa.GaussianBlur(sigma=(0,2))),\n    iaa.Sometimes(0.3,\n                 iaa.Grayscale(alpha=(0.0, 1.0))),\n    iaa.Sometimes(0.3,\n                 iaa.SigmoidContrast(gain=(3, 6), cutoff=(0.4, 0.6))),\n    iaa.Sometimes(0.3,\n        iaa.CoarseDropout((0.0, 0.05), size_percent=(0.05, 0.6), per_channel=0.5))\n])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"images = [get_input_image(training_folder+image_filename) for image_filename in samples_df[:10].image_id.values]\nlabels = samples_df[:10].label.values\naugmented_images = augmentation.augment_images(images=images)\n\nsample_number = len(augmented_images)\nfig = plt.figure(figsize = (20,sample_number))\nfor i in range(0,sample_number):\n    ax = fig.add_subplot(2, 5, i+1)\n    ax.imshow(augmented_images[i])\n    ax.set_title(str(labels[i]))\nplt.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def custom_generator(image_path_list, folder, batch_size=16, training_mode=True):\n    \n    while True:\n        for start in range(0, len(image_path_list), batch_size):\n            X_batch = []\n            Y_batch = []\n            end = min(start + batch_size, training_item_count)\n            \n            image_list = [get_input_image(folder+\"/\"+image_path) for image_path in image_path_list[start:end]]\n            groundtruth_image_list = [get_groundtruth_mask(image) for image in image_list]\n            \n            #only apply augmentation during training\n            if training_mode:\n                image_list = augmentation.augment_images(images=image_list)\n\n            X_batch = np.array(image_list)/255.\n            Y_batch = np.array(groundtruth_image_list)\n\n            yield X_batch, Y_batch","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"batch_size = 8\ncallbacks = [ReduceLROnPlateau(monitor='val_loss', patience=1, verbose=1, factor=0.5),\n             EarlyStopping(monitor='val_loss', patience=3),\n             ModelCheckpoint(filepath='best_model.h5', monitor='val_loss', save_best_only=True)]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"unet_model = unet.build_model(512,\n                              channels=3,\n                              num_classes=2,\n                              layer_depth=4,\n                              filters_root=64,\n                              padding=\"same\")\n\nunet.finalize_model(unet_model,\n                    loss=tf.keras.losses.BinaryCrossentropy(),\n                    metrics=[tf.keras.metrics.BinaryAccuracy()],\n                    auc=False,\n                    learning_rate=1e-4)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = unet_model.fit_generator(custom_generator(training_df[\"image_id\"], training_folder, batch_size=batch_size, training_mode=True),\n                  steps_per_epoch = int(len(training_df)/batch_size),\n                  epochs = 20, \n                  validation_data=custom_generator(validation_df[\"image_id\"], training_folder, batch_size=batch_size),\n                  validation_steps=int(len(validation_df)/batch_size),\n                  callbacks=callbacks)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"unet_model.load_weights(\"best_model.h5\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualising results from our segmentation model"},{"metadata":{},"cell_type":"markdown","source":"Once the model is trained, I can compare the results against the computer vision technique. It will be especially interesting to see what happens on diseased leaves."},{"metadata":{"trusted":true},"cell_type":"code","source":"!cat ../input/cassava-leaf-disease-classification/label_num_to_disease_map.json","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"training_folder = '../input/cassava-leaf-disease-classification/train_images/'\ndisease_df = pd.read_csv(\"../input/cassava-leaf-disease-classification/train.csv\")\ndisease_df = disease_df.query(\"label!=4\")","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"First, I proceed to a visual assessment on healthy Cassava leaves."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 20))\nitems = 5\nfor idx, image_id in enumerate(validation_df.image_id[:items]):\n    img_path = training_folder+image_id\n    img = get_input_image(img_path)\n    ax = plt.subplot(items, 3, idx*3 + 1)\n    ax.set_title(\"original\")\n    plt.imshow(img)\n    \n    CIVE_band = get_CIVE_band(img)\n    masked_img = apply_CIVE_mask(img, CIVE_band)\n    ax = plt.subplot(items, 3, idx*3 + 2)\n    ax.set_title(\"CIVE+Otsu\")\n    plt.imshow(masked_img)\n\n    mask = unet_model.predict(np.array([img]))\n    soil_mask = np.where(np.transpose(mask[0],(2,0,1))[0]<=0.5,1,0).astype(np.uint8)\n    masked_img = cv2.bitwise_and(img,img,mask = soil_mask)\n    ax = plt.subplot(items, 3, idx*3 + 3)\n    ax.set_title(\"Unet\")\n    plt.imshow(masked_img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"And then on diseased leaves."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(20, 20))\nitems = 5\nfor idx, image_id in enumerate(disease_df.image_id[10:10+items]):\n    img_path = training_folder+image_id\n    img = get_input_image(img_path)\n    ax = plt.subplot(items, 3, idx*3 + 1)\n    ax.set_title(\"original\")\n    plt.imshow(img)\n    \n    CIVE_band = get_CIVE_band(img)\n    masked_img = apply_CIVE_mask(img, CIVE_band)\n    ax = plt.subplot(items, 3, idx*3 + 2)\n    ax.set_title(\"CIVE+Otsu\")\n    plt.imshow(masked_img)\n\n    mask = unet_model.predict(np.array([img]))\n    soil_mask = np.where(np.transpose(mask[0],(2,0,1))[0]<=0.5,1,0).astype(np.uint8)\n    masked_img = cv2.bitwise_and(img,img,mask = soil_mask)\n    ax = plt.subplot(items, 3, idx*3 + 3)\n    ax.set_title(\"Unet\")\n    plt.imshow(masked_img)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It is complicated to make any finite conclusion purely from visual assessment and we do not have any groundtruth. Nonetheless, the Unet does seem to perform a bit better than the original computer vision methods, especially on leaves in the background and in messy pictures. I could imagine this model being used as a preprocessing method for the Cassava leaf disease classification problem, and I would be very keen to hear from anyone trying this out."},{"metadata":{},"cell_type":"markdown","source":"### Thanks for reading this notebook! If you found this notebook helpful, please give it an upvote. It is always greatly appreciated!"}],"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}