{
 "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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "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": 101,
   "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": 103,
   "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": 104,
   "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": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage import convolve\n",
    "\n",
    "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": [
    "## 7.5 Separable Gaussian Convolution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "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": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f715efdd080>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15392 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": 128,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f714a8e35f8>"
      ]
     },
     "execution_count": 128,
     "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": 53,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f714aa439b0>]"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 12832 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 = 10\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": "markdown",
   "metadata": {},
   "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": 94,
   "metadata": {},
   "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": 105,
   "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": [
    "### 7.7 Comparison of theory and practice\n",
    "#### 1 Compare the analytical first-order derivatives in $x$ and $y$ of $f$ with your discrete implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f714a14ad30>"
      ]
     },
     "execution_count": 124,
     "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": 121,
   "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": 121,
     "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": "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
}
