{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9e1e5739",
   "metadata": {},
   "source": [
    "# Exercise 4: Image classification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3115cfe3",
   "metadata": {},
   "source": [
    "## 0. Activate your virtual environment\n",
    "1. Since the Noto update, it is not possible anymore to run pip commands without having first activated your personal virtual environment. Follow the guidelines provided in `🗀 / Documentation / 10_Envs_and_kernels.ipynb` to create your venv\n",
    "2. Once your environment has been created, you can select it for the current notebook on the top right of the window: replace `Python` by the `Py3 Pint` environment that your have created."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1936ed2",
   "metadata": {},
   "source": [
    "## 1. Install dependencies, if needed\n",
    "\n",
    "Open a Terminal by going on `File -> New -> Terminal`, then activate the environment that you have created with `my_venvs_activate pint_env`. Finally, run the command below in the terminal to install the dependencies:\n",
    "```bash\n",
    "pip install -q pillow matplotlib scikit-image scikit-learn gdown\n",
    "```\n",
    "At this point, you can come back to this Jupyter Notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49444b14",
   "metadata": {},
   "source": [
    "## 2. Extract regions and features\n",
    "\n",
    "- Download the images by running the next cell "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "73056b12",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gdown\n",
    "import os\n",
    "import zipfile\n",
    "\n",
    "if not os.path.exists(\"images_week4\"):\n",
    "    id = \"1n8wf8K7LK_F8q1SUBQB3HX7eWtumySfc\"\n",
    "    gdown.download(id=id, output=\"images_week4.zip\", quiet=True)\n",
    "    with zipfile.ZipFile(\"images_week4.zip\", 'r') as zip_ref:\n",
    "        zip_ref.extractall()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e212ac1",
   "metadata": {},
   "source": [
    "Alternatively, you can download it directly from [here](https://drive.google.com/file/d/1n8wf8K7LK_F8q1SUBQB3HX7eWtumySfc/view?usp=sharing) via the browser and unzip in this folder.\n",
    "\n",
    "- The zip file contains the following subfolders:\n",
    "    - images\n",
    "    - gt: Ground truth (GT) labels for each pixel of the image, labeled with the following classes\n",
    "        - 0. Impervious\n",
    "        - 1. Building\n",
    "        - 2. Low vegetation\n",
    "        - 3. Tree\n",
    "- Run the code block of functions that extract features (it is the same code you implemented in the past exercise)\n",
    "- Extract and save regions for each image using SLIC (use the parameter n_segments=1000)\n",
    "- Extract and save features for the regions of each image"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c17d7a84",
   "metadata": {},
   "source": [
    "warmup: let's load the images and display them again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5015c05",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # for plotting\n",
    "from skimage.io import imsave, imread # to read and save images\n",
    "import os # to list and search folders\n",
    "\n",
    "# TODO: list all files within the \"images\" subfolder and within the \"gt\" subfolder\n",
    "images_list = os.listdir(\"todo\")\n",
    "ground_truth_list = os.listdir(\"todo\")\n",
    "\n",
    "for image_filename, ground_truth_filename in zip(images_list, ground_truth_list):\n",
    "    fig, axs = plt.subplots(1,2, figsize=(6,4))\n",
    "    \n",
    "    # first subplot\n",
    "    ax = axs[0]\n",
    "    \n",
    "    # TODO: specify the path to this image file\n",
    "    image = imread(f\"todo\")\n",
    "    \n",
    "    ax.imshow(image)\n",
    "    \n",
    "    # second subplot\n",
    "    ax = axs[1]\n",
    "    \n",
    "    # TODO: specify the path to this ground truth file\n",
    "    image = imread(f\"todo\")\n",
    "    \n",
    "    ax.imshow(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d77aac7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The following functions are the same you implemented the past week\n",
    "# You do not need to edit the functions of this code block (just run it)\n",
    "def convert_to_shape_pixels_by_bands(data):\n",
    "    num_dimensions = len(data.shape)\n",
    "    assert(num_dimensions == 2 or num_dimensions == 3)\n",
    "    if num_dimensions == 3:\n",
    "        num_bands = data.shape[2]\n",
    "        return data.reshape((-1, num_bands))\n",
    "    else:\n",
    "        return data\n",
    "\n",
    "def compute_average_feature(data):\n",
    "    # If needed convert data to the shape (num_pixels x num_bands)\n",
    "    data_2d = convert_to_shape_pixels_by_bands(data)\n",
    "    num_bands = data_2d.shape[1]\n",
    "    avg_features = np.mean(data_2d, axis=0)\n",
    "    return avg_features\n",
    "\n",
    "def compute_standard_deviation_feature(data):\n",
    "    # If needed convert data to the shape (num_pixels x num_bands)\n",
    "    data_2d = convert_to_shape_pixels_by_bands(data)\n",
    "    avg_features = np.std(data_2d, axis=0)\n",
    "    return avg_features\n",
    "\n",
    "def compute_histogram_feature(data, num_bins=10):\n",
    "    # If needed convert data to the shape (num_pixels x num_bands)\n",
    "    data_2d = convert_to_shape_pixels_by_bands(data)\n",
    "    num_bands = data_2d.shape[1]\n",
    "    hist_features = np.zeros((num_bands, num_bins)).astype(np.float32)\n",
    "    for b in range(num_bands):\n",
    "        # Compute the histogram for each band \n",
    "        #       use the function np.histogram(array, bins=num_bins)\n",
    "        hist, boundaries = np.histogram(data_2d[:, b], bins=num_bins)\n",
    "        hist_features[b, :] = hist\n",
    "    # Return a 1D array containing all the values\n",
    "    return hist_features.flatten()\n",
    "\n",
    "def compute_image_features_from_regions(image, segmentation_map):\n",
    "    region_id_list = np.unique(segmentation_map)\n",
    "    all_features = []\n",
    "    for region_id in region_id_list:\n",
    "        # Obtain pixel values of each regions, with shape (num_pixels x num_bands)\n",
    "        pixel_values = image[segmentation_map==region_id]\n",
    "        # Compute the average, standard deviation and histogram features\n",
    "        #       and concatenated them unsing the function (np.concatenate)\n",
    "        avg = compute_average_feature(pixel_values)\n",
    "        features = compute_standard_deviation_feature(pixel_values)\n",
    "        hist_features = compute_histogram_feature(pixel_values)\n",
    "        features = np.concatenate([avg, features, hist_features])\n",
    "        # Add concatenated features to the variable all_features\n",
    "        all_features.append(features)\n",
    "    # convert list to numpy array of shape: (num_regions x num_bands)\n",
    "    return np.array(all_features).astype(np.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "090c8a98",
   "metadata": {},
   "source": [
    "Let's extract regions and their features (average, standard deviation and histogram of each band)\n",
    "\n",
    "1. use the [slic](https://scikit-image.org/docs/dev/api/skimage.segmentation.html?highlight=slic#skimage.segmentation.slic) algorithm from skimage\n",
    "2. use `compute_image_features_from_regions` to extract features\n",
    "3. save numpy arrays with `np.save(filename.npy)` and images with `imsave(filename.tif)`\n",
    "\n",
    "Tip: you can replace substrings in `text` with `text.replace(\"a\",\"b\")` this will be useful when renaming files (e.g. `images_1.tif` to `regions_1.tif`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e1e6357",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from skimage.segmentation import slic\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "for image_filename in images_list:\n",
    "    \n",
    "    # TODO: Read image\n",
    "    image_path = \n",
    "    \n",
    "    # TODO: Segment image using SLIC (use the parameter n_segments=1000)\n",
    "    segmented_image = \n",
    "    \n",
    "    # Save SLIC regions with imsave (e.g. in images_week4/regions/regions_*.tif)\n",
    "    # rename images_1.tif -> regions_1.tif\n",
    "    regions_filename = image_filename.replace(\"image\",\"regions\")\n",
    "    regions_path = f\"images_week4/regions/{regions_filename}\"\n",
    "    imsave(regions_path, segmented_image.astype(np.uint32))\n",
    "    print(\"regions saved in \" + regions_path)\n",
    "    \n",
    "    # TODO: Compute features with compute_image_features_from_regions\n",
    "    region_features = \n",
    "    \n",
    "    # Save features\n",
    "    # rename images_1.tif -> features_1.npy\n",
    "    features_filename = image_filename.replace(\"image\",\"features\").replace(\".tif\",\".npy\")\n",
    "    region_features_path = f\"images_week4/features/{features_filename}\"\n",
    "    np.save(region_features_path, region_features)\n",
    "    print(f\"features saved in {region_features_path}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "002675af",
   "metadata": {},
   "source": [
    "## 3. Create training dataset\n",
    "- Define which images will be in the training set, and which ones the test set\n",
    "- Read the images from the training set and compute the ideal label for each region (the label that has the largest intersection with each region)\n",
    "- Create an array of targets (groud truth), per region, that combines the targets of all the training images\n",
    "- Create an array of features that combines the features of all the training images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc361d6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: split images randomly into training and testing partitions\n",
    "\n",
    "# fixing seed to have same results\n",
    "np.random.seed(0)\n",
    "\n",
    "# list of all image numbers\n",
    "all_image_numbers = [1, 2, 3, 4, 5, 6]\n",
    "\n",
    "# TODO: choose 3 entries randomly for the training set\n",
    "train_image_numbers = \n",
    "\n",
    "# TODO: select other 3 for the test set. Make sure there is no overlap\n",
    "test_image_numbers = \n",
    "        \n",
    "# optional convert list into array\n",
    "test_image_numbers = np.array(test_image_numbers)\n",
    "\n",
    "print(train_image_numbers)\n",
    "print(test_image_numbers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ea79971",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_label_per_region(segmented_image, label_map):\n",
    "    \"\"\"\n",
    "    Returns a 1D numpy array that contains the label for each region, shape: (num_regions)\n",
    "            For each region, we obtain the label that has the largest intersection with it\n",
    "    \"\"\"\n",
    "    region_id_list = np.unique(segmented_image)\n",
    "    label_id_list = np.unique(label_map)\n",
    "    region_labels = []\n",
    "    for region_id in region_id_list:\n",
    "        mask_region = segmented_image == region_id\n",
    "        \n",
    "        intersection_per_label = []\n",
    "        for label_id in label_id_list:\n",
    "            mask_label = label_map == label_id\n",
    "            # Compute intersection of each region with each label\n",
    "            intersection = np.sum(mask_region * mask_label)\n",
    "            intersection_per_label.append(intersection)\n",
    "        \n",
    "        intersection_per_label = np.array(intersection_per_label)\n",
    "        # Obtain the index of the label with largest intersection\n",
    "        selected_label = np.argmax(intersection_per_label)\n",
    "        region_labels.append(selected_label)\n",
    "    \n",
    "    return np.array(region_labels).astype(np.uint32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67331b1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm.auto import tqdm # progress bar\n",
    "\n",
    "# Create arrays of training targets and features \n",
    "all_train_region_features = []\n",
    "all_train_region_labels = []\n",
    "\n",
    "for image_number in tqdm(train_image_numbers):\n",
    "    \n",
    "    # TODO: read slic-segmented image (named regions_{image_number}.tif)\n",
    "    segmented_image = \n",
    "    \n",
    "    # TODO: read ground truth image\n",
    "    gt_image = \n",
    "    \n",
    "    # TODO: Get labels per region using the function \"get_label_per_region\" defined above\n",
    "    region_labels = \n",
    "    \n",
    "    # Add current region labels to the variable all_train_region_labels\n",
    "    all_train_region_labels.append(region_labels)\n",
    "    \n",
    "    # TODO: read features (features_{image_number}.npy files) using the function: np.load(file_path.npy)\n",
    "    region_features = \n",
    "    \n",
    "    # Add current region features to the variable all_train_region_features\n",
    "    all_train_region_features.append(region_features)\n",
    "\n",
    "# Tranforming the list all_train_region_labels in an array of shape: (num_all_regions)\n",
    "train_labels = np.concatenate(all_train_region_labels)\n",
    "print(\"train_labels shape \" + str(train_labels.shape))\n",
    "# Tranforming the list all_train_region_features in an array of shape: (num_all_regions, num_features)\n",
    "train_features = np.concatenate(all_train_region_features)\n",
    "print(\"train_features shape \" + str(train_features.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "081a289d",
   "metadata": {},
   "source": [
    "## 4. Normalize features\n",
    "\n",
    "Normalize features by substracting the mean and dividing by the standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "002678b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_per_feature = np.mean(train_features, axis=0)\n",
    "std_per_feature = np.std(train_features, axis=0)\n",
    "# TODO: normalize features by substracting the mean and dividing by the standard deviation\n",
    "norm_train_features = "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50857f70",
   "metadata": {},
   "source": [
    "## 5. Train Random Forest classifier\n",
    "\n",
    "Train model using the ***sklearn*** library. More information about this classifier [here](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23aebb08",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "# TODO: create random forest classifier with the parameter random_state=10\n",
    "classifier = \n",
    "# TODO: fit the model with the normalized features\n",
    "classifier.fit(\"to fill\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62cf8c9b",
   "metadata": {},
   "source": [
    "## 6. Predict classification map for test images\n",
    "\n",
    "For each image in the test set\n",
    "- Load features\n",
    "- Predict labels by region\n",
    "- Predict labels for each pixel\n",
    "- Save the predictions on disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "341b69b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for image_number in test_image_numbers:\n",
    "    # TODO: read segmented image\n",
    "    segmented_image = \n",
    "    \n",
    "    # TODO: read features using the function: np.load(file_path.npy)\n",
    "    region_features = \n",
    "    # TODO: normalize features by substracting the mean and dividing by the standard deviation (of the train set) \n",
    "    norm_region_features = \n",
    "    \n",
    "    # TODO: predict label of regions \n",
    "    label_predictions = classifier.predict(\"todo\")\n",
    "    \n",
    "    # Compute label predictions per pixel\n",
    "    prediction_map = np.zeros(segmented_image.shape).astype(np.uint8)\n",
    "    region_id_list = np.unique(segmented_image)\n",
    "    for idx, region_id in enumerate(region_id_list):\n",
    "        # TODO: set the label predicted for a region to all the pixels of that region\n",
    "        prediction_map[segmented_image==region_id] = \n",
    "        \n",
    "    # Save prediction map image\n",
    "    prediction_map_path = f\"images_week4/prediction_map_{image_number}.tif\"\n",
    "    imsave(prediction_map_path, prediction_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4d61bb9",
   "metadata": {},
   "source": [
    "## 7. Visualize predictions\n",
    "\n",
    "Visualize one of the classification predictions save one disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c61147a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def display_label_image(label_map):\n",
    "    colors = np.array([[255, 255, 255], # Impervious: white\n",
    "                       [0, 0, 255],     # Building: blue\n",
    "                       [0, 255, 255],   # Low vegetation: cyan\n",
    "                       [0, 255, 0]])    # Tree: green\n",
    "    colors = colors.astype(np.uint8)\n",
    "    color_map = np.zeros((label_map.shape[0], label_map.shape[1], 3)).astype(np.uint8)\n",
    "    for label_id in range(colors.shape[0]):\n",
    "        color_map[label_map == label_id] = colors[label_id, :]\n",
    "    plt.figure(figsize = (11,11)) \n",
    "    plt.imshow(color_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf6eac9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take first test image\n",
    "test_image_number = test_image_numbers[0]\n",
    "\n",
    "# TODO: display one of the prediction map images save on disk (using the function display_label_image, defined above)\n",
    "example_prediction = \n",
    "display_label_image(example_prediction)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "367f6d86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Display the ground truth corresponding to the image you visualized above\n",
    "example_gt = \n",
    "display_label_image(example_gt)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe500bf9",
   "metadata": {},
   "source": [
    "## 8. Answer the question and discuss with your neighbor\n",
    "\n",
    "> Comparing the prediction and the ground truth, which types of errors do you find in the predictions?\n",
    "\n",
    "Tip: google for Type 1/Type 2 errors if you think it is too straightforward."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b17f004b",
   "metadata": {},
   "source": [
    "## Congratulations. You Finished! See you next week"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "IPEO",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}