VecCalc ImProcessing (FoAM).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vector Calculus and image processing\n",
    "\n",
    "**Objectives:** In this notebook you will\n",
    "- see differential operators acting on images\n",
    "- learn about a simple partial differential equation\n",
    "- discretize and numerically approximate the solutions of this PDE\n",
    "- make a connection to energy minimization and Gradient descent\n",
    "- apply this to image denoising and reconstruction\n",
    "\n",
    "The following box you have to run once at the beginning of each session in order to load all packages that are needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ArgumentError",
     "evalue": "ArgumentError: Package Images not found in current path:\n- Run `import Pkg; Pkg.add(\"Images\")` to install the Images package.\n",
     "output_type": "error",
     "traceback": [
      "ArgumentError: Package Images not found in current path:\n- Run `import Pkg; Pkg.add(\"Images\")` to install the Images package.\n",
      "",
      "Stacktrace:",
      " [1] require(::Module, ::Symbol) at .\\loading.jl:823",
      " [2] top-level scope at In[1]:1"
     ]
    }
   ],
   "source": [
    "using Images, Gadfly, TestImages, ImageView, Colors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Images as scalar fields\n",
    "\n",
    "Gray scale images are scalar fields $f\\colon \\mathbb{R}^2\\to \\mathbb{R}$. A number, usually between 0 (totally black) and 255 (totally white) (8bits), encodes the gray level at a certain position. Thus one may also think of images as landscapes. A digital image is just a matrix. More precisely, a rectangle of certain size is discretized and given via an $M\\times N$ matrix of pixel positions. The entries of the matrix encode the gray levels of the pixels. So a typical digital gray scale image is a function $f\\colon\\{1,\\ldots,M\\}\\times\\{1,\\ldots,N\\}\\to\\{0,\\ldots,255\\}$ where $f(m,n)$ is the gray level of the pixel at position $(m,n)$. Let us first read a jpeg image and tranform it into a matrix of gray level entries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "nbpresent": {
     "id": "c5194dd5-2c9d-4998-97b7-c246184f0703"
    }
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: load not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: load not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[2]:1"
     ]
    }
   ],
   "source": [
    "#im = convert(Array{Float64}, testimage(\"cameraman\"))\n",
    "im = convert(Array{Float64}, Gray.(load(\"huggingzebras.jpg\")))\n",
    "M, N = size(im)\n",
    "Gray.(im)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# The Gradient\n",
    "\n",
    "The main question on zebras is: Are they black with white stripes, or are they white with black stripes? But mathematics has nothing to say about this! Instead we do the following:\n",
    "\n",
    "Discrete derivatives are computed as finite difference approximations of derivatives. For instance, as partial derivatives in the coordinate directions one might consider\n",
    "\n",
    "$$ \\delta_1f(n_1,n_2) := f(n_1+1,n_2) - f(n_1,n_2) $$\n",
    "$$ \\delta_2f(n_1,n_2) := f(n_1,n_2+1) - f(n_1,n_2) $$\n",
    "\n",
    "One has to be careful at the boundary of the domain. You could for example consider periodic boundary conditions which corresponds to computing the indices $n_i+1$ modulo $N$, i.e. you just let $n_i+1=1$ for $n_i=N$ and $n_i-1=N$ for $n_i=1$. In our below implementation we won't even do that, but just disregard the boundary and possible artifacts created there.\n",
    "\n",
    "A discrete gradient is\n",
    "\n",
    "$$\\nabla f(n):= \\begin{bmatrix}\\delta_1f(n)\\\\ \\delta_2f(n) \\\\ \\end{bmatrix} \\in \\mathbb{R}^{2}$$\n",
    "\n",
    "i.e. the gradient maps images to vector fileds, $\\nabla\\colon \\mathbb{R}^{N\\times M}\\to \\mathbb{R}^{N\\times M\\times 2}$.\n",
    "\n",
    "Let's compute and look at the partial derivatives in the two coordinate directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: M not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: M not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[3]:1"
     ]
    }
   ],
   "source": [
    "gx = zeros(M,N)\n",
    "gy = zeros(M,N)\n",
    "j = 1:(M-1)\n",
    "k = 1:(N-1)\n",
    "\n",
    "gx[j,k] = im[j,k.+1] - im[j,k]\n",
    "gy[j,k] = im[j.+1,k] - im[j,k]\n",
    "\n",
    "#Display matrix as grayscale image\n",
    "Gray.([gx gy])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Denote the Euclidean length of the gradient at a particular $(n_1,n_2)$ by\n",
    "\n",
    "$$\n",
    "|| \\nabla f(n_1,n_2) || := \\sqrt{(\\delta_1f(n_1,n_2))^2+(\\delta_2f(n_1,n_2))^2}\n",
    "$$\n",
    "\n",
    "This is dispayed next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: gx not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: gx not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[4]:1"
     ]
    }
   ],
   "source": [
    "Gr_mag = sqrt.(gx.^2 + gy.^2)\n",
    "Gray.(Gr_mag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The range of the gradient magnitude is huge. And we only have 256 gray levels. So the small gradient magnitudes are not even visible. We can look at this on a log scale. Then for example the tree in the background becomes visible. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View Gr_mag after applying a logarithm ... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "944dc4b0-5545-4874-bc32-ce6dc16d4b10"
    }
   },
   "source": [
    "# Gradient descent\n",
    "\n",
    "We know that the gradient points into the direction of the steepest ascent and that its length measures how steep the ascent in this direction is. If we are anywhere in a lanscape of a scalar field, we can compute the gradient and walk a step into this direction towards a maximum or a step into the opposite direction towards a minimum.\n",
    "Note that, in general, at every point the direction and magnitude of the gradient changes.\n",
    "\n",
    "\n",
    "**Example:**  Consider the quadratic dwell\n",
    "\n",
    "$$ g(x,y) = \\frac{1}{2}(\\alpha x^2 +y^2 )$$\n",
    "\n",
    "for $x,y\\in\\mathbb{R}$. Of course, we know that the minimum of this function is at $(0,0)$. Let's see what a gradient descent does. Since\n",
    "we have a formula for $g$ we can compute the gradient analytically\n",
    "\n",
    "$$ \\nabla g (x,y) = \\begin{bmatrix}\\alpha x\\\\ y \\\\ \\end{bmatrix} $$\n",
    "\n",
    "The idea is to compute the gradient at the point we are currently at, and to walk a small step into the opposite direction. Then, compute the gradient at the new position and again walk a step into precisely to opposite of this direction. \n",
    "That is, starting from some point $p^{(0)}=(x^{(0)},y^{(0)})^\\top$, iteratively compute a sequence of points $p^{(k)}=(x^{(k)},y^{(k)})^\\top$ as follows\n",
    "\n",
    "$$ p^{(k+1)} = p^{(k)} - s \\nabla g (p^{(k)}) $$\n",
    "\n",
    "\n",
    "But what step size $s>0$ to choose? Essentially the choices fall into three categories: the Good, the Bad, and the Ugly. If the stepsize is too small, here say $s=0.02$, the gradient descent will be very slow. Choosen in the right ballpark, the convergence is faster, say for $s=0.18$. If the step size is chosen too large, then the gradient descent may not converge to anything at all, say for $s=0.22$. It is a question adressed in numerical real analysis how to choose the step size for a particular function at hand. Describe what happens for different values of $s$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: Geom not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: Geom not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[6]:4"
     ]
    }
   ],
   "source": [
    "X = Array(-20:0.1:20); Y = X #discretised domain for the function\n",
    "alpha = 10\n",
    "Z = 0.5*((alpha * (X.^2)) .+ Y'.^2) #compute all values z=g(x,y)\n",
    "p = plot(x=X, y=Y, z=Z, Geom.contour(levels=40))\n",
    "\n",
    "numiters = 30 #number of iterations in gradient descent\n",
    "s = 0.18;     #choose a stepsize\n",
    "\n",
    "# Create empty array to store the iterates\n",
    "X = Array{Float64}(undef, numiters) \n",
    "Y = Array{Float64}(undef, numiters)\n",
    "X[1] = 4; Y[1] = 19.8 #choose a starting point\n",
    "\n",
    "# Gradient descent\n",
    "for k = 2:numiters\n",
    "    Xold, Yold = X[k-1], Y[k-1]\n",
    "    X[k], Y[k] = Xold - s*alpha*Xold, Yold - s*Yold\n",
    "end\n",
    "push!(p, layer(x=X, y=Y, Geom.point, Geom.path))\n",
    "draw(SVG(700px, 500px), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "122607dd-b961-4dfe-aab1-1e1491005bda"
    }
   },
   "source": [
    "The iterates may approach the minimum rather slow if the problem has bad conditions. Say, if $\\alpha = 100$ then the quadratic is much like a skateboard halfpipe. With step size $s=0.02$ we converge very slow because of serious zigzagging. But for bigger $s$ we diverge. There are ways to overcome such issues, for instance by changing the step sizes throughout the iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: Geom not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: Geom not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[7]:4"
     ]
    }
   ],
   "source": [
    "X = Array(-20:0.1:20); Y=X #discretised domain for the function\n",
    "alpha = 100\n",
    "Z = 0.5*((alpha * (X.^2)) .+ Y'.^2) #compute all values z=g(x,y)\n",
    "p = plot(x=X, y=Y, z=Z, Geom.contour(levels=40))\n",
    "\n",
    "numiters = 40 #number of iterations in gradient descent\n",
    "s = 0.02;     #choose a stepsize\n",
    "X = Array{Float64}(undef, numiters)\n",
    "Y = Array{Float64}(undef, numiters)\n",
    "X[1]=4; Y[1] = 19.8 #choose a starting point\n",
    "\n",
    "# Gradient descent\n",
    "for k=2:numiters\n",
    "    Xold, Yold = X[k-1], Y[k-1]\n",
    "    X[k], Y[k] = Xold - s*alpha*Xold, Yold - s*Yold\n",
    "end\n",
    "push!(p, layer(x=X, y=Y, Geom.point, Geom.path))\n",
    "draw(SVG(700px, 500px), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You probably have observed that in a complicating landscape gradient descent will in general not lead us to a global minimum: What if there are many local minima or maxima? One can get trapped in one of them. The choice of starting point can then determine to what point we converge.\n",
    "Such issues cannot happen if the function you are dealing with are \"convex\". Also one has to think about whether a minimum exists at all.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "83ac1886-9598-49d9-88de-fff6898509b4"
    }
   },
   "source": [
    "# Divergence and Laplace operator\n",
    "\n",
    "Above, we discretized the partial derivative by a simple forward difference, e.g. \n",
    "$ \\delta_1f(n_1,n_2) = f(n_1+1,n_2) - f(n_1,n_2) $. Now, is this the approximation of the partial derivative at $(n_1,n_2)$ or at $(n_1+1,n_2)$? Well, it's not clear. Maybe it is an approximation at $(n_1+1/2,n_2)$? But there are no half pixels on the digital image. Moreover, we might as well have discretized the partial derivative with  backwards differences:\n",
    "\n",
    "$$ \\delta_1^*f(n_1,n_2) := f(n_1,n_2) - f(n_1-1,n_2) $$\n",
    "$$ \\delta_2^*f(n_1,n_2) := f(n_1,n_2) - f(n_1,n_2-1) $$\n",
    "\n",
    "To compute a second derivative, we don't want to apply the forward difference twice. It seems more reasonable to do a forward and then a backwards difference, in order not to move to far away from $(n_1,n_2)$.\n",
    "\n",
    "\n",
    " - **Exercise:** Let $f,g\\in\\mathbb{R}^N$ and let $\\delta f$ denote the forward and $\\delta^* f$ the backward difference. Assuming periodic boundary conditions, show that for the Euclidean inner product we have\n",
    "$$\n",
    "\\langle \\delta f, g \\rangle = -\\,\n",
    "\\langle f,\\delta^* g\\rangle\n",
    "$$\n",
    "This can be interpreted as a discrete version of the integration by parts formula\n",
    "$$\n",
    "\\int_0^1 f'g = - \\int_0^1 fg'\n",
    "$$\n",
    "for smooth functions $f,g$ with $f(0)=f(1)$ and $g(0)=g(1)$.\n",
    "\n",
    "\n",
    "\n",
    "Define the discrete divergence $\\text{div}\\colon \\mathbb{R}^{N\\times M\\times 2}\\to \\mathbb{R}^{N\\times M}$ using the backwards differences,\n",
    "\n",
    "$$ \\text{div}(v)(n) := \\delta_1^*v_1(n)+\\delta_2^*v_2(n)$$.\n",
    "\n",
    "Notice that this definition implies a discrete version of Green's theorem:\n",
    "$$\n",
    "\\langle \\nabla f, v  \\rangle_{\\mathbb{R}^{N\\times M\\times 2}}  \n",
    "= - \\, \\langle f, \\text{div}(v)  \\rangle_{\\mathbb{R}^{N\\times M}}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- **Exercise:** Let $x\\in\\mathbb{R}^n$ be a vector and $B\\in\\mathbb{R}^{n\\times n}$ be a matrix. Let $B^\\top$ be the matrix derived from $B$ by letting the $k$-th row of $B$ become the $k$-th column of $B^\\top$. $B^\\top$ is called the transpose of $B$. Show that\n",
    "$$\n",
    "\\langle x,Bx\\rangle = \\langle B^{\\top}x,x\\rangle\n",
    "$$\n",
    "With this in mind it is reasonable to think of the gradient and the negative transpose of the divergence\n",
    "$$\n",
    "\\nabla^{\\top}=-\\text{div}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the Laplace operator is defined as the sum of the second order partial derivatives\n",
    "\n",
    "$$ \\Delta f = \\text{div}(\\nabla f)$$\n",
    "\n",
    "For our discrete images that means\n",
    "\n",
    "$$\n",
    "\\Delta f(n_1,n_2) = f(n_1+1,n_2) + f(n_1,n_2+1) + f(n_1-1,n_2) + f(n_1,n_2-1) - 4f(n_1,n_2) .\n",
    "$$\n",
    "\n",
    " - **Exercise:** This indeed approximates the continuous second derivative. For a smooth $g\\colon\\mathbb{R}\\to\\mathbb{R}$ you can try to show (using second order Taylor or otherwise) that\n",
    "$$\n",
    "\\lim_{h\\to 0 } \\frac{f'(x+h)-2f(x)+f'(x-h)}{h^2} = f''(x).\n",
    "$$\n",
    "Again this shows that our discrete Laplace approximates $\\frac{\\partial^2f}{\\partial x_1^2}+\\frac{\\partial^2f}{\\partial x_2^2}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integral operators\n",
    "\n",
    "\n",
    "Consider the total \"energy\" of the gradient vector field\n",
    "\n",
    "$$J(f):=\\frac{1}{2}\\|\\nabla f\\|^2_2 := \\frac{1}{2}\\sum_n\\|\\nabla f(n)\\|^2.$$\n",
    "\n",
    "The operator $J$ maps a function to a number, here the total energy of the gradient vector field, i.e. it measures the smoothness of $J$. The more uniformly smooth $f$, the smaller $J(f)$ will be. If $f$ was defined on a continuous rather than a discrete domain we would replace the sum by an integral.\n",
    "\n",
    "A related operator is the \"total variation\" defined as the summation of the individual gradient lengths\n",
    "\n",
    "$$J_{TV}(f):= \\|\\nabla f\\|_1 := \\sum_n\\|\\nabla f(n)\\|.$$\n",
    "\n",
    "\n",
    " - **Exercise:** What is the total variation and its geometric interpretation for the following image?\n",
    " \n",
    " $$\n",
    " \\begin{pmatrix}\n",
    " 0&0&0&0&0\\\\\n",
    " 0&1&1&0&0\\\\\n",
    " 0&1&1&1&0\\\\\n",
    " 0&1&1&0&0\\\\\n",
    " 0&0&0&0&0\\\\\n",
    " \\end{pmatrix}\n",
    " $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient flow PDE and Energy minimization\n",
    "\n",
    "We now apply energy minimization to de-noise images. Consider the following image and its noisy version. Data aquired in real life with sensors or cameras is always corrupted by a certain amount of noise. Here we artificially create a noisy image and would like to recover the original image. One idea is to denoise by reducing the gradient energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "nbpresent": {
     "id": "ce446957-d62d-4374-85f5-4f37783ebaac"
    }
   },
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: load not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: load not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[8]:1"
     ]
    }
   ],
   "source": [
    "im = convert(Array{Float64}, Gray.(load(\"chomsky.jpg\")))\n",
    "im = im[51:306, 1:256]\n",
    "#im = convert(Array{Float64}, testimage(\"cameraman\"))\n",
    "M, N = size(im)\n",
    "# Add some noise of variance sigma\n",
    "sigma = 0.06  \n",
    "randmat = randn(size(im));\n",
    "im_noisy = im + sigma*randmat\n",
    "\n",
    "Gray.([im im_noisy])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "b97ebe18-6e07-4bfc-829d-d11d56c55c17"
    }
   },
   "source": [
    "## Gradient flow\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "The gradient of the scalar field $J$ is a vector field. Locally up to first order it describes the variation of $J$, i.e.\n",
    "\n",
    "$$\n",
    "J(f+\\epsilon) = J(f) + \\nabla J(f) \\cdot \\epsilon + o(\\|\\epsilon\\|)\n",
    "$$\n",
    "\n",
    "\n",
    "Gradient descent should bring us towards a minimizer of $J$, i.e. we attempt to iterate\n",
    "\n",
    "$$\n",
    "f^{(k+1)}=f^{(k)}-\\tau\\nabla J(f^{(k)}).\n",
    "$$\n",
    "\n",
    "The latter equation can be seen as a finite difference approximation of the gradient flow partial differential equation (PDE) for a function $f=f(t,x)$ that depends on a continuous time parameter $t\\geq 0$, namely\n",
    "\n",
    "$$\n",
    "\\frac{\\partial f}{\\partial t} = - \\nabla J(f) \\quad\\text{ with initial values }\\quad f(0,\\cdot)=f\n",
    "$$\n",
    "\n",
    "The discretization follows choosing the discrete time steps $t_k=k\\tau.$\n",
    " \n",
    " \n",
    " \n",
    " - **Example:** As a simple case to compare with, here is an ordinary differential equation (ODE): Suppose $g(t)$ is a function of one variable $t\\in\\mathbb{R}$ only, and consider the ordinary differential equation $g'=kg$ with initial value $g(0)=g_0$. This simple equation can be directly solved for $g$ via separation of the variables as follows\n",
    " $$\n",
    " \\frac{dg}{dt}=kg \\quad\\Longrightarrow\\quad \\frac{dg}{g}=kdt\n",
    " $$\n",
    " Integrating yields $\\ln g=kt+c$ and thus $g(t)=Ce^{kt}$ where $g>0$. The constant $k$ is determined by the initial condition $g_0=g(0)=Ce^0=C$. By differentiating you can check that $g(t)=g_0e^{kt}$ indeed solves the differential equation. This was an ODE because $g$ depends only on one variable. A PDE is an equation which gives a relationship between the partial derivatives of a function of several independent variables and the gradient flow above is an example of a PDE.\n",
    " \n",
    " - **Exercise:** The equation of the previous example belongs to a more general class of ODE's, those with *seperate variables*. One can attempt to solve them by directly integrating as follows: If $g(y)\\neq 0$ then\n",
    "$$\n",
    "y'=\\frac{dy}{dx}=f(x)g(y)\\quad\\Longrightarrow\n",
    "\\quad \\frac{dy}{g(y)}=f(x)\\,dx\n",
    "\\quad\\Longrightarrow\n",
    "\\quad \\int\\frac{dy}{g(y)}=\\int f(x)\\,dx\n",
    "$$\n",
    "leads to an implicit formula for $y$. In addition one may have solutions of the form $y=y_0$ namely if for some constant $y_0$ the equation $g(y_0)=0$ holds.\n",
    "Find the solution to $y'=-\\frac{x}{y}$ (where $y\\neq 0$).\n",
    "\n",
    "\n",
    "\n",
    "## The heat flow\n",
    "\n",
    "- **Exercise:** Let $x\\in\\mathbb{R}^n$ be a vector and $B\\in\\mathbb{R}^{n\\times n}$ be a matrix. Show that\n",
    "$$\n",
    "\\nabla \\langle x,Bx\\rangle = B^{\\top}x + Bx.\n",
    "$$\n",
    "(Hint: Show $\\partial_k \\langle x,Bx\\rangle=(B^{\\top}x)_k + (Bx)_k$) Use this to show that\n",
    "$$\n",
    "\\nabla \\frac{1}{2}\\|Bx\\|^2=B^\\top Bx\n",
    "$$\n",
    "(Hint: What is $(B^\\top B)^\\top$)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "For $J$ being the total gradient energy we therefore get \n",
    "\n",
    "\n",
    "$$ \\nabla J(f) = -\\Delta f$$\n",
    "\n",
    "Thus in this case the gradient flow results precisely in a very famous PDE from physics: The heat equation\n",
    "\n",
    "$$\n",
    "\\frac{\\partial f}{\\partial t} = \\Delta f \\quad\\text{ with }\\quad f(0,\\cdot)=f\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[9]:4"
     ]
    }
   ],
   "source": [
    "# Parameters\n",
    "tau = 1/4 \n",
    "niter = 300\n",
    "\n",
    "M, N =size(im);\n",
    "Laplace = zeros(size(im))\n",
    "m = 2:(M-1)\n",
    "n = 2:(N-1)\n",
    "\n",
    "# Gradient descent\n",
    "f = im;     \n",
    "for i = 1:niter\n",
    "    Laplace[n,m] = f[n.-1,m] + f[n.+1,m] + f[n,m.-1] + f[n,m.+1] - 4*f[n,m];\n",
    "    f = f + tau*Laplace;\n",
    "end\n",
    "\n",
    "Gray.([im f])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The total variation flow\n",
    "\n",
    "Observe the effect of the heat flow. Noise would get blurred out and disappear, but so would be all important information in the image. The heat diffusion is isotropic and uniform. It does not regard any edges in the image. We would like to slow down the diffusion across the edges. That way one could hope to remove noise from an image without blurring edges.\n",
    "\n",
    "We could use the magnitude of the gradient to produce anisotropic diffusion that preserves edges. Instead of using\n",
    "$\\Delta f=\\text{div}(\\nabla f)$ how about we only diffuse in the direction of the gradient if the gradient magnitude is small, i.e. not across strong edges. We could try to achieve this by using\n",
    "\n",
    "$$\n",
    "\\text{div}\\left(\\frac{\\nabla f}{\\|\\nabla f\\|}\\right)\n",
    "$$\n",
    "for all points for which $\\|\\nabla f\\|\\neq 0$. As it happens, this is the gradient flow for the total variation. That is, one can show,\n",
    "$$\n",
    "(\\nabla J_{TV}(f))(n) = - \\text{div}\\left(\\frac{\\nabla f}{\\|\\nabla f\\|}\\right)(n)\n",
    "$$\n",
    "whenever $\\nabla f(n)\\neq 0$.\n",
    "\n",
    "There is still the problem that the flow is unstable: Division by zero is not defined and wherever the gradient is very small, $(\\nabla J_{TV}(f))(n)$ blows up. In many natural images, there will be regions where the gradient is zero or very small.\n",
    "To get out of this trouble, consider a smoothed version of the total variation\n",
    "\n",
    "$$J_{TV}^{\\epsilon}(f)=\\sum_n \\sqrt{\\epsilon^2+\\|\\nabla f(n)\\|^2}$$\n",
    "\n",
    "\n",
    "- ** Exercise:** To see what is going on here sketch for $s\\in[-1,1]$ the functions $g(s)=|s|$ and $g_{\\epsilon}(s)=\\sqrt{\\epsilon^2+s^2}$. What happens to $g_{\\epsilon}(s)$ when $\\epsilon\\to 0$ ?\n",
    "\n",
    "\n",
    "We want to use this stabilized TV energy, i.e.\n",
    "$$ \\nabla J_{TV}^{\\epsilon}(f)=-\\text{div}\\left(\\frac{\\nabla f}{\\sqrt{\\epsilon^2+\\|\\nabla f\\|^2}}\\right)$$\n",
    "\n",
    "Choosing a smaller $\\epsilon$ makes the flow closer to a minimization of the total variation, but makes the computation unstable due to divisions by very small numbers. Try different $\\epsilon$ in the following visiulaization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[10]:4"
     ]
    }
   ],
   "source": [
    "epsilon = 0.01;\n",
    "\n",
    "f = im;\n",
    "M, N = size(im)\n",
    "Laplace = zeros(size(im))\n",
    "j = 1:(M-1)\n",
    "k = 1:(N-1)\n",
    "m = 2:(M-1)\n",
    "n = 2:(N-1)\n",
    "\n",
    "fdx = f[j.+1,k] - f[j,k]\n",
    "fdy = f[j,k.+1] - f[j,k]\n",
    "Gnorm = sqrt.(epsilon^2 .+ fdx.^2 .+ fdy.^2);\n",
    "fdx = fdx ./ Gnorm\n",
    "fdy = fdy ./ Gnorm;\n",
    "div = fdx[n,m] - fdx[n.-1,m] + fdy[n,m] - fdy[n,m.-1];\n",
    "    \n",
    "Gray.(div)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The smoothed total variation flow now is\n",
    "\n",
    "$$\n",
    "\\frac{\\partial f}{\\partial t} = \\text{div}\\left(\\frac{\\nabla f}{\\sqrt{\\epsilon^2+\\|\\nabla f\\|^2}}\\right) \n",
    "\\quad\\text{ with } \\quad f(0,\\cdot)=f.\n",
    "$$\n",
    "\n",
    "\n",
    "In practice the flow is again being computed via gradient descent. For the smoothed total variation flow to converge, one needs to impose $\\tau<\\epsilon/4$. Thus being more faithful to the TV energy requires smaller time steps and thus yields a slower algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[11]:5"
     ]
    }
   ],
   "source": [
    "# Parameters:\n",
    "epsilon= 0.006\n",
    "tau = epsilon/5\n",
    "niter = 400\n",
    "\n",
    "M, N = size(im);\n",
    "fdx = zeros(size(im))\n",
    "fdy = zeros(size(im))\n",
    "div = zeros(size(im))\n",
    "\n",
    "j = 1:(M-1); k = 1:(N-1); m = 2:(M-1); n=2:(N-1)\n",
    "\n",
    "# Gradient descent:\n",
    "ftv = im\n",
    "for i=1:niter\n",
    "    fdx[j,k] = ftv[j.+1,k] - ftv[j,k]\n",
    "    fdy[j,k] = ftv[j,k.+1] - ftv[j,k]\n",
    "    Gnorm = sqrt.( epsilon^2 .+ fdx.^2 .+ fdy.^2)\n",
    "    fdx = fdx ./ Gnorm\n",
    "    fdy = fdy ./ Gnorm\n",
    "    div[n,m] = fdx[n,m] - fdx[n.-1,m] + fdy[n,m] - fdy[n,m.-1]\n",
    "    ftv += tau*div;\n",
    "end\n",
    "Gray.([im ftv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try denoising using gradient flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[12]:5"
     ]
    }
   ],
   "source": [
    "# Parameters:\n",
    "epsilon = 0.0015\n",
    "tau = epsilon/5\n",
    "niter = 250\n",
    "\n",
    "M,N = size(im)\n",
    "fdx = zeros(size(im))\n",
    "fdy = zeros(size(im))\n",
    "div = zeros(size(im))\n",
    "j = 1:(M-1); k = 1:(N-1); m = 2:(M-1); n = 2:(N-1)\n",
    "\n",
    "# Gradient descent:\n",
    "ftv = im_noisy\n",
    "for i = 1:niter\n",
    "    fdx[j,k] = ftv[j.+1,k] - ftv[j,k]\n",
    "    fdy[j,k] = ftv[j,k.+1] - ftv[j,k]\n",
    "    Gnorm = sqrt.(epsilon^2 .+ fdx.^2 .+ fdy.^2)\n",
    "    fdx = fdx ./ Gnorm\n",
    "    fdy = fdy ./ Gnorm\n",
    "    div[n,m] = fdx[n,m] - fdx[n.-1,m] + fdy[n,m] - fdy[n,m.-1]\n",
    "    ftv += tau*div\n",
    "end\n",
    "\n",
    "Gray.([im_noisy ftv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDE flow for inpainting\n",
    "\n",
    "\n",
    "Data aquired from measurement devices like digital cameras or ccd's in microscopes or telescopes will always be corrupted, for example by random noise. Another possible corruption could be missing data, i.e. one only has scattered data of a function and needs to interpolate the function for the points where one didn't aquire measurements. Of course, for a totally arbitrary function the values of the function at those points could have been virtually anything and there is no hope of finding the true but unknown values. However, if we are dealing with natural images, we may reasonably assume they have small overall gradient energy or small total variation. This holds in particular for example for many medical images. \n",
    "\n",
    "Again, let us artificially create such a situation by randomly deleting a percentage of the pixels and suppose that this actually was the data we have to work with. How can we reconstruct the original from this corrupted image?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: M not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: M not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[13]:3"
     ]
    }
   ],
   "source": [
    "# Randomly knock off a certain percentage of the pixels\n",
    "per = 0.8 #fraction of pixels to knock out\n",
    "mask = rand(M,N)\n",
    "mask[ mask.> per ] .= 1\n",
    "mask[ mask.< 1 ]   .= 0\n",
    "f = (mask).*im;\n",
    "\n",
    "\n",
    "Gray.([im f])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A key idea again is that whatever we fill into the missing pixels, the result should overall \n",
    "be rather smooth. That is, we are using our prior assumption that we are dealing with a \n",
    "natural image. We discovered earlier that the heat flow will converge to a minimizer for the overall gradient energy. So how about trying a heat flow. But: The heat flow will also affect the points where we actually know the true data. So in each iteration, let's just put the given data back. Then the heat diffusion will only flow into the missing pixels and leave our known pixels unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[14]:3"
     ]
    }
   ],
   "source": [
    "tau = 1/4\n",
    "niter = 1000\n",
    "\n",
    "M,N = size(im)\n",
    "Laplace = zeros(size(im))\n",
    "m = 2:(N-1); n=2:(M-1)\n",
    "\n",
    "fheat = f\n",
    "for i = 1:niter\n",
    "    Laplace[n,m] = fheat[n.-1,m] + fheat[n.+1,m] + fheat[n,m.-1] + fheat[n,m.+1] - 4*fheat[n,m];\n",
    "    fheat += tau*Laplace;\n",
    "    fheat = fheat.*(mask.==0) + f; \n",
    "end\n",
    "\n",
    "Gray.([im f fheat])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us try the same with the total variation flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[15]:4"
     ]
    }
   ],
   "source": [
    "epsilon = 0.1\n",
    "tau = epsilon/5\n",
    "niter = 400\n",
    "\n",
    "M,N = size(im)\n",
    "fdx = zeros(size(im));  fdy = zeros(size(im)); div = zeros(size(im))\n",
    "j = 1:(N-1); k = 1:(M-1); m = 2:(N-1); n = 2:(M-1)\n",
    "\n",
    "#Gradient descent:\n",
    "ftv = f;\n",
    "for i = 1:niter\n",
    "    fdx[j,k] = ftv[j.+1,k] - ftv[j,k]; fdy[j,k] = ftv[j,k.+1] - ftv[j,k]\n",
    "    Gnorm = sqrt.( epsilon^2 .+ fdx.^2 .+ fdy.^2)\n",
    "    fdx = fdx ./ Gnorm; fdy = fdy ./ Gnorm\n",
    "    div[n,m] = fdx[n,m] - fdx[n.-1,m] + fdy[n,m] - fdy[n,m.-1]\n",
    "    ftv += tau*div\n",
    "    ftv = ftv.*(mask.==0) + f\n",
    "end\n",
    "\n",
    "Gray.([f fheat ftv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regularization for simultaneous denoising and impainting\n",
    "\n",
    "Let's get a little more realistic. A certain amount of noise will always be there. So suppose we have to inpaint, but the measurements we have are also noisy. Now we would not want to put the noisy data back in each step of the iteration. Because we also want to denoise. Let's first create some artificial data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "UndefVarError",
     "evalue": "UndefVarError: M not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: M not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope at In[16]:8"
     ]
    }
   ],
   "source": [
    "# Add some noise\n",
    "sigma    = 0.06  #Variance of the noise\n",
    "randmat  = randn(size(im));\n",
    "im_noisy = im + sigma*randmat\n",
    "\n",
    "# Randomly knock off a certain percentage of the pixels\n",
    "per = 0.15 #fraction of pixels to knock out\n",
    "mask = rand(M,N)\n",
    "mask[ mask.> per ] .= 1\n",
    "mask[ mask.< 1 ]   .= 0\n",
    "f = (mask).*im_noisy;\n",
    "\n",
    "\n",
    "Gray.([im f])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "e617ddf2-60ff-45f3-a292-349f871b8fa6"
    }
   },
   "source": [
    "Let us try to model the situation we are in. First, to make our life easier let us assume our image was\n",
    "just a single column\n",
    "of pixels. So we are observing some vector $f\\in\\mathbb{R}^N$, one column of $N$ pixels. Assume further that our observation contains a certain small amount of random noise, i.e. there is a noise vector $\\eta\\in\\mathbb{R}^N$. \n",
    "Finally, assume some pixels are just totally missing. We model this with a $N\\times N$ matrix $P$ that is diagonal (all off-diagonal entries are zero) and that on the diagonal only has zeros and ones. If say 40% of the diagonal entries are zero, we are missing 40% of the data in our observation, e.g.\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "1&0&0&0&0 \\\\\n",
    "0&0&0&0&0 \\\\\n",
    "0&0&1&0&0 \\\\\n",
    "0&0&0&1&0 \\\\\n",
    "0&0&0&0&0 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u_1\\\\\n",
    "u_2\\\\\n",
    "u_3\\\\\n",
    "u_4\\\\\n",
    "u_5\\\\\n",
    "\\end{pmatrix}\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "u_1\\\\\n",
    "0\\\\\n",
    "u_3\\\\\n",
    "u_4\\\\\n",
    "0\\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "\n",
    "Our model now reads\n",
    "\n",
    "$$\n",
    "f=P u + \\eta\n",
    "$$\n",
    "\n",
    "where $u$ is the unknow image we want to recover from observing the \"corrupted\" image $f$.\n",
    "\n",
    "One strategy is to solve the following minimization problem\n",
    "\n",
    "$$ \\min\\limits_{u} \\frac{1}{2}\\|f-P u\\|^2+ \\lambda J(u) $$\n",
    "\n",
    "Here $\\frac{1}{2}\\|f-P u\\|^2$ is called the data fidelity term. The intuition is that the noise $\\eta$ is relatively small, thus the Euclidean distance between $f$ and $Pu$ should be small, but not zero. The second term $\\lambda J(u) $\n",
    "is called the regularization term. In principle anything could have been in the places where the pixels are missing. However we are assuming that the unknown image was a natural image. Thus it makes sense to look for $u$ with small gradient energy $J(u)$ and $\\lambda$ is a parameter with which one can control how much a large gradient energy will be penalized. Other priors, say $J_{TV}^{\\epsilon}$, can also be used as regularization.\n",
    "\n",
    "- **Exercise:** We attempt to find $u$ via gradient descent, starting from some initial value $f^{(0)}$. Show that this leads to the iteration \n",
    "$$ f^{(k+1)} =f^{(k)} - \\tau(P^\\top(P f^{(k)}-f)\n",
    "-\\lambda\\nabla J(f^{(k)})) $$\n",
    "\n",
    "Note that in other applications there might be another degradation operator instead of $P$. For instance some operator that models how the image might be blurred due to camera shake, motion, or optics.\n",
    "\n",
    "An important question is how to choose $\\lambda$. Choosen too small not all noise may be removed, while choosing too large may result in destroying important feature of the image."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient energy regularization\n",
    "\n",
    "There is an important small detail in the implementation of the gradient descent. Images are in fact matrices, not vectors. Thus the projection $P$ is not a diagonal matrix, but some matrix (not diagonal) with entries being either 0 or 1 and which is multiplied entrywise (sic) with the image $u$. A diagonal matrix coincides with its transpose matrix. Thus in the implementation we still keep $P^\\top=P$, even though the matrix is not diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[17]:4"
     ]
    }
   ],
   "source": [
    "tau = 0.01\n",
    "lambda = 3\n",
    "niter = 500\n",
    "\n",
    "M,N = size(im)\n",
    "Laplace = zeros(size(im))\n",
    "m = 2:(N-1); n = 2:(M-1)\n",
    "\n",
    "fheat = f\n",
    "for i = 1:niter\n",
    "    Laplace[n,m] = fheat[n.-1,m] + fheat[n.+1,m] + fheat[n,m.-1] + fheat[n,m.+1] - 4*fheat[n,m]\n",
    "    fheat = fheat - tau*mask.*(mask.*fheat-f) + lambda*tau*Laplace\n",
    "end\n",
    "\n",
    "Gray.([f fheat])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total variation regularization\n",
    "\n",
    "Now we compare gradient energy regularization and TV regularization, i.e. instead of the regularization $\\lambda J(u)$ we now use $\\lambda J_{TV}^{(\\epsilon)}(u)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "ename": "BoundsError",
     "evalue": "BoundsError: attempt to access ()\n  at index [1]",
     "output_type": "error",
     "traceback": [
      "BoundsError: attempt to access ()\n  at index [1]",
      "",
      "Stacktrace:",
      " [1] indexed_iterate at .\\tuple.jl:63 [inlined] (repeats 2 times)",
      " [2] top-level scope at In[18]:5"
     ]
    }
   ],
   "source": [
    "epsilon = 0.2\n",
    "tau = epsilon/5\n",
    "niter = 500\n",
    "lambda  = 1\n",
    "\n",
    "M,N = size(im)\n",
    "fdx = zeros(size(im))\n",
    "fdy = zeros(size(im))\n",
    "div = zeros(size(im))\n",
    "j = 1:(N-1); k = 1:(M-1); m = 2:(N-1); n=2:(M-1)\n",
    "\n",
    "# Gradient descent:\n",
    "ftv = f\n",
    "for i = 1:niter\n",
    "    fdx[j,k] = ftv[j.+1,k] - ftv[j,k]\n",
    "    fdy[j,k] = ftv[j,k.+1] - ftv[j,k]\n",
    "    Gnorm = sqrt.(epsilon^2 .+ fdx.^2 .+ fdy.^2)\n",
    "    fdx = fdx ./ Gnorm\n",
    "    fdy = fdy ./ Gnorm\n",
    "    div[n,m] = fdx[n,m] - fdx[n.-1,m] + fdy[n,m] - fdy[n,m.-1]\n",
    "    ftv = ftv - tau*mask.*(mask.*ftv-f) + lambda*tau*div\n",
    "end\n",
    "\n",
    "Gray.([f fheat ftv])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References and further reading\n",
    "\n",
    "\n",
    "Numerical Analysis is a vast and vastly important subject. In the computer everything is discrete! Apparently, although this is a subject in its own, one needs a firm knowledge of the basics of classical real analysis.\n",
    "\n",
    "Gradient descent is as old as the mountains but still one of the most important numerical algorithms.  For instance, it is the crucial piece of mathematics underlying the error backpropagation step when training neural networks in machine learning. Methods on how to accelerate its convergence therefore are subject of current research.  \n",
    "\n",
    "The gradient energy sum/integral is connected to the name [Sergei Sobolev][3]. \n",
    "The total variation has been introduced to the area of signal processing by Rudin, [Osher][5] and Fetami. \n",
    "We will come back to the heat equation and its exact solution when we talk about the adventurous life and mathematics of the great French mathematician [Joseph Fourier][4], who started the field of harmonic analysis.\n",
    "\n",
    "Large parts of the material of this notebook are based on the text: Stephane Mallat, *A wavelet tour of signal processing* (Academic Press 2008), which is available as electronic copy through the [NUS library][1].\n",
    "\n",
    "\n",
    "[1]: http://linc.nus.edu.sg/search/\n",
    "[3]: https://en.wikipedia.org/wiki/Sergei_Sobolev\n",
    "[4]: https://en.wikipedia.org/wiki/Joseph_Fourier\n",
    "[5]: https://en.wikipedia.org/wiki/Stanley_Osher\n",
    "[6]: https://en.wikipedia.org/wiki/Unsharp_masking\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.1.0",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.0"
  },
  "nbpresent": {
   "slides": {
    "15923f4e-95da-4fe4-9c79-44ae9bce1463": {
     "id": "15923f4e-95da-4fe4-9c79-44ae9bce1463",
     "prev": "dc30df78-51cb-4108-b563-ad9752a623df",
     "regions": {
      "9b76d039-b5ae-4e6f-a647-77678264e88d": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "ffbe619b-be35-4211-9bf6-e499ad28d46d",
        "part": "whole"
       },
       "id": "9b76d039-b5ae-4e6f-a647-77678264e88d"
      }
     }
    },
    "19c72999-07cf-4c05-a738-b0d220ad480f": {
     "id": "19c72999-07cf-4c05-a738-b0d220ad480f",
     "prev": "ed0598a3-d2b7-4c06-8db8-2893b98acff1",
     "regions": {
      "ce6560ed-da7e-4465-a95b-c8f0fc8687a9": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "468c33ea-96a8-4f32-97f4-6991bda75ca1",
        "part": "whole"
       },
       "id": "ce6560ed-da7e-4465-a95b-c8f0fc8687a9"
      }
     }
    },
    "23c75bcc-2eeb-4022-bd7b-463a95d6162e": {
     "id": "23c75bcc-2eeb-4022-bd7b-463a95d6162e",
     "prev": "5b91d1b6-5514-40be-a36e-660707d03faa",
     "regions": {
      "da59bf4b-b4c2-4c58-b19c-c5e8d4b414b8": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "bd0284b6-c309-462e-8dea-13cd329f58b1",
        "part": "whole"
       },
       "id": "da59bf4b-b4c2-4c58-b19c-c5e8d4b414b8"
      }
     }
    },
    "2496634c-1b34-4268-905b-31464f96e9fb": {
     "id": "2496634c-1b34-4268-905b-31464f96e9fb",
     "prev": "bdc6095c-08fa-4e80-8df9-58762d0d8cde",
     "regions": {
      "9ec26f82-c2e0-4865-987f-e188df515c10": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "944dc4b0-5545-4874-bc32-ce6dc16d4b10",
        "part": "whole"
       },
       "id": "9ec26f82-c2e0-4865-987f-e188df515c10"
      }
     }
    },
    "3e4f95bf-8594-478f-9a1f-faacef2871a6": {
     "id": "3e4f95bf-8594-478f-9a1f-faacef2871a6",
     "prev": null,
     "regions": {
      "ff098ecf-72a1-4022-aeed-4f3d84f21362": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "8b6886e8-311d-4086-85bc-4c7ff294a9c5",
        "part": "whole"
       },
       "id": "ff098ecf-72a1-4022-aeed-4f3d84f21362"
      }
     }
    },
    "4d0d0eb4-9145-47e4-a3a8-fd7bfba87751": {
     "id": "4d0d0eb4-9145-47e4-a3a8-fd7bfba87751",
     "prev": "f480d3c6-ce84-4cde-978a-c3aa95471748",
     "regions": {
      "2ac103c3-4570-446c-b9df-98542b1455d1": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "413022cf-24fa-4d0b-ae72-e882eb484059",
        "part": "whole"
       },
       "id": "2ac103c3-4570-446c-b9df-98542b1455d1"
      }
     }
    },
    "5b91d1b6-5514-40be-a36e-660707d03faa": {
     "id": "5b91d1b6-5514-40be-a36e-660707d03faa",
     "prev": "d22d1508-aba6-4a26-b8c3-11429d9e69a9",
     "regions": {
      "798be81b-4d18-4ad7-855a-974e70f9cf62": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "e572123f-e8d4-41d0-b78e-a81a3158c340",
        "part": "whole"
       },
       "id": "798be81b-4d18-4ad7-855a-974e70f9cf62"
      }
     }
    },
    "5f5f9f2b-ddda-45a1-9d8b-9991d1ea2dbd": {
     "id": "5f5f9f2b-ddda-45a1-9d8b-9991d1ea2dbd",
     "prev": "a4e95e08-31ec-4a46-aac5-98210773b9cd",
     "regions": {
      "69cd456e-7510-4f0f-ada4-6914344a6a81": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "c5194dd5-2c9d-4998-97b7-c246184f0703",
        "part": "whole"
       },
       "id": "69cd456e-7510-4f0f-ada4-6914344a6a81"
      }
     }
    },
    "64c6eff1-2da8-43bb-90b3-619a5a651e9a": {
     "id": "64c6eff1-2da8-43bb-90b3-619a5a651e9a",
     "prev": "d1e8fca8-673c-4766-9f67-e077714cdc49",
     "regions": {
      "57a409d2-5184-4c78-b4cd-4c0fa2a67978": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "83ac1886-9598-49d9-88de-fff6898509b4",
        "part": "whole"
       },
       "id": "57a409d2-5184-4c78-b4cd-4c0fa2a67978"
      }
     }
    },
    "727c0763-759f-45bb-b581-95415116b9d5": {
     "id": "727c0763-759f-45bb-b581-95415116b9d5",
     "prev": "b721c89e-c27e-4a14-b513-423b03cd874e",
     "regions": {
      "a6df0843-9ff2-4a84-83e9-41352deec7ed": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "4e426eef-6ef7-4067-b420-6aabd5b87be9",
        "part": "whole"
       },
       "id": "a6df0843-9ff2-4a84-83e9-41352deec7ed"
      }
     }
    },
    "a4e95e08-31ec-4a46-aac5-98210773b9cd": {
     "id": "a4e95e08-31ec-4a46-aac5-98210773b9cd",
     "prev": "3e4f95bf-8594-478f-9a1f-faacef2871a6",
     "regions": {
      "4c4fa8c4-02a4-428d-92b0-f90fee1981e8": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "a3959ad8-a955-4972-b11f-8e3cdd87b8c5",
        "part": "whole"
       },
       "id": "4c4fa8c4-02a4-428d-92b0-f90fee1981e8"
      }
     }
    },
    "b721c89e-c27e-4a14-b513-423b03cd874e": {
     "id": "b721c89e-c27e-4a14-b513-423b03cd874e",
     "prev": "d9abbce4-b3dc-4c82-bb7e-a0cf6fb37a3f",
     "regions": {
      "33c95330-14a1-411e-9e4a-3edee19a4ad6": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "f16ed3fa-09d9-49dc-9be2-50fc8ffa4fbc",
        "part": "whole"
       },
       "id": "33c95330-14a1-411e-9e4a-3edee19a4ad6"
      }
     }
    },
    "bdc6095c-08fa-4e80-8df9-58762d0d8cde": {
     "id": "bdc6095c-08fa-4e80-8df9-58762d0d8cde",
     "prev": "5f5f9f2b-ddda-45a1-9d8b-9991d1ea2dbd",
     "regions": {
      "d50e7f91-7bf5-4a25-b292-9f692267d199": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.02138245347512345,
        "y": 0.10303835928598547
       },
       "content": {
        "cell": "4740a16f-2f05-4540-8cdf-5bb002db8236",
        "part": "whole"
       },
       "id": "d50e7f91-7bf5-4a25-b292-9f692267d199"
      }
     }
    },
    "d1e8fca8-673c-4766-9f67-e077714cdc49": {
     "id": "d1e8fca8-673c-4766-9f67-e077714cdc49",
     "prev": "d5ee9599-4d7e-47e1-b7d3-653851c580e3",
     "regions": {
      "68852ed9-5a66-4aa5-a1bc-fbd612b60014": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "a76027d7-fc9d-4533-8592-c70fe9be8263",
        "part": "whole"
       },
       "id": "68852ed9-5a66-4aa5-a1bc-fbd612b60014"
      }
     }
    },
    "d2101dc4-db3f-4a1f-bcbf-43d15c54efde": {
     "id": "d2101dc4-db3f-4a1f-bcbf-43d15c54efde",
     "prev": "2496634c-1b34-4268-905b-31464f96e9fb",
     "regions": {
      "aa84db64-2e2f-4aaf-ad0e-5657e14b8e8b": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "d00ebc9a-5c1b-4b66-b9ed-93b2f10aafa1",
        "part": "whole"
       },
       "id": "aa84db64-2e2f-4aaf-ad0e-5657e14b8e8b"
      }
     }
    },
    "d22d1508-aba6-4a26-b8c3-11429d9e69a9": {
     "id": "d22d1508-aba6-4a26-b8c3-11429d9e69a9",
     "prev": "64c6eff1-2da8-43bb-90b3-619a5a651e9a",
     "regions": {
      "bd8c5c3d-8cb1-46ef-829f-cbd1140f4ccc": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "2471b46f-451d-4599-b11e-a94e1a5540cd",
        "part": "whole"
       },
       "id": "bd8c5c3d-8cb1-46ef-829f-cbd1140f4ccc"
      }
     }
    },
    "d5ee9599-4d7e-47e1-b7d3-653851c580e3": {
     "id": "d5ee9599-4d7e-47e1-b7d3-653851c580e3",
     "prev": "d2101dc4-db3f-4a1f-bcbf-43d15c54efde",
     "regions": {
      "5eb41739-0640-4792-bcdf-697234786518": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "122607dd-b961-4dfe-aab1-1e1491005bda",
        "part": "whole"
       },
       "id": "5eb41739-0640-4792-bcdf-697234786518"
      }
     }
    },
    "d9abbce4-b3dc-4c82-bb7e-a0cf6fb37a3f": {
     "id": "d9abbce4-b3dc-4c82-bb7e-a0cf6fb37a3f",
     "prev": "19c72999-07cf-4c05-a738-b0d220ad480f",
     "regions": {
      "a4a5c7fa-98e2-498c-bcca-af1073032a4a": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "e617ddf2-60ff-45f3-a292-349f871b8fa6",
        "part": "whole"
       },
       "id": "a4a5c7fa-98e2-498c-bcca-af1073032a4a"
      }
     }
    },
    "dc30df78-51cb-4108-b563-ad9752a623df": {
     "id": "dc30df78-51cb-4108-b563-ad9752a623df",
     "prev": "4d0d0eb4-9145-47e4-a3a8-fd7bfba87751",
     "regions": {
      "6fc1ce7f-0208-4a31-b356-c76df944d66e": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "ce446957-d62d-4374-85f5-4f37783ebaac",
        "part": "whole"
       },
       "id": "6fc1ce7f-0208-4a31-b356-c76df944d66e"
      }
     }
    },
    "ed0598a3-d2b7-4c06-8db8-2893b98acff1": {
     "id": "ed0598a3-d2b7-4c06-8db8-2893b98acff1",
     "prev": "15923f4e-95da-4fe4-9c79-44ae9bce1463",
     "regions": {
      "e2775394-474c-4d9b-8c17-6e217f3b988b": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "033f18a7-1991-4b91-9797-577650effb58",
        "part": "whole"
       },
       "id": "e2775394-474c-4d9b-8c17-6e217f3b988b"
      }
     }
    },
    "f480d3c6-ce84-4cde-978a-c3aa95471748": {
     "id": "f480d3c6-ce84-4cde-978a-c3aa95471748",
     "prev": "23c75bcc-2eeb-4022-bd7b-463a95d6162e",
     "regions": {
      "efb36e6a-cad1-452c-9ac7-b4f51f322c49": {
       "attrs": {
        "height": 0.8,
        "width": 0.8,
        "x": 0.1,
        "y": 0.1
       },
       "content": {
        "cell": "b97ebe18-6e07-4bfc-829d-d11d56c55c17",
        "part": "whole"
       },
       "id": "efb36e6a-cad1-452c-9ac7-b4f51f322c49"
      }
     }
    }
   },
   "themes": {}
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}