opencv_practical.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reality Capture & Precision 3D Scanning: Introduction to OpenCV in python\n",
    "#### David Griffiths &nbsp;&nbsp;&nbsp; david.griffiths.16@ucl.ac.uk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Google colab compatability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Only run if using Google Colab\n",
    "!git clone https://github.com/3Dimaging-ucl/ucl_cege0092 data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What is an image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "\n",
    "img = np.random.randint(low=0, high=255, size=(100,100, 3))\n",
    "plt.imshow(img), plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic OpenCV functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Import OpenCV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "# Set image directory '/data/images' for google colab\n",
    "image_dir = 'images/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Load image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = image_dir+'aerial_small.tif'\n",
    "\n",
    "# Load image in BGR colorspace\n",
    "img_color = cv2.imread(file, 1)\n",
    "#Load image in grayscale\n",
    "img_gray = cv2.imread(file, 0)\n",
    "# Convert BGR to RGB\n",
    "img_color_RGB = cv2.cvtColor(img_color, cv2.COLOR_BGR2RGB)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the images we can create a small helper function. Here we take a list of images, their respective titles, the size of the figure and the subplot shape as input arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_images(images, titles, figsize, shape, ticks=False):\n",
    "    ''' Image plotter helper function.\n",
    "        \n",
    "        parameters:\n",
    "            images  : list\n",
    "            titles  : titles\n",
    "            figsize : tuple\n",
    "            shape   : tuple\n",
    "                      rows x cols\n",
    "            ticks   : Bool\n",
    "                      Set True to show axis ticks\n",
    "        returns:\n",
    "            None\n",
    "    '''\n",
    "    plt.figure(figsize=figsize)\n",
    "    \n",
    "    for i in range(len(images)):\n",
    "        plt.subplot(shape[0],shape[1],i+1),plt.imshow(images[i],'gray')\n",
    "        plt.title(titles[i])\n",
    "        if ticks == False:\n",
    "            plt.xticks([]),plt.yticks([])\n",
    "    plt.show()\n",
    "    \n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot images using helper function\n",
    "plot_images([img_color, img_gray, img_color_RGB], \n",
    "            ['BGR Color Image', 'Grayscale Image', 'RGB Color Image'], \n",
    "            (20,20),\n",
    "            (1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Read and edit pixels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenCV images are stored in the form of numpy arrays. The numpy array has the dimensions (image height x image width x number of channels). The images can therefore be treated as standard numpy $n$ dimension arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print ('One channel grayscale image: ', img_gray.shape)\n",
    "print ('3 Chanel BGR colorspace image: ', img_color.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access a specific pixel value basic numpy array indexing can be used. For example, to return the intensity of pixel x=500, y= 1000 you can use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel = 500, 750\n",
    "\n",
    "print (img_gray[pixel])\n",
    "print (img_color[pixel])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how the colour image has 3 channels and therefore three pixel intensities (Blue (0), Green (1), Red (2)). Further indexing can reutrn a specific intensity i.e:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "blue = img_color[pixel][0]\n",
    "green = img_color[pixel][1]\n",
    "red = img_color[pixel][2]\n",
    "print (blue, green, red)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Task**: How would we create a simple 'for' loop to loop through every pixel in a given image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try task here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Extra:** How could we use numpy indexing to remove the need for a 'for' loop?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try task here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Challenge:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How could we apply the same principles for a RGB image with 3 channels? Load in the image `newyork.jpg` and mask out all pixels except the park using thresholding. Think... what makes these pixels different from the others. \n",
    "\n",
    "***Hint*** masked pixels have a value of 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread(image_dir+'newyork.jpg', 1)\n",
    "img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n",
    "plt.figure(figsize=(20,20)), plt.imshow(img), plt.show()\n",
    "\n",
    "# Try the rest of the challenge here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometric Transformations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When carrying out geometric transformations of an image it is important to think of the image as a standard matrix. This way typical linear algebra can be used to perform transformations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we will import a new image which which will allow us to easily see the result of the transformations.\n",
    "\n",
    "img = cv2.imread(image_dir+'circle.jpg', 1)\n",
    "img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n",
    "plt.imshow(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenCV offers a resize function which can be called using cv2.resize. \n",
    "# We will go through a few ways to use this function.\n",
    "\n",
    "rows, cols, channels = img.shape\n",
    "\n",
    "# Stretch collumns by a factor of 2\n",
    "img_resize = cv2.resize(img, (cols*2, rows), interpolation = cv2.INTER_CUBIC)\n",
    "\n",
    "# Notice the interpolation method here. OpenCV offers multiple interpolation functions.\n",
    "# In general it is reccommended:\n",
    "# Shrinking: cv2.INTER_AREA\n",
    "# Zooming: cv2.INTER_CUBIC (slow) or cv2.INTER_LINEAR\n",
    "# Default (general): cv2.INTER_LINEAR\n",
    "\n",
    "# Half the size of the image\n",
    "img_resize = cv2.resize(img, (int(cols/2), int(rows/2)), interpolation = cv2.INTER_AREA)\n",
    "\n",
    "# Set to specified size\n",
    "img_resize = cv2.resize(img, (180, 120),interpolation = cv2.INTER_LINEAR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(img_resize)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Translation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenCV provides a function (`cv2.warpAffine`) to carry out affine transformations. The main requirements for this are the original image you wish to apply the transformation and the transformation matrix $M$. We must therefore first compute the transformation matrix $M$. If you know the shift in ($x$, $y$) that you wish to compute you can let:  \n",
    "\n",
    "$M= \\begin{bmatrix} 1 & 0 & ∆x \\\\ 0 & 1 & ∆y \\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.float32([[1,0,25],[0,1,50]])\n",
    "translation = cv2.warpAffine(img, M, (cols, rows))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(translation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly a rotation matrix can be computed using an $M$ matrix following:\n",
    "\n",
    "$M= \\begin{bmatrix} \\cos(\\theta) & -\\sin(\\theta) \\\\ \\sin(\\theta) & \\cos(\\theta) \\end{bmatrix}$\n",
    "\n",
    "OpenCV provides a tool for creating these matrix with the function `cv2.getRotationMatrix2D`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create rotation matrix for 25 degrees rotation in respect to the image center\n",
    "M = cv2.getRotationMatrix2D((cols/2, rows/2), 25, 1)\n",
    "\n",
    "# Create rotation matrix for 25 degrees rotation with different roation center\n",
    "M_offset = cv2.getRotationMatrix2D((250, 250), 25, 1)\n",
    "\n",
    "transformation = cv2.warpAffine(img, M, (cols, rows))\n",
    "transformation_offset = cv2.warpAffine(img, M_offset, (cols, rows))\n",
    "\n",
    "print ('M matrix \\n', M, '\\n\\n')\n",
    "print ('M offset matrix \\n', M_offset)\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.subplot(221), plt.imshow(transformation)\n",
    "plt.subplot(222), plt.imshow(transformation_offset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Affine transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In affine transformations all parallel lines will remain parallel after the transformation. OpenCV provides a function `cv2.getAffineTransformation` to determine the transformation matrix. The minimum number of points for an affine transformation is 3. To determine the $M$ matrix we must therefore make 2 3x1 matrix with the first containing the original co-ordinates of the points, the second will contain the location of the points after the transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "originalPts = np.float32([[50, 50], [100, 150], [180, 120]])\n",
    "transformedPts = np.float32([[60, 75],[110, 190], [175, 100]])\n",
    "\n",
    "# Determine M matrix\n",
    "M = cv2.getAffineTransform(originalPts, transformedPts)\n",
    "\n",
    "# Pass image and M matrix in cv2.warpAffine function\n",
    "transformation = cv2.warpAffine(img, M, (cols, rows))\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.subplot(221), plt.imshow(img)\n",
    "plt.subplot(222), plt.imshow(transformation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Perspective transformation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perspective transformation is similar to affine however instead of preserving parallel lines it preserves collinearity and incidence. As perspective offers an extra degree of translation 4 points are required to determine the $M$ matrix. By providing the co-ordinates in the same manner as the affine transformation, the matrix can be computed using the function `cv2.getPerspectiveTransform`. The transormation can then be carried out using `cv2.warpPerspective`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "originalPts = np.float32([[50, 50], [100, 150], [180, 120], [30, 30]])\n",
    "transformedPts = np.float32([[50, 53], [98, 155], [175, 120], [27, 31]])\n",
    "\n",
    "# Determine M matrix\n",
    "M = cv2.getPerspectiveTransform(originalPts, transformedPts)\n",
    "\n",
    "# Pass image and M matrix in cv2.warpAffine function\n",
    "transformation = cv2.warpPerspective(img, M, (cols, rows))\n",
    "\n",
    "plt.figure(figsize=(7,7))\n",
    "plt.subplot(221), plt.imshow(img)\n",
    "plt.subplot(222), plt.imshow(transformation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Image Thresholding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic thresholding in OpenCV is achieved through the `cv2.threshold` function. There are multiple options to determine how the threshold can be computed, but first we will create a simple binary threshold. \n",
    "\n",
    "**TASK:** Using the code we created earlier how could this be adapted into a simple binary classifier for a grayscale image with a critical value of 150?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First lets load and display a digital elevation model (DEM)\n",
    "img = cv2.imread(image_dir+'lidar_small.tif', 0)\n",
    "plt.imshow(img, 'gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try task here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Extra**: How could we use numpy indexing to remove the need for a 'for' loop?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try task here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now carry this out using OpenCV. The function cv2.THRESH_BINARY defines the threshold style. The availble styles are:\n",
    "\n",
    "`cv2.THRESH_BINARY`\n",
    "\n",
    "`cv2.THRESH_BINARY_INV`\n",
    "\n",
    "`cv2.THRESH_TRUNC`\n",
    "\n",
    "`cv2.THRESH_TOZERO`\n",
    "\n",
    "`cv2.THRESH_TOZERO_INV`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret, thresh = cv2.threshold(img, 50, 255, cv2.THRESH_BINARY)\n",
    "ret, thresh_inv = cv2.threshold(img, 50, 255, cv2.THRESH_BINARY_INV)\n",
    "\n",
    "plot_images([img, thresh, thresh_inv],\n",
    "            ['Original', 'Binary', 'Binary Inverse'],\n",
    "            (15,15),\n",
    "            (1, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adaptive_thresh = cv2.adaptiveThreshold(img,255, \\\n",
    "                            cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY,15,7)\n",
    "\n",
    "adaptive_thresh_g = cv2.adaptiveThreshold(img,255, \\\n",
    "                                cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY,15,7)\n",
    "\n",
    "plot_images([adaptive_thresh, adaptive_thresh_g],\n",
    "            ['Adaptive Thresh', 'Adaptive Thresh w/ Gaussian'],\n",
    "            (10,10),\n",
    "            (1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Otsu's Binarisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otsu's binirisation is a method to automatically determine the threshold for a bimodal image. This attempts to find the optimum threshold value by converting the image to a histogram and finding the value in between the two main peaks. This method is therefore only suitable for bimodal images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ret, otsu = cv2.threshold(img,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "plot_images([img, thresh, adaptive_thresh, otsu],\n",
    "            ['Original', 'Binary Threshold', 'Adaptive Threshold', 'Otsu Binarisation'],\n",
    "            (20, 20),\n",
    "            (1, 4))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convolution Filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/conv2d_padding.gif\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convolutional filters processes images by passing a matrix (kernel) over each pixel within an image. The dot product of the weighted kernel matrix and the overlaying pixels of the image is used to create a new pixel value at the given point. Each pass is often referred to as a convolution.  This can be used for common image processing such as blurring, sharpening and edge detection. OpenCV provides many pre-built functions for kernels convolutions as well as the `cv2.filter2D()` function which can be used to pass a custom-made filter over an image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bluring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bluring (or smoothing) is typically used in image processing for noise reduction and smoothing. Here we achieve this by convolving a low-pass filter kerel over an image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Basic Blurring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = cv2.imread(image_dir+'noise.jpg', 1)\n",
    "img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)\n",
    "img_gs = cv2.imread(image_dir+'noise.jpg', 0)\n",
    "\n",
    "# A simple averaging filter could be made by passing first creating a kerel\n",
    "kernel = np.ones((5,5), np.float32)/25\n",
    "# Using the filter2D function the image can be averaged\n",
    "averaged_img = cv2.filter2D(img, -1, kernel)\n",
    "\n",
    "#Alternatively we can use a pre-made OpenCV blur\n",
    "mean_blur = cv2.blur(img_gs,(5,5))\n",
    "\n",
    "#Median blurs are good for removing salt and pepper effects\n",
    "\n",
    "img_sp = cv2.imread(image_dir+'salt.jpg', 0)\n",
    "\n",
    "median_blur = cv2.medianBlur(img_sp, 5)\n",
    "\n",
    "plot_images([img, mean_blur, img_sp, median_blur],\n",
    "            ['Cat Image', 'Mean Blur', 'Salt and Pepper Image', 'Median Blur'],\n",
    "            (10, 10),\n",
    "            (2, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Gaussian Blurring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gaussian blur uses the gaussian function to calculate each pixels transformation. The equation used to define the Guassian function in 2-dimensions ($x$ and $y$) is defined as: $G_{(x)} = \\frac{1}{\\sqrt{2\\Pi\\sigma^2}}e^{-\\frac{x^2}{2\\sigma^2}}$. \n",
    "\n",
    "The resulting 3x3 kernel used is: $G~(x,y)= \\begin{bmatrix} 1 & 2 & 1 \\\\ 2 & 4 & 2 \\\\ 1 & 2 & 1 \\end{bmatrix} * \\frac{1}{16}$\n",
    "\n",
    "As blurring filters smooth areas with a high variance of intensity they have a strong response to edges. This results in the smoothing of an edge. If our image processing problem is concerned with the edges this is not ideal. Bilateral filters preserve edges by considering the neighbouring pixels and weights assainged to them. The weights have two components, the first is the same weighting as used in the Gaussian kernel above. The second takes into account the difference in intensity between the neighbouring pixels and the evaluated one. A more detailed description of this can be [found here](http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/MANDUCHI1/Bilateral_Filtering.html).\n",
    "\n",
    "If we are primarilly looking at denoising, OpenCV offers more specialised features for such tasks. Take a look at `cv2.fastNlMeansDenoising()` and `cv2.fastNlMeansDenoisingColored()` if you are dealing with images that contain a lot of noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gaus_blur = cv2.GaussianBlur(img, (7, 7), 0)\n",
    "\n",
    "plot_images([img, gaus_blur],\n",
    "            ['Original', 'Gaussian Blurred'],\n",
    "            (15, 15),\n",
    "            (1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Edge Detection and Image Gradients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Edge detection is one of the most fundamental operations used in image processing. Edge detection can be used for feature extraction/segmentation as well as limiting operations to pixels with geometric interest. Here we will look at some of the most common edge detectors discussed in this module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/edge_detection2.png\" >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sobel Operator "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sobel operator is a kernal convolution that is used to find edges within an image. The sobel operator is a discrete differentiation operator which computes an approximation to the derivative of an image that is serperate in the x and y direction. This therefore responds strongly to areas where there is a high variation of pixel intensities. By default this is achieved by passing the following 3x3 kernels over a given image:\n",
    "\n",
    "$G~x= \\begin{bmatrix} -1 & 0 & 1 \\\\ -2 & 0 & 2 \\\\ -1 & 0 & 1 \\end{bmatrix}$ &nbsp;&nbsp;&nbsp; $G~y= \\begin{bmatrix} -1 & -2 & -1 \\\\ 0 & 0 & 0 \\\\ 1 & 2 & 1 \\end{bmatrix}$\n",
    "\n",
    "Using the kernel size argument '`ksize`' the size of the kernel can be increased. If `ksize=-1` the Sharr derivative is used instead of the sobel. When using a 3x3 kernel size OpenCV reccommends using the Sharr kernal which is as follows:\n",
    "\n",
    "$G~x= \\begin{bmatrix} -3 & 0 & 3 \\\\ -10 & 0 & 10 \\\\ -3 & 0 & 3 \\end{bmatrix}$ &nbsp;&nbsp;&nbsp; $G~y= \\begin{bmatrix} -3 & -10 & -3 \\\\ 0 & 0 & 0 \\\\ 3 & 10 & 3 \\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = image_dir+'lidar_small.tif'\n",
    "img = cv2.imread(file, 0)\n",
    "\n",
    "sobelx64f = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=3)\n",
    "abs_sobel64f = np.absolute(sobelx64f)\n",
    "sobelx_8u = np.uint8(abs_sobel64f)\n",
    "\n",
    "sobely64f = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=3)\n",
    "abs_sobel64f = np.absolute(sobely64f)\n",
    "sobely_8u = np.uint8(abs_sobel64f)\n",
    "\n",
    "plot_images([img, sobelx_8u, sobely_8u],\n",
    "            ['Original', 'Sobel X', 'Sobel Y'],\n",
    "            (20, 20),\n",
    "            (1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Whilst the sobel is useful for finding edges in the $x$ and $y$ directions independently, rarely would we use these as our final output. It would be more common to combine the two images which contains all edges. The below function shows us how given an image as input we can return an edge magnitude image for both $x$ any $y$ directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getGradientMagnitude(im):\n",
    "    ddepth = cv2.CV_32F\n",
    "    dx = cv2.Sobel(im, ddepth, 1, 0)\n",
    "    dy = cv2.Sobel(im, ddepth, 0, 1)\n",
    "    dxabs = cv2.convertScaleAbs(dx)\n",
    "    dyabs = cv2.convertScaleAbs(dy)\n",
    "    mag = cv2.addWeighted(dxabs, 0.5, dyabs, 0.5, 0)\n",
    "    return mag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenCV offer a function to return the kernel coefficients\n",
    "\n",
    "sobel7x_deriv = cv2.getDerivKernels(1, 0, 5)\n",
    "sobel7y_deriv = cv2.getDerivKernels(0, 1, 5)\n",
    "\n",
    "# To display this we need to compute the outer product of the two vectors\n",
    "\n",
    "print ('Sobel 7x7 kernel x direction:\\n', np.outer(sobel7x_deriv[1], sobel7x_deriv[0]))\n",
    "print ('\\nSobel 7x7 kernel y direction:\\n', np.outer(sobel7y_deriv[1], sobel7y_deriv[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Laplacian Edge Detector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The laplacian edge detector computes the residual of an image using the equation, $\\frac{\\partial^2src}{\\partial x^2} + \\frac{\\partial^2src}{\\partial y^2}$ where the derivatives for $x$ and $y$ are found using sobel in the $x$ and $y$ direction respectively. Due to this relationship when we call `cv2.Laplacian` the sobel operator is run internally. When the kernel size `ksize=1` the following filter is applied.\n",
    "\n",
    "$G= \\begin{bmatrix} 0 & 1 & 0 \\\\ 1 & -4 & 1 \\\\ 0 & 1 & 0 \\end{bmatrix}$\n",
    "\n",
    "Alternatively we could also adapt this to consider diagonals using the following:\n",
    "\n",
    "$G~x= \\begin{bmatrix} 1 & 1 & 1 \\\\ 1 & -8 & 1 \\\\ 1 & 1 & 1 \\end{bmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scaling can be used to scale the returned pixel intensity\n",
    "\n",
    "laplacian = cv2.Laplacian(img, cv2.CV_8U, ksize=1, scale=1) \n",
    "abs_lap64f = np.absolute(laplacian)\n",
    "lap_8u = np.uint8(abs_lap64f)\n",
    "\n",
    "laplacian_s5 = cv2.Laplacian(img, cv2.CV_8U, ksize=1, scale=5)\n",
    "abs_lap64f_s5 = np.absolute(laplacian_s5)\n",
    "lap_8u_s5 = np.uint8(abs_lap64f_s5)\n",
    "\n",
    "plot_images([lap_8u, lap_8u_s5],\n",
    "            ['Laplacian', 'Laplacian Scales x5'],\n",
    "            (15, 15),\n",
    "            (1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Canny Edge Detector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Developed by John F. Canny in 1986, the canny edge detection algorithms extracts edges using a 5 stage approach:\n",
    "\n",
    "* Noise removal and smoothing using Gaussian filter\n",
    "* Gradient intensity representations for the image computed (Sobel)\n",
    "* Apply non-maximum suppression to remove false edges\n",
    "* Thresholding by applying an upper and lower boundary on the gradient values\n",
    "* Hysteresis thresholding to track along connected edges\n",
    "\n",
    "In OpenCV all of these steps are wrapped into 1 function `cv2.Canny()`. The function takes 4 main arguments:\n",
    "\n",
    "* Lower threshold for hysteresis thresholding (Double)\n",
    "* Upper threshold for hysteresis thresholding (Double)\n",
    "* Kernel size for sobel operator (Int)\n",
    "* Whether or not to use $L_{2}$ norm for image gradient instead of the default $L_{1}$ (Bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "canny_lowthresh = cv2.Canny(img, 25,75,5,L2gradient=True)\n",
    "canny_highthresh = cv2.Canny(img,100,250,5,L2gradient=True)\n",
    "\n",
    "plot_images([canny_lowthresh, canny_highthresh],\n",
    "            ['Low threshold', 'High threshold'],\n",
    "            (15, 15),\n",
    "            (1, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Morphological Filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last topic we will look at is morphological transformations. As the name suggests we are going to make physical changes to the shape of images. This is usually performed on binary images although it is possible to carry this out on RGB images if appropriate measures are taken. The most common forms of morphological transformations are:\n",
    "\n",
    "* Dilation\n",
    "* Erosion\n",
    "* Opening\n",
    "* Closing\n",
    "\n",
    "However, OpenCV offers support for a larger number of transformations. Also, take a look at [scikit-image](http://scikit-image.org) where you can find many unique and novel transformations. Here we will cover the use of the erosion filter to remove non-terrain objects (i.e. buildings, vegetation, cars etc.) from Digital Surface Models (DSMs) to generate Digital Terrain Models (DTMs). A good introduction to all the filters mentioned above can be found [here](https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_imgproc/py_morphological_ops/py_morphological_ops.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = image_dir+'lidar_small.tif'\n",
    "img = cv2.imread(file, 0)\n",
    "\n",
    "# Define our erosion kernel\n",
    "kernel = np.ones((3,3),np.uint8)\n",
    "\n",
    "# Perform erosion\n",
    "erosion = cv2.erode(img,kernel,iterations = 1)\n",
    "\n",
    "# Smooth resulting DTM with a Gaussian Blur\n",
    "erosion_g = cv2.GaussianBlur(erosion, (3, 3), 1)\n",
    "\n",
    "# Alternatively we could using the dilation operator on the eroded\n",
    "# image to remove any small remaining above ground artifacts. This\n",
    "# is very similar to the 'Opening' operator\n",
    "\n",
    "dilation = cv2.dilate(erosion,kernel,iterations = 1)\n",
    "\n",
    "plot_images([img, erosion_g, dilation], \n",
    "            ['Original', 'Erosion', 'Erosion -> Dilation'], \n",
    "            (10,10), \n",
    "            (1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have a satisfactory DTM we can then use this to remove the terrain from the image returning only the above ground objects. We can achieve this by using the `cv2.subtract()` or `cv2.absdiff()` functions. As the difference will never be perfect (i.e. all - ground = above) we are usually left will a very small value for the ground. To remove this we can either perform a manual threshold, or alternatively, use an otsu binarisation. Otsu is a good model here as the image will likely be bimodal (ground and above ground)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subtract DTM from DSM\n",
    "res = cv2.subtract(img, erosion)\n",
    "\n",
    "# Otsu binarisation\n",
    "ret, otsu = cv2.threshold(res,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
    "\n",
    "plot_images([img, res, otsu], \n",
    "            ['Original', 'Subtracted', 'Otsu'], \n",
    "            (20,20), \n",
    "            (1, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extras..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contours"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Until now we have been carrying out segmentation of images but have not had any real way of handling the segmented pixels. The easiest way to handle an extracted feature is to find the contours of the image. Contours are just a curve that connects continous points along a given boundary. To accurately find contours in OpenCV we need to first convert the image into a binary image. This can be achieved using the techniques discussed above (i.e. thresholding, edge detection and morphological filters). By using these tools we create a white feature of interest on a black background, if this is achieved contours can be easily drawn along these boundaries (preferrably white objects on black background). OpenCV returns countours in a list of coordinates that describe the shape of the contour it represents. If the contour is large the list will be large. To reduce this a contour approximation method can be used which only stores the important points (i.e. direction edges). If the contour is a very complicated shape this list may still be very long. This is very common when the contour line is derived from an edge detector and is very noisy. To correct for this we can apply a contour approximation which fits the contour to a similar (and more simple) shape (see the **Douglas-Peucker** algorithm for more information).\n",
    "\n",
    "Let's take a look at how we can extract the contours for all the above-ground objects we just created, and draw them back onto the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find contours\n",
    "image, contours, hierarchy = cv2.findContours(otsu, cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)\n",
    "\n",
    "# Draw contours onto image\n",
    "original_img = img.copy()\n",
    "img = cv2.cvtColor(img, cv2.COLOR_GRAY2RGB)\n",
    "img = cv2.drawContours(img, contours, -1, (0, 0, 255), 3)\n",
    "\n",
    "plot_images([original_img, img], ['Original', 'Contour Image'], (20,20), (1,2))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}