02_features_exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After segments have been created for a given image, it is pivotal to characterize them in a meaningful manner. The power of OBIA partly relies on the fact that a variety of descriptors can be used to characterize segments beyond their spectral properties. Quite often shape, textural or contextual features are beneficial for a classification of the segments. The following notebook serves to explore some basic possibilities to calculate such features.\n",
    "\n",
    "Note that term \"features\" as it used currently, is somewhat ambiguous or imprecise. From the point of the computer vision domain, what we are actually interested in are rather descriptors than features. To be precise: Features are commonly understood to be properties/parts of the image that are of particular interest (-> edges, corners, lines, keypoints, etc), whereas descriptors are used to characterize the features or the image as a whole (-> local vs. global descriptors). However, in the context of machine learning features are input variables used to characterize samples and since the descriptors are often serving this purpose, the exact differentiation of both terms will not be considered further in the following.\n",
    "\n",
    "When it comes to calculating features in python, a way needs to be found to aggregate the information of all pixels for a given segment according to the definition of the specific feature of interest. Depending on the chosen framework, there are different ways to do so. Using the data cube friendly (rio-)xarray package, one can rely on pandas-like groupby operations, for example. Usually somewhat more performant is the `regionprops()` function from [skimage](https://scikit-image.org/docs/stable/api/skimage.measure.html#skimage.measure.regionprops_table). It renders the need for manual grouping obsolete and comes along with some already pre-implemented features. Therefore, we will mainly focus on this one below."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to run this on google colab\n",
    "if 'google.colab' in str(get_ipython()):\n",
    "    import os\n",
    "    repo_dir = \"obia_tutorials\"\n",
    "    marker_file = os.path.join(repo_dir, \".setup_done\")    \n",
    "    if not os.path.exists(marker_file):\n",
    "        !git clone https://github.com/fkroeber/obia_tutorials.git\n",
    "        !pip install -r obia_tutorials/requirements.txt\n",
    "        with open(marker_file, 'w') as f:\n",
    "            f.write(\"Setup completed\")\n",
    "    if not os.getcwd().endswith(repo_dir):\n",
    "        os.chdir(repo_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import matplotlib.colors as mcolors\n",
    "import matplotlib.pyplot as plt\n",
    "import mahotas as mh\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import plotly.express as px\n",
    "import rioxarray as rxr\n",
    "import xarray as xr\n",
    "import time\n",
    "\n",
    "from matplotlib import rcParams\n",
    "from scipy import stats\n",
    "from skimage import graph\n",
    "from skimage.color import label2rgb\n",
    "from skimage.measure import regionprops, regionprops_table\n",
    "from skimage.morphology import binary_dilation\n",
    "from skimage.segmentation import mark_boundaries, slic\n",
    "from skimage.util import map_array\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.model_selection import GridSearchCV, RepeatedKFold\n",
    "from tqdm import tqdm\n",
    "\n",
    "# set general figsize\n",
    "rcParams['figure.figsize'] = (7.5, 7.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exemplary data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For demonstration purposes, we will use a 4-band R,G,B,NIR image segmented with SLIC. The segmentation parameters below were chosen by trial-and-error until objects of interest (e.g. trees and their shadows) appeared to be segmented in a meaningful manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read image\n",
    "img_path = \"sample_data/ortho_subset_I.tif\"\n",
    "ds = xr.open_dataset(img_path)\n",
    "\n",
    "# extract arrays\n",
    "arr = np.transpose(np.array(ds.to_array()[0]), (1,2,0))\n",
    "rgb = arr.astype(np.uint8)[...,[0,1,2]]\n",
    "ndvi = (arr[...,3] - arr[...,0])/(arr[...,3] + arr[...,0])\n",
    "\n",
    "# perform segmentation\n",
    "segments = slic(arr, n_segments=3000, compactness=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 1: Print the shape of the image and the number of segments.</span>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise rgb & segmentation\n",
    "fig, axs = plt.subplots(ncols=3, figsize=(15,5), constrained_layout=True)\n",
    "axs[0].imshow(rgb)\n",
    "axs[1].imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "axs[2].imshow(label2rgb(segments, rgb, kind='avg').astype(int))\n",
    "\n",
    "axs[0].set_title(\"original rgb\")\n",
    "axs[1].set_title(\"segmentation boundaries\")\n",
    "axs[2].set_title(\"mean rgb per segment\")\n",
    "\n",
    "for ax in axs:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this areal image is a very high-resolution image and quite large (1500 x 1500 pxls). This allows to investigate the scalability of the calculation of the proposed feature descriptors. To this end, a decorator for timing the execution of some subsequently defined functions is defined. If you want to get the timing more accurately, you can put the magic `%%timeit` at the beginning of a chuck to repeat the cell execution multiple times and get an average. For a rough estimate, however, the function decorator below should be sufficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# decorator for timing\n",
    "def timeit(method):\n",
    "    def timed(*args, **kw):\n",
    "        tstart = time.perf_counter()\n",
    "        result = method(*args, **kw)\n",
    "        tend = time.perf_counter()\n",
    "        total_ms = (tend - tstart)\n",
    "        print(f\"Execution time: {total_ms:.2f}s\")\n",
    "        return result   \n",
    "    return timed"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Types of feature descriptors"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### I. Spectral features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Very intuitive and fast to calculate are the measures of central tendency and variation, the most popular of which are the mean spectral value for a segment and its standard deviation. The calculation of the mean value is pretty straightforward and was already carried out implicitly by using the function `label2rgb` above. For non-visualisation related purposes, the pre-implemented mean calculation provided by `regionprops` can be used as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define mean calculation function\n",
    "@timeit\n",
    "def calc_mean(seg_arr, img_arr):\n",
    "    spec_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        intensity_image = img_arr,\n",
    "        properties = [\"label\", \"intensity_mean\"]\n",
    "        )\n",
    "    return pd.DataFrame(spec_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "means = calc_mean(segments, np.dstack([rgb, ndvi]))\n",
    "display(means)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the calculation of the standard deviation we can write a custom extra_property function to be passed to regionprops. These functions need to take a regionmask and (optionally) an intensity mask as input and define a function that transforms these 2D array into an output value for a given segment. It is important to note that for many intensity-based calculations, the regionmask (which is a binary array) needs to be used to first filter all intensity values that actually belong to the segment for the given intensity mask. Otherwise, the calculation is performed on the intensity array that represents the bounding box rather than the mask of the segment. The example for the standard deviation demonstrates the basic usage of the extra_properties.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 2: Implement the calculation of the standard deviation as an extra_property below.</span>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define stdv function\n",
    "def std(regionmask, intensity_img):\n",
    "    # code goes here \n",
    "    return std\n",
    "    \n",
    "@timeit\n",
    "def calc_std(seg_arr, img_arr):\n",
    "    spec_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        intensity_image = img_arr,\n",
    "        properties = [\"label\",],\n",
    "        extra_properties=(std,)\n",
    "        )\n",
    "    return pd.DataFrame(spec_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "stds = calc_std(segments, np.dstack([rgb, ndvi]))\n",
    "display(stds)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that at the moment, regionprops doesn't support reducing the information across bands. Instead, all the intensity-specific functions are applied channel-wise. Since band-wise information is often desired, this is not too much of a problem. In the example above, calculating the std across r,g,b bands and nir wouldn't be very reasonable due to the different scale of the single variables. For the sake of completeness, however, the following shows one way of calculating the standard deviation across bands (using only rgb). For this purpose, a transformation to xarrays followed by grouping according to the segmentation layer is used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define stdv across bands function\n",
    "@timeit\n",
    "def calc_std_3d(seg_arr, img_arr):\n",
    "    seg_arr = xr.DataArray(seg_arr)\n",
    "    img_arr = xr.DataArray(img_arr)\n",
    "    stds = img_arr.groupby(seg_arr).std(dim=...)\n",
    "    return pd.DataFrame({\"label\": stds.group, \"std\": np.array(stds)})\n",
    "\n",
    "# perform calculation & get timing\n",
    "stds_3d = calc_std_3d(segments, rgb)\n",
    "display(stds_3d)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the results there is a helper function called [skimage.util.map_array](https://scikit-image.org/docs/stable/api/skimage.util.html#skimage.util.map_array). \n",
    "\n",
    "**<span style=\"color:orange\">Task 3: Take a look at the documentation of the function map_array to visualise the results that we've calculated above.</span>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping NDVI vals back to image segments\n",
    "\n",
    "# code goes here\n",
    "# mapped_mean = ...\n",
    "# mapped_std = ...\n",
    "\n",
    "# visualise results\n",
    "fig, axs = plt.subplots(ncols=4, figsize=(20,10), constrained_layout=True)\n",
    "axs[0].imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "axs[1].imshow(ndvi, cmap=\"YlGn\")\n",
    "axs[2].imshow(mapped_mean, cmap=\"YlGn\")\n",
    "axs[3].imshow(mapped_std)\n",
    "\n",
    "axs[0].set_title(\"rgb & boundaries\")\n",
    "axs[1].set_title(\"ndvi per pixel\")\n",
    "axs[2].set_title(\"mean ndvi per segment\")\n",
    "axs[3].set_title(\"std ndvi per segment\")\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mapping the obtained values back to the segments array shows that - as expected - pure road segments have very low NDVI values and variations. However, segments that include several objects, such as cars parked near trees, stand out as inhomogeneous. Vegetation segments have strong average NDVI values and a medium to high variability. The water surface appears as a very homogenous region with hardly any differences in mean and variation of NDVI among segments belonging to this group. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### II. Shape features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Starting with some basic shape features, one may want to calculate the area of the bounding box of segments to get an indication for their extent. Given the bounding box one can further calculate features such as rectangularity defined as the ratio between a segment's area and its bbox.        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 4: Implement the calculation of the rectangularity.</span>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define rectangularity\n",
    "def rectangularity(regionmask):\n",
    "    # code goes here\n",
    "    \n",
    "@timeit\n",
    "def calc_rect_shapes(seg_arr):\n",
    "    shp_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        properties = [\"label\", \"area_bbox\"],\n",
    "        extra_properties=(rectangularity,)\n",
    "        )\n",
    "    return pd.DataFrame(shp_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "rects = calc_rect_shapes(segments)\n",
    "display(rects)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that rectangularity as defined above isn't rotation-invariant, i.e. it value changes depending on the orientation of a segment. Since many natural objects can have various orientations while still being the same object (e.g. parked cars), it is therefore rarely a meaningful feature. \n",
    "\n",
    "Conceptually stronger are features such as compactness or circularity, which relates the area of an object to its perimeter (called shape-index in eCognition). Another sometimes useful feature is solidity defined as the ratio of pixels in a segment to pixels of its convex hull. The later one is also pre-implemented in regionprops."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define compactness\n",
    "@timeit\n",
    "def calc_shapes(seg_arr, solidity=True):\n",
    "    if solidity:\n",
    "        props = [\"label\", \"solidity\", \"area\", \"perimeter\"]\n",
    "    else:\n",
    "        props = [\"label\", \"area\", \"perimeter\"]\n",
    "    shp_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        properties = props,\n",
    "        )\n",
    "    compactness = 4*np.pi*shp_feats[\"area\"]/(shp_feats[\"perimeter\"]**2)\n",
    "    shp_feats[\"compactness\"] = compactness\n",
    "    shp_feats.pop(\"area\")\n",
    "    shp_feats.pop(\"perimeter\")\n",
    "    return pd.DataFrame(shp_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "shapes = calc_shapes(segments)\n",
    "display(shapes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mapping vals back to image segments\n",
    "mapped_compactness = map_array(segments, np.array(shapes[\"label\"]), np.array(shapes[\"compactness\"]))\n",
    "mapped_solidity = map_array(segments, np.array(shapes[\"label\"]), np.array(shapes[\"solidity\"]))\n",
    "\n",
    "# visualise stdvs\n",
    "fig, axs = plt.subplots(ncols=3, figsize=(15,10), constrained_layout=True)\n",
    "axs[0].imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "axs[1].imshow(mapped_compactness)\n",
    "axs[2].imshow(mapped_solidity)\n",
    "\n",
    "axs[0].set_title(\"rgb & boundaries\")\n",
    "axs[1].set_title(\"compactness\")\n",
    "axs[2].set_title(\"solidity\")\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 5: Which feature descriptor (compactness vs. solidity) seems favourable to you? Do they differ in terms of computational costs? Run the calculation with and without solidity to evaluate the difference.</span>**\n",
    "\n",
    "<i>Answer goes here</i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Talking about scalability, the family of the so-called [image moments](https://en.wikipedia.org/wiki/Image_moment) should be mentioned as an additional set of computationally favorable descriptors. As weighted averages of the image pixels' intensities they provide information about the orientation and symmetry of a segment. One important formulation of image moments are the Hu moments. They extend the basic central moment definitions in a way that the resulting set of features is invariant to certain transformations. These transformations include translation (movement), rotation, and scaling (size changes), which makes them powerful feature descriptors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 6: Insert the corresponding property for the Hu moments below.</span>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timeit\n",
    "def calc_moments(seg_arr):\n",
    "    shp_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        # code goes here\n",
    "        # properties = [\"label\", ...],\n",
    "        )\n",
    "    return pd.DataFrame(shp_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "moments = calc_moments(segments)\n",
    "display(moments)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise hu moments\n",
    "fig, axs = plt.subplots(ncols=4, nrows=2, figsize=(20,10), constrained_layout=True)\n",
    "\n",
    "# mapping vals back to image segments\n",
    "for i, ax in enumerate(axs.flat):\n",
    "    if i == 0:\n",
    "        ax.imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "        ax.set_title(\"rgb & boundaries\")\n",
    "    else:\n",
    "        mapped_moments = map_array(segments, np.array(moments[\"label\"]), np.array(moments.iloc[:,i]))\n",
    "        ax.imshow(mapped_moments)\n",
    "        ax.set_title(moments.iloc[:,i].name)\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The interpretation of the single moments is somewhat difficult due to their missing direct correspondence to physical aspects of the shapes. They are rather constructed in a mathematical way to satisfy the above-mentioned invariance properties. One can generalize that the 1st Hu moment (moments_hu_0) represents the overall symmetry of the object around its centre of mass. It's sensitive to elongation, meaning that elongated shapes will have a higher value for this moment. This sensitivity to elongation is also inherent to the 2nd Hu moment which is related to the object's symmetry along the major and minor axes (note that the shape descriptors `axis_major_length` and `axis_minor_length` provided by skimage are therefore based on image moment calculation under the hood). Higher order moments often relate to other aspects of the segments' shapes such as skewness and kurtosis. As we are dealing with a superpixel segmentation, the shapes are rather symmetrical and therefore less prominent in the corresponding visualizations.  "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### III. Textural features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Textural features are commonly calculated globally, i.e. for the whole image. However, there is nothing that prevents us from calculating them for each segment individually. Famous texture descriptors are the Haralick features which are based on the 2D histogram distributions of the intensity differences between neighboring pixels. These so-called grey-level co-occurrence matrices (GLCM) can be used to derive a set of descriptors such as the angular second moment, contrast, correlation,...\n",
    "\n",
    "Computing the Haralick features locally for segments requires some pre-processing by buffering the segments to avoid errors in case of thin segments that don't have any neighbor in a certain direction. Also note that the Haralick features are based on single band greyscale images and require them to be non-negative integers - if we want to calculate them for the NDVI, we thus have to rescale the values first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timeit\n",
    "def calc_haralick_feats(seg_arr, img_arr):\n",
    "    \"\"\"\n",
    "    Calculates the 13 stable Haralick features for image segments,\n",
    "    based on mahotas as more efficient & comprehensive implementation compared to skimage\n",
    "    \"\"\"\n",
    "    seg_arr = seg_arr.astype(np.uint16)\n",
    "    img_arr = img_arr.astype(np.uint8)\n",
    "    regions = regionprops(seg_arr)\n",
    "    region_props = {}\n",
    "    for region in tqdm(regions):\n",
    "        # create dilated region mask to avoid boundary effects\n",
    "        mask = seg_arr == region.label\n",
    "        mask = binary_dilation(mask)\n",
    "        # clip image to bbox\n",
    "        [rows, cols] = np.where(mask)\n",
    "        row1, row2 = min(rows), max(rows)+1\n",
    "        col1, col2 = min(cols), max(cols)+1\n",
    "        mask_clipped = mask[row1:row2, col1:col2]\n",
    "        img_clipped = img_arr[row1:row2, col1:col2]\n",
    "        # fill background with zeros\n",
    "        clip = np.where(mask_clipped, img_clipped, 0)\n",
    "        # calculate texture features for different orientations & average\n",
    "        feats = mh.features.haralick(clip, ignore_zeros=True).mean(0)\n",
    "        region_props[region.label] = feats\n",
    "    # summarise results\n",
    "    haralick_df = pd.DataFrame(region_props).T\n",
    "    haralick_cols = [f\"haralick-{x}\" for x in np.arange(haralick_df.shape[1])] \n",
    "    haralick_df.columns = haralick_cols\n",
    "    return haralick_df\n",
    "\n",
    "# perform calculation & get timing\n",
    "haralicks = calc_haralick_feats(segments, (100*(ndvi+1)).astype(np.uint8))\n",
    "display(haralicks)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 7: Why did this calculation take so long? Can you think of other textural descriptors that may be computationally cheaper?</span>**\n",
    "\n",
    "<i>Answer goes here</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us take a look at some of the results and see if it was worth it to perform this expensive calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise results\n",
    "mapped_asm = map_array(segments, np.unique(segments), np.array(haralicks[\"haralick-0\"]))\n",
    "mapped_contrast = map_array(segments, np.unique(segments), np.array(haralicks[\"haralick-1\"]))\n",
    "mapped_correlation = map_array(segments, np.unique(segments), np.array(haralicks[\"haralick-2\"]))\n",
    "\n",
    "fig, axs = plt.subplots(ncols=4, figsize=(20,10), constrained_layout=True)\n",
    "axs[0].imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "axs[1].imshow(mapped_asm)\n",
    "axs[2].imshow(mapped_contrast)\n",
    "axs[3].imshow(mapped_correlation)\n",
    "\n",
    "axs[0].set_title(\"rgb & boundaries\")\n",
    "axs[1].set_title(\"angular second moment\")\n",
    "axs[2].set_title(\"contrast\")\n",
    "axs[3].set_title(\"correlation\")\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:orange\">Task 8: Describe and interpret the results briefly. Do these Haralick features over some additional information compared to the previously calculated ones?</span>**\n",
    "\n",
    "<i>Answer goes here</i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "If we want to find a cheaper substitution for this textural descriptor, we may look for first-order metrics. They operate on the greyscale values and their distribution directly rather than analyzing similarities between pixel pairs. One of these first-order metrics, the variance, was already introduced in the section of spectral features. It inherently conveys information on the magnitude of greyscale differences for a given segment. Another popular distribution-based metric for the degree of disorder of values is the entropy measure as shown subsequently. Note that the entropy measure is based on the bins of the greyscale histogram and therefore either requires a common scaling of input bands or a binning schema that fits the specific band of interest.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy_ndvi(regionmask, intensity_img):\n",
    "    vals = intensity_img[regionmask]\n",
    "    arr = stats.relfreq(vals, 100, defaultreallimits=(-1,1))[0]\n",
    "    return stats.entropy(arr)\n",
    "\n",
    "@timeit\n",
    "def calc_entropy(seg_arr, img_arr):\n",
    "    entropy_feats = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        intensity_image = img_arr, \n",
    "        properties = [\"label\"],\n",
    "        extra_properties = (entropy_ndvi,)\n",
    "        )\n",
    "    return pd.DataFrame(entropy_feats)  \n",
    "\n",
    "# perform calculation & get timing\n",
    "entropy = calc_entropy(segments, ndvi)\n",
    "display(entropy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise results\n",
    "mapped_entropy = map_array(segments, np.unique(segments), np.array(entropy[\"entropy_ndvi\"]))\n",
    "\n",
    "fig, axs = plt.subplots(ncols=2, figsize=(7.5,7.5), constrained_layout=True)\n",
    "axs[0].imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"thick\"))\n",
    "axs[1].imshow(mapped_entropy)\n",
    "\n",
    "axs[0].set_title(\"rgb & boundaries\")\n",
    "axs[1].set_title(\"entropy\")\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### IV. Contextural features"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Contextual features concern the neighborhood of a given segment and are particularly interesting so they allow the characterization of segments with a larger effective field of view. To calculate such features the representation of a segmentation as a so-called Region Adjacency Graph (RAG) can be leveraged. In this graph-based representation, the segments are represented by nodes connected to their neighbors via edges. Skimage provides RAG functionality based on the underlying Networkx library.  \n",
    "\n",
    "To demonstrate the usage of contextual features we will focus on the spectral value of the neighbors. However, the following examples can easily be extended to other features of the neighboring segments (such as shape features).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timeit\n",
    "def calc_mean_std_neighbours(seg_arr, img_arr):\n",
    "    # calculate region adjacency graph\n",
    "    rag = graph.RAG(seg_arr)\n",
    "    # get properties of neighbours that are of interest\n",
    "    neighbor_attrs = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        intensity_image = img_arr,\n",
    "        properties = [\"label\", \"mean_intensity\"],\n",
    "        )\n",
    "    neighbor_attrs = pd.DataFrame(neighbor_attrs)\n",
    "    neighbor_attrs.set_index(\"label\", inplace=True)\n",
    "    # accumulate data of neighbouring segments\n",
    "    neighbor_data = []\n",
    "    for node, neighbors in rag.adjacency():\n",
    "        nn_attrs = neighbor_attrs.loc[list(map(int, neighbors))]\n",
    "        nn_means = nn_attrs.mean()\n",
    "        nn_std = nn_attrs.std()\n",
    "        neighbor_data.append([node, len(neighbors), *nn_means, *nn_std])\n",
    "    # compile results\n",
    "    neighbor_attrs_df = pd.DataFrame(neighbor_data).set_index(0)\n",
    "    neighbor_feats = neighbor_attrs_df.fillna(0)\n",
    "    neighbor_feats.index.names = ['label']\n",
    "    n_bands_img_arr = img_arr.shape[-1] if len(img_arr.shape) == 3 else 1\n",
    "    neighbor_feats.columns = [\n",
    "        'neighbours_n',\n",
    "        *[f'neighbours_mean_{i}' for i in np.arange(n_bands_img_arr)],\n",
    "        *[f'neighbours_std_{i}' for i in np.arange(n_bands_img_arr)],\n",
    "        ]\n",
    "    return neighbor_feats.reset_index()\n",
    "\n",
    "# perform calculation & get timing\n",
    "neighb_vals = calc_mean_std_neighbours(segments, rgb)\n",
    "display(neighb_vals)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The real power of a graph-based representation, however, lies in the beauty of graph-based mathematics. They allow more efficient calculations, which render the need of the nasty for-loop used above. In particular, a 2D adjacency matrix describing the neighboring relationship between nodes can be used to calculate the information on neighbors means and stdvs in a faster way by using numpy dot product calculations. Note that the result is the same as above.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@timeit\n",
    "def calc_mean_std_neighbours(seg_arr, img_arr):\n",
    "    # create region adjacency graph & matrix\n",
    "    rag = graph.RAG(seg_arr)\n",
    "    adj_matrix = nx.adjacency_matrix(rag)\n",
    "    # get properties of neighbours that are of interest\n",
    "    neighbor_attrs = regionprops_table(\n",
    "        label_image = seg_arr, \n",
    "        intensity_image = img_arr,\n",
    "        properties = [\"label\", \"mean_intensity\"],\n",
    "        )\n",
    "    neighbor_attrs = pd.DataFrame(neighbor_attrs)\n",
    "    neighbor_attrs.set_index(\"label\", inplace=True)\n",
    "    # get number of neighbours\n",
    "    num_neighbors = adj_matrix.sum(axis=1) \n",
    "    # convert neighbor_attrs to a numpy array\n",
    "    neighbor_attrs_array = neighbor_attrs.reindex(list(rag.nodes)).to_numpy()\n",
    "    # get mean vals\n",
    "    neighbor_sum = adj_matrix.dot(neighbor_attrs_array)\n",
    "    neighbor_mean = neighbor_sum/(num_neighbors.reshape(-1,1))\n",
    "    # get std vals\n",
    "    neighbor_attrs_sq_sum = adj_matrix.dot(np.square(neighbor_attrs_array))\n",
    "    num_neighbors_bessel = np.maximum(num_neighbors - 1, 1)[:, None]\n",
    "    neighbor_attrs_sq_sum_bessel = np.divide(neighbor_attrs_sq_sum, num_neighbors_bessel, where=num_neighbors_bessel != 0)\n",
    "    neighbor_var = (\n",
    "        neighbor_attrs_sq_sum_bessel - \n",
    "        np.divide(np.square(neighbor_sum), num_neighbors[:, None] * num_neighbors_bessel, where=num_neighbors_bessel != 0)\n",
    "        )\n",
    "    neighbor_std = np.sqrt(neighbor_var)\n",
    "    # combine the results\n",
    "    neighbor_data = np.column_stack([rag.nodes, num_neighbors, neighbor_mean, neighbor_std])\n",
    "    neighbor_attrs_df = pd.DataFrame(neighbor_data).set_index(0)\n",
    "    neighbor_feats = neighbor_attrs_df.fillna(0)\n",
    "    neighbor_feats.index.names = ['label']\n",
    "    n_bands_img_arr = img_arr.shape[-1] if len(img_arr.shape) == 3 else 1\n",
    "    neighbor_feats.columns = [\n",
    "        'neighbours_n',\n",
    "        *[f'neighbours_mean_{i}' for i in np.arange(n_bands_img_arr)],\n",
    "        *[f'neighbours_std_{i}' for i in np.arange(n_bands_img_arr)],\n",
    "        ]\n",
    "    return neighbor_feats.reset_index()\n",
    "\n",
    "# perform calculation & get timing\n",
    "neighb_vals = calc_mean_std_neighbours(segments, rgb)\n",
    "display(neighb_vals)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that there are also some contextual features that doesn't require the calculation of a RAG in the first place. Among these descriptors, you can find the [Euler number](https://sydney4.medium.com/the-euler-number-of-a-binary-image-1d3cc9e57caa).\n",
    "\n",
    "**<span style=\"color:orange\">Task 9: Implement the calculation of the Euler number below. Evaluate its value in the given case.</span>**  \n",
    "\n",
    "<i>Answer goes here</i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define stdv across bands function\n",
    "@timeit\n",
    "def calc_euler(seg_arr):\n",
    "    # code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synergy & Use case"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting everything together, you can now calculate a set of more than 20 features to describe the segments holistically. The expressiveness of these features shall be demonstrated below by using them as input for a random forest based classification of objects. In particular, we will try to create a binary classifier to identify the trees in the image. To this end, some samples in the original image are labelled and used together with some negative samples (representing the trees shadows, water, parking lots, etc) as input to a random forest. To evaluate the classifier, we will use a different subset spatially next to the one used for training, segment it analogously to the first one and apply the trained random forest to retrieve the predictions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate set of various features\n",
    "@timeit\n",
    "def calc_all_feats(seg_arr, rgbndvi_arr):\n",
    "    ndvi_arr = rgbndvi_arr[...,-1]\n",
    "    spec_feats = pd.DataFrame(\n",
    "        regionprops_table(\n",
    "            label_image = seg_arr, \n",
    "            intensity_image = rgbndvi_arr,\n",
    "            properties = [\"label\", \"intensity_mean\"]\n",
    "        )\n",
    "    )\n",
    "    shp_feats = pd.DataFrame(\n",
    "        regionprops_table(\n",
    "            label_image = seg_arr, \n",
    "            properties = [\"moments_hu\"],\n",
    "        )\n",
    "    )\n",
    "    text_feats_I = pd.DataFrame(\n",
    "        regionprops_table(\n",
    "            label_image = seg_arr, \n",
    "            intensity_image = rgbndvi_arr,\n",
    "            properties = [],\n",
    "            extra_properties=(std,)\n",
    "        )\n",
    "    )\n",
    "    text_feats_II = calc_entropy(seg_arr, ndvi_arr).iloc[:,1:]\n",
    "    context_feats = calc_mean_std_neighbours(seg_arr, rgbndvi_arr).iloc[:,1:]\n",
    "    all_feats = pd.concat([spec_feats, shp_feats, text_feats_I, text_feats_II, context_feats], axis=1)\n",
    "    return all_feats  \n",
    "\n",
    "# perform calculation & get timing\n",
    "all_feats_df = calc_all_feats(segments, np.dstack([rgb,ndvi]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# samples for both classes\n",
    "trees_idxs = [\n",
    "    1204, 1186, 1134, 1229, 1329, 1201, 1140, 1218, 1351, 962, 830, 886, 926, 1373, 1405, 1517, 1547, 1461, 1526, 1406, 1358, 1460, 1556, 735, 833, 838, 802, 513, 616, 565, 438, 258, 375, 334, 284, 204, 250, 155, 183, 241, 344, 359, 412, 328, 424, 506, 468, 554, 459, 579, 560, 769, 890, 791, 798, 811, 771, 846, 1076, 784, 919, 915, 133, 246, 356, 405, 425, 474, 166, 158, 34, 200, 209, 245, 55, 2003, 1984, 1910, 1813, 1873, 1923, 2296, 2220, 2143, 2113, 2153, 2319, 2275, 2284, 2422, 2451, 2398, 2460, 2455, 2208, 2218, 2199, 2036, 1986, 1924\n",
    "]\n",
    "\n",
    "other_idxs = [\n",
    "    2647, 2653, 2717, 2784, 2110, 2061, 2122, 1988, 1943, 2018, 1861, 1896, 1909, 1900, 1829, 1758, 1719, 1703, 1794, 1812, 1749, 1581, 1580, 1522, 1451, 1438, 1479, 1408, 1511, 1631, 1869, 1811, 1978, 2104, 2187, 2299, 2400, 2732, 2642, 2567, 2348, 2369, 2401, 2281, 1850, 1117, 197, 1002, 872, 216, 432, 418, 357, 382, 419, 336, 295, 351, 460, 644, 708, 800, 902, 1005, 1063, 174, 824, 831, 780, 1291, 1412, 1537, 1497, 1570, 1605, 1491, 1595, 1652, 1768, 1807, 1887, 1898, 1968, 1573, 1540, 1536, 1439, 1388, 1694, 2640 \n",
    "]\n",
    "\n",
    "# map sample_idxs to 2D grid\n",
    "trees_bool = [2 if x in trees_idxs else 0 for x in np.array(means[\"label\"])]\n",
    "other_bool = [1 if x in other_idxs else 0 for x in np.array(means[\"label\"])]\n",
    "both_labels = [t if t else o for t,o in zip(trees_bool, other_bool)]\n",
    "\n",
    "samples_df = pd.DataFrame({\n",
    "    \"id\": np.array(means[\"label\"]),\n",
    "    \"label\": both_labels\n",
    "})\n",
    "\n",
    "mapped_samples = map_array(\n",
    "    segments, \n",
    "    np.array(samples_df[\"id\"]), \n",
    "    np.array(samples_df[\"label\"])\n",
    "    )\n",
    "\n",
    "# visualise drawn samples\n",
    "background_c, others_c, trees_c = '#7f8282', '#00008b', '#00FF00'\n",
    "binary_cmap = mcolors.ListedColormap([background_c, others_c, trees_c])\n",
    "plt.imshow(mark_boundaries(rgb, segments, (1,0,0), mode=\"outer\"))\n",
    "plt.imshow(mapped_samples, cmap=binary_cmap, alpha=0.5)\n",
    "plt.axis(\"off\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract feature data for all samples\n",
    "tree_samples = all_feats_df[all_feats_df[\"label\"].isin(trees_idxs)]\n",
    "tree_samples.insert(0, \"gt\", 1)\n",
    "other_samples = all_feats_df[all_feats_df[\"label\"].isin(other_idxs)]\n",
    "other_samples.insert(0, \"gt\", 0)\n",
    "feats_samples = pd.concat([tree_samples, other_samples]).drop(columns=\"label\")\n",
    "feats_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate random forest & grid search hyperparameter tuning\n",
    "rf = RandomForestClassifier(\n",
    "    n_estimators=100,\n",
    "    random_state=42,\n",
    "    oob_score=True\n",
    ")\n",
    "grid_search = GridSearchCV(\n",
    "    rf,\n",
    "    {\"max_features\": np.arange(0.05, 0.55, 0.05)},\n",
    "    scoring=\"accuracy\",\n",
    "    cv=RepeatedKFold(n_splits=3, n_repeats=5, random_state=42),\n",
    "    n_jobs=-1,\n",
    "    verbose=1\n",
    ")\n",
    "\n",
    "# fit random forest & evaluate quantitatively\n",
    "X = feats_samples.iloc[:,1:]\n",
    "y = feats_samples.iloc[:,0]\n",
    "grid_search.fit(X, y)\n",
    "preds = grid_search.predict(all_feats_df.iloc[:,1:])\n",
    "\n",
    "print(f\"Out-of-bag accuarcy RF: {grid_search.best_estimator_.oob_score_:.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read & segment test image\n",
    "img_path_test = \"sample_data/ortho_subset_II.tif\"\n",
    "ds_test = xr.open_dataset(img_path_test)\n",
    "arr_test = np.transpose(np.array(ds_test.to_array()[0]), (1,2,0))\n",
    "rgb_test = arr_test.astype(np.uint8)[...,[0,1,2]]\n",
    "ndvi_test = (arr_test[...,3] - arr_test[...,0])/(arr_test[...,3] + arr_test[...,0])\n",
    "segments_test = slic(arr_test, n_segments=3000, compactness=0.3)\n",
    "\n",
    "# calculate features for test set\n",
    "feats_df_test = calc_all_feats(segments_test, np.dstack([rgb_test, ndvi_test]))\n",
    "\n",
    "# get predictions for test set\n",
    "preds_test = grid_search.predict(feats_df_test.iloc[:,1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualise final results\n",
    "mapped_preds = map_array(segments, np.unique(segments), preds)\n",
    "mapped_preds_test = map_array(segments_test, np.unique(segments_test), preds_test)\n",
    "\n",
    "fig, axs = plt.subplots(ncols=4, figsize=(20,10), constrained_layout=True)\n",
    "axs[0].imshow(rgb)\n",
    "axs[1].imshow(mapped_preds, cmap=binary_cmap, interpolation=\"nearest\")\n",
    "axs[2].imshow(rgb_test)\n",
    "axs[3].imshow(mapped_preds_test, cmap=binary_cmap, interpolation=\"nearest\")\n",
    "\n",
    "axs[0].set_title(\"train image\")\n",
    "axs[1].set_title(\"trees prediction\")\n",
    "axs[2].set_title(\"test image\")\n",
    "axs[3].set_title(\"trees prediction\")\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set_axis_off();"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the model's predictions are far from perfect but not too bad for a model that has been trained on just a handful of samples and features (less than 200 and 25, respectively). Note also that a pixel-wise classification would never achieve such accuracies and that deep learning approaches are more accurate but usually need hundreds or thousands of training images - our small model here runs on 1 sparsely annotated training image ;)        "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion & Further thoughts"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial has showcased possibilities to calculate various feature descriptors. While it can be generalized that some features are commonly more expressive than others, the exact choice of the descriptors depends on the segmentation and the use case at hand. Other important considerations include the... \n",
    "* interpretability of features - some having a very intuitive meaning while others are less interpretable    \n",
    "* scalability of calculations - esp. in the era of big data relevant as soon as we want to apply models to larger amounts of data\n",
    "\n",
    "Further directions that may be interesting to explore are larger scale experiments on which feature sets appear to be most expressive for certain use cases. Training machine learning models such as random forests on different feature sets for different datasets and identifying the smallest, computationally most efficient feature set that allows to achieve high accuracies for a given task may be of broader interest especially in the context of the rise of deep learning models. While training deep learning models requires a lot of sample data and takes a relatively long time, the inference time (i.e. the time to apply the model to unseen samples) is usually very low. The development of a competitive non-deep learning based OBIA method therefore requires more than ever a focus on small, meaningful and efficiently computable feature sets in order to keep up.    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "obia_tutorials",
   "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.10.1"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "03ca474701badda27714bc2f8463045a15615715ede98ee8a907c42ed944622a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}