{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "35KpYPTtCWQd"
   },
   "source": [
    "# Local Structure\n",
    "\n",
    "**Students:** Sam van Kampen, Robin Wacanno\n",
    "\n",
    "**Student numbers:** 11874716, 11741163"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "F5HhYy_lCWQg"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ajLf0x63CWQq"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.ndimage import convolve\n",
    "from skimage import data\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "vkQiD4NBCWQx"
   },
   "source": [
    "### 7.3 Analytical Local Structure\n",
    "#### 1 Calculate analytically the following derivatives of this function: fx, fy, fxx, fxy and fyy and give the results in your report. Note that fx denotes the first order partial derivative in x, etc.\n",
    "Let $$f(x, y) = A\\sin(Vx)+B\\cos(Wy)$$\n",
    "so:\n",
    "$$f_x = VA\\cos(Vx)$$\n",
    "$$f_y = -WB\\sin(Wy)$$\n",
    "$$f_{xx} = -V^2A\\sin(Vx)$$\n",
    "$$f_{xy} = 0$$\n",
    "$$f_{yy} = -W^2B\\cos(Wy)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "n5mzdNkKCWQy"
   },
   "source": [
    "#### 2 Display X and Y as images to see what is done. Read the documentation of imshow to understand the extent argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 216
    },
    "colab_type": "code",
    "id": "a0c8m7DSCWQ0",
    "outputId": "4d4d2646-0109-49c1-ef5c-695d8e7d38bc"
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 5536 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.arange(-100,101); y = np.arange(-100,101)\n",
    "\n",
    "X, Y = np.meshgrid(x,y)\n",
    "\n",
    "plt.clf()\n",
    "\n",
    "ax = plt.subplot(121)\n",
    "ax.set_title('X')\n",
    "ax.imshow(X, cmap=plt.cm.gray, extent=(-100,100,100,-100))\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "ax.set_title('Y')\n",
    "ax.imshow(Y, cmap=plt.cm.gray, extent=(-100,100,100,-100))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "hhqsv6RnCWQ8"
   },
   "source": [
    "Here, `X` and `Y` are coordinate matrices based on the coordinate vectors `x` and `y`. \n",
    " \n",
    "The `extend` parameter in the `imshow` function defines the bounding box which will be filled by the drawn image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A0Jo24n4CWQ-"
   },
   "source": [
    "#### 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 216
    },
    "colab_type": "code",
    "id": "qyi515l3CWRA",
    "outputId": "d7885ae5-e42c-48db-f04d-531445a4fe73"
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 6512 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 1; B = 2; V = 6*np.pi/201; W = 4*np.pi/201\n",
    "F = A*np.sin(V*X) + B*np.cos(W*Y)\n",
    "\n",
    "Fx = V*A*np.cos(V*X)\n",
    "Fy = -W*B*np.sin(W*Y)\n",
    "\n",
    "plt.clf();\n",
    "\n",
    "ax = plt.subplot(121)\n",
    "ax.set_title('Fx')\n",
    "ax.imshow(Fx, cmap=plt.cm.gray, extent=(-100,100,100,-100))\n",
    "\n",
    "ax = plt.subplot(122)\n",
    "ax.set_title('Fy')\n",
    "ax.imshow(Fy, cmap=plt.cm.gray, extent=(-100,100,100,-100))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 269
    },
    "colab_type": "code",
    "id": "0Zw_Xq6HCWRH",
    "outputId": "e763373e-23b4-4332-f236-656d4e3c797e"
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 63488 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xx = np.arange(-100, 101, 10)\n",
    "yy = np.arange(-100, 101, 10)\n",
    "\n",
    "XX, YY = np.meshgrid(xx, yy)\n",
    "\n",
    "FFx = Fx[::10, ::10]\n",
    "FFy = Fy[::10, ::10]\n",
    "\n",
    "plt.clf()\n",
    "\n",
    "plt.imshow(F, cmap=plt.cm.gray, extent=(-100,100,100,-100), origin='upper')\n",
    "\n",
    "plt.quiver( xx, yy, FFx, FFy, color='red', angles='xy' )\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "B5qkpfiQDWfl"
   },
   "source": [
    "The quiver plot above clearly shows arrows pointing towards the larthest increase of luminance values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ARbYdZHeCWRN"
   },
   "source": [
    "### 7.4 Gaussian Convolution\n",
    "#### 1 Write a function Gauss(s) that returns a 2D sampled Gaussian function with scale s in a numpy array. You should make sure that an appropriate choice is made for the size of the sampling grid. That choice depends on the scale. Furthermore you should make sure that the sum of all kernel values is equal to one. The origin should be in the middle of the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "PFxeQPLpCWRP"
   },
   "outputs": [],
   "source": [
    "def Gauss(s, n_samp=0):\n",
    "    x, y = np.linspace(-3*s, 3*s, s*2+1+n_samp), np.linspace(-3*s, 3*s, s*2+1+n_samp)\n",
    "    X, Y = np.meshgrid(x,y)\n",
    "    k = (1/(2*np.pi*s**2)*np.exp(-(X**2+Y**2)/(2*s**2)))\n",
    "    sum_k = np.sum(k)\n",
    "    norm_k = k / sum_k\n",
    "    return norm_k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Vsg5-Kh1DmXP"
   },
   "source": [
    "The amount of samples in the function above is calculated using the expression `s*2+1`. This rather small amount of samples is used since a faster increase would cause `MemoryErrors` in further exercises. The domain used to sample the Gaussian is $[-3s, 3s]$ in for both the $x$ and $y$ axes. This truncation point is chosen since 99% of the probability density is inside this domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "qsgx6iFCCWRU"
   },
   "source": [
    "#### 2 In your report you should present a plot of the kernel. In case you have a modern version of matplotlib you might have 3D plotting capabilities (look in the documentation on the web). Else display the array of kernel values for a small value of s (say s=1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 266
    },
    "colab_type": "code",
    "id": "XriiteZkCWRV",
    "outputId": "c13bd540-e8d8-43c8-fd68-c0d9ebadff99"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 68032 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = 1\n",
    "\n",
    "x, y = np.linspace(-3*s, 3*s, s*2+1 + 50), np.linspace(-3*s, 3*s, s*2+1 + 50)\n",
    "X, Y = np.meshgrid(x,y)\n",
    "\n",
    "plt.clf()\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(X, Y, Gauss(1, 50), cmap=plt.cm.YlGnBu_r)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Esvqp4VnFW2x"
   },
   "source": [
    "The plot above shows the the expected graph based on an interpolated Gaussian function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "xcZA3aIOCWRc"
   },
   "source": [
    "#### 3 Given the Gaussian kernel you should be able to do the Gaussian convolution using `convolve(F,Gauss(s),mode='nearest')`. Time the runtime of this function for several values of s (say s=1,2,3,5,7,9,11,15,19) and plot the time versus the scale. Show the plot in your report. What is the order of computational complexity in terms of the scale?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 393
    },
    "colab_type": "code",
    "id": "GEg-uULSCWRd",
    "outputId": "6b176a63-479d-4bb0-e7a2-391e77f6031a"
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13772 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ss = [1, 2, 3, 5, 7, 9, 11, 15, 19]\n",
    "ts = []\n",
    "n = 10\n",
    "f = data.imread(\"data/lena.png\", as_gray=True)\n",
    "\n",
    "for s in ss:\n",
    "    t = time.time()\n",
    "    for _ in range(n):\n",
    "        convolve(f, Gauss(s), mode='nearest')\n",
    "    ts.append((time.time() - t)/n)\n",
    "\n",
    "plt.clf()\n",
    "\n",
    "plt.plot(ss, ts)\n",
    "plt.xlabel(\"s\")\n",
    "plt.ylabel(\"runtime (s)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "juoGT9V0CWRk"
   },
   "source": [
    "The plot shown above illustrates the quadratic complexity increase as a function of scale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "VGUc_ZYFCWRm"
   },
   "source": [
    "### 7.5 Separable Gaussian Convolution\n",
    "\n",
    "#### 1 Write the Gauss1 function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-DiUBSVeCWRn"
   },
   "outputs": [],
   "source": [
    "def Gauss1(s):\n",
    "    \"\"\"\n",
    "    1D sampled Gaussian function. Samples in 2s+1 points.\n",
    "    Returns a list of x values and GS(x) values.\n",
    "    \"\"\"\n",
    "    x = np.linspace(-3*s, 3*s, 2*s+1)\n",
    "    sampled_gauss = np.exp(-pow(x/s, 2)/2)/(s*np.sqrt(2*np.pi))\n",
    "    # Even normaliseren...\n",
    "    sampled_gauss /= np.sum(sampled_gauss)\n",
    "    return x, sampled_gauss\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll see if the function works by plotting a curve for s = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "G5HO1PwaCWRs"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f531490da20>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15640 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(*Gauss1(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, in 2D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "dn8EliX1CWRx"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f531486bef0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 61316 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "s = 10\n",
    "\n",
    "xs, sampled_gauss = Gauss1(s)\n",
    "ys = xs\n",
    "X, Y = np.meshgrid(xs, ys)\n",
    "\n",
    "# To get the values for G(x, y), we can use vector-vector multiplication.\n",
    "g2d = sampled_gauss.reshape((-1, 1)) @ sampled_gauss.reshape((1, -1))\n",
    "ax.plot_surface(X, Y, g2d, cmap=plt.cm.YlGnBu_r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2 Calculate the Gaussian convolution using the separability property and time the 2D Gaussian convolution as a function of the scale s. What is the order of computational complexity now?\n",
    "\n",
    "We'll measure this the same way we measured it for the 2-D Gaussian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2DJAcpSWCWR4"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5314830ac8>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13684 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from skimage import data\n",
    "import time\n",
    "from scipy.ndimage import convolve1d\n",
    "ss = [1, 2, 3, 5, 7, 9, 11, 15, 19]\n",
    "ts = []\n",
    "f = data.imread(\"data/lena.png\", as_gray=True)\n",
    "\n",
    "ntests = 100\n",
    "\n",
    "for s in ss:\n",
    "    t = time.time()\n",
    "    for _ in range(ntests):\n",
    "        g = Gauss1(s)[1]\n",
    "        convolve1d(convolve1d(f, g, mode='nearest', axis=0), g, mode='nearest', axis=1)\n",
    "    ts.append((time.time() - t)/ntests)\n",
    "    \n",
    "plt.plot(ss, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "03xZjPsdCWR8"
   },
   "source": [
    "As can be seen from the plot above, the computational complexity is now linear."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8iebg2oqCWR-"
   },
   "source": [
    "#### 3 (bonus) In Matlab the function convolve automagically checks whether a 2D weight matrix W is separable. Use your Gauss function from a previous exercise to generate a Gaussian kernel. Calculate the SVD of W. What are the singular values? Then what are the 1D convolution kernels that make the 2D convolution separable?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "2XQhJ3BYCWSC"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.08489383] [ True]\n"
     ]
    }
   ],
   "source": [
    "from numpy.linalg import svd\n",
    "\n",
    "W = Gauss(10)\n",
    "u, s, vt = svd(W)\n",
    "nonzero_singular_values = s[abs(s) > 1e-10]\n",
    "print(nonzero_singular_values, nonzero_singular_values == s[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that only the first singular value is non-zero (well, the others are very, very close to zero). This indicates separability, because convolution with a matrix X is the same as convolving with $USV^T$, and since $S[i,j]$ is 0 for $(i, j) \\ne (0, 0)$, `U[:,0] @ S[0,0] @ V^T[:,0]`.\n",
    "\n",
    "The 1-D kernels that make up the 2-D convolution are, then, the first column of $U$ and $V^T$, scaled by $S[0, 0]$ (e.g. by multiplying both by $\\sqrt{S[0,0]}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "HPn_USkxCWSJ"
   },
   "source": [
    "### 7.6 Gaussian Derivatives\n",
    "#### 1 Show analytically (that is using math) that all derivatives of the 2D Gaussian function are separable functions as well.\n",
    "Gegeven is:\n",
    "\\begin{align}\n",
    "\\frac{\\partial}{\\partial x\\partial y} (f * G^s) & = f * \\frac{\\partial}{\\partial x\\partial y} G^s\\\\\n",
    " & = f * \\frac{\\partial}{\\partial x} G^s(x) \\cdot \\frac{\\partial}{\\partial y} G^s(y)\n",
    "\\end{align}\n",
    "Dit is separable want:\n",
    "$$y \\notin var(G^s(x)) \\text{ en } x \\notin var(G^s(y))$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "u4c_fxfUCWSK"
   },
   "source": [
    "#### 2 Write a function `gD(F, s, iorder, jorder)` that returns the Gaussian derivative convolution of image F. The order of differentation is iorder for the first direction (top to bottom!) and jorder for the second direction (left to right). Only differentation up to order 2 is needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "oe4BUegECWSL"
   },
   "outputs": [],
   "source": [
    "def gD(F, s, iorder, jorder):\n",
    "    x, y = np.linspace(-3*s, 3*s, s*2+1), np.linspace(-3*s, 3*s, s*2+1)\n",
    "    if iorder == 0:\n",
    "        yd = 1\n",
    "    elif iorder == 1:\n",
    "        yd = -(y/s**2)\n",
    "    else:\n",
    "        yd = (y**2-s**2)/s**4\n",
    "    if jorder == 0:\n",
    "        xd = 1\n",
    "    elif jorder == 1:\n",
    "        xd = -(x/s**2)\n",
    "    else:\n",
    "        xd = (x**2-s**2)/s**4\n",
    "    kx = xd * Gauss1(s)[1]\n",
    "    ky = yd * Gauss1(s)[1]\n",
    "    return convolve1d(convolve1d(F, ky, mode='nearest', axis=0), kx, mode='nearest', axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BV47EqutCWSQ"
   },
   "source": [
    "#### 3 Visualize the Gaussian derivative of the cameraman image up to order 2 in your report (i.e. the derivatives f (zero order), fx and fy (first order) and fxx, fxy and fyy (second order).\n",
    "Using the gD function you can calculate the discrete first order derivative in x direction with dFx = gD(F, s, 0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-0fxYS3zCWSR"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f5314722080>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 17896 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = data.imread(\"data/cameraman.png\", as_gray=True)\n",
    "s = 1\n",
    "\n",
    "plt.clf()\n",
    "plt.title(\"Fx\")\n",
    "plt.imshow(gD(f, s, 0, 1), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "9UUQ3NuHCWSV"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f53146f7d30>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 21412 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.title(\"Fy\")\n",
    "plt.imshow(gD(f, s, 1, 0), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "ZJPRltZxCWSb"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f53146d2c18>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 85856 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.title(\"Fxx\")\n",
    "plt.imshow(gD(f, s, 0, 2), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "scsCR_MRCWSf"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f5314623eb8>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13280 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.title(\"Fxy\")\n",
    "plt.imshow(gD(f, s, 1, 1), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "A6Npx7tnKobB"
   },
   "source": [
    "### 7.7 Comparison of theory and practice\n",
    "\n",
    "#### 1 Compare the analytical first-order derivatives in $x$ and $y$ of $f$ with your discrete implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "uqissvcYKxvW"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f53145adbe0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 6432 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute the discrete derivatives (with Gaussian smoothing)\n",
    "dFx = gD(F, 1, 0, 1)\n",
    "dFy = gD(F, 1, 1, 0)\n",
    "\n",
    "# Plot the x-derivative\n",
    "p1 = plt.subplot(2, 2, 1)\n",
    "p1.imshow(dFx, cmap=plt.cm.gray)\n",
    "p2 = plt.subplot(2, 2, 2)\n",
    "p2.imshow(Fx, cmap=plt.cm.gray)\n",
    "\n",
    "# Plot the y-derivative\n",
    "p3 = plt.subplot(2, 2, 3)\n",
    "p3.imshow(dFy, cmap=plt.cm.gray)\n",
    "p4 = plt.subplot(2, 2, 4)\n",
    "p4.imshow(Fy, cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, the images are quite similar (discrete numerical derivatives on LHS, analytical on RHS)\n",
    "\n",
    "#### 2 Calculate `Fx/dFx` and convince yourself that it is nearly constant in every pixel. Where does it deviate from the constant value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541],\n",
       "       [30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541],\n",
       "       [30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541],\n",
       "       ...,\n",
       "       [30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541],\n",
       "       [30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541],\n",
       "       [30.78516541, 15.3586906 , 15.3586906 , ..., 15.3586906 ,\n",
       "        15.3586906 , 30.78516541]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fx/dFx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, it is constant in all values except for at the left and right edges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "S83sT6h6CWSn"
   },
   "source": [
    "### 7.8 Canny Edge Detector\n",
    "#### 1 Implement the Canny edge detector using the gD function that you have implemented in previous sections.\n",
    "\n",
    "The result of the function E = canny(F, s) should be an image E where the value in a pixel equals the gradient norm (in that pixel) in case it is an edge pixel and zero otherwise.\n",
    "\n",
    "To find the zero crossings in an image you can write a simple function that looks in each 3x3 neighborhood of a pixel whether there are both negative and positive values on opposite sides of the central pixel. This is not a perfect solution but it does the job quite well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "_9gTS48qCWSp"
   },
   "outputs": [],
   "source": [
    "def canny(F, s, cutoff=10):\n",
    "    F_c = F.copy()\n",
    "    \n",
    "    Fx = gD(F, s, 0, 1)\n",
    "    Fy = gD(F, s, 1, 0)\n",
    "    Fxx = gD(F, s, 0, 2)\n",
    "    Fyy = gD(F, s, 2, 0)\n",
    "    Fxy = gD(F, s, 1, 1)\n",
    "    \n",
    "    ylen, xlen = F.shape\n",
    "    \n",
    "    for y, row in enumerate(F):\n",
    "        for x, pix in enumerate(row):\n",
    "            first_greater = np.sqrt(Fx[x,y]**2 + Fy[x,y]**2) > cutoff\n",
    "            \n",
    "            contains_negative = False\n",
    "            contains_positive = False\n",
    "            \n",
    "            for xa in range(max(0, x - 1), min(xlen, x + 1)):\n",
    "                for ya in range(max(0, y - 1), min(ylen, y + 1)):\n",
    "                    second_val = (Fx[xa, ya]**2)*Fxx[xa, ya] + 2*Fx[xa,ya]*Fy[xa,ya]*Fxy[xa,ya] + (Fy[xa,ya]**2)*Fyy[xa,ya]\n",
    "                    if (second_val <= 0):\n",
    "                        contains_negative = True\n",
    "                    elif (second_val >= 0):\n",
    "                        contains_positive = True\n",
    "            \n",
    "            if (first_greater and contains_negative and contains_positive):\n",
    "                F_c[x, y] = 255\n",
    "            else:\n",
    "                F_c[x, y] = 0\n",
    "            \n",
    "    \n",
    "    return F_c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uYChIgPNCWSv"
   },
   "source": [
    "#### 2 Test your canny function on the camereman image. Show the results in your report for several scales s used in the derivatives calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "8EKxegzRCWSv"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f53148de470>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 10576 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(canny(f, 1), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "3BKe_ZeUCWSy"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f5316273320>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13784 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(canny(f, 2), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f5316293b70>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 11180 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(canny(f, 5), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "Local Structure.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
