{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local Structures\n",
    "\n",
    "**Students:** Sam van Kampen, Robin Wacanno\n",
    "\n",
    "**Student numbers:** 11874716, 11741163"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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": {},
   "source": [
    "### 7.3\n",
    "#### 1\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": {},
   "source": [
    "#### 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "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": {},
   "source": [
    "Here, X \n",
    "# TODO!!\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": {},
   "source": [
    "#### 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "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": {},
   "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": {},
   "source": [
    "### 7.4\n",
    "#### 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gauss(s):\n",
    "    x, y = np.linspace(-3*s, 3*s, s*2+1), np.linspace(-3*s, 3*s, s*2+1)\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": {},
   "source": [
    "#### 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32052 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), np.linspace(-3*s, 3*s, s*2+1)\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), cmap=plt.cm.YlGnBu_r)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13748 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": {},
   "source": [
    "The plot shown above illustrates the quadratic complexity increase as a function of scale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Gauss1(s):\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3fc9ce5860>]"
      ]
     },
     "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": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f3fc9d30390>"
      ]
     },
     "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",
    "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": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3fc9b3c940>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 13656 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": {},
   "source": [
    "As can be seen from the plot above, the computational complexity is now linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.6\n",
    "#### 1\n",
    "Gegeven is:\n",
    "\\begin{align}\n",
    "\\frac{\\partial}{\\partial x\\partial y} (f\\circledast G^s) & = f\\circledast \\frac{\\partial}{\\partial x\\partial y} G^s\\\\\n",
    " & = f\\circledast \\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": {},
   "source": [
    "#### 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "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": {},
   "source": [
    "#### 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc3c93400>"
      ]
     },
     "execution_count": 43,
     "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": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc3c66dd8>"
      ]
     },
     "execution_count": 44,
     "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": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc3c13e80>"
      ]
     },
     "execution_count": 45,
     "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": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc3b9c710>"
      ]
     },
     "execution_count": 46,
     "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": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc3b78320>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 85860 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.clf()\n",
    "plt.title(\"Fyy\")\n",
    "plt.imshow(gD(f, s, 2, 0), cmap=plt.cm.gray)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.8\n",
    "#### 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "def grad_norm(F, s, x, y):\n",
    "    return np.linalg.norm([gD(F, s, 0, 1), gD(F, s, 1, 0)])\n",
    "\n",
    "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": {},
   "source": [
    "#### 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f3fc36a03c8>"
      ]
     },
     "execution_count": 80,
     "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": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array([[169, 168, 174, ..., 194, 191, 193],\n",
       "       [174, 175, 176, ..., 188, 186, 186],\n",
       "       [171, 178, 172, ..., 195, 188, 191],\n",
       "       ...,\n",
       "       [178, 183, 183, ..., 144, 186, 174],\n",
       "       [187, 180, 157, ..., 148, 183, 152],\n",
       "       [169, 170, 160, ..., 160, 198, 165]], dtype=uint8)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
