{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1.1 \n",
    "\n",
    "## Import Libraries\n",
    "Python requires importing libraries and functions you need to access specific tools like science (scipy), linear algebra (numpy), and graphics (matplotlib). These libraries can be installed using the ```pip``` command line tool. Alternatively you can install an python distribution like [Anaconda](https://www.continuum.io/downloads) or [Canopy](https://www.enthought.com/products/canopy/)  which have these and many other standard package pre-installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt # plotting\n",
    "from skimage.io import imread # read in images\n",
    "import numpy as np # linear algebra / matrices\n",
    "# make the notebook interactive\n",
    "from ipywidgets import interact, interactive, fixed \n",
    "import ipywidgets as widgets #add new widgets\n",
    "from IPython.display import display\n",
    "class idict(dict):\n",
    "    def __init__(self,*args,**kwargs) : dict.__init__(self,*args,**kwargs) \n",
    "    def __str__(self): return 'ImageDictionary'\n",
    "    def __repr__(self): return 'ImageDictionary'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Images\n",
    "Make sure you extract the ```matlab.zip``` file to the same directory as this notebook so there is a ```data/``` directory (or fix the paths after the ```imread``` command"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=imread('data/scroll.tif')\n",
    "b=imread('data/wood.tif')\n",
    "c=imread('data/asphalt_gray.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Showing images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# setup the plotting environment\n",
    "plt.imshow(a, cmap = 'gray'); # show a single image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subplots\n",
    "Here we show multiple subplots within a single figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15,5))\n",
    "ax1.imshow(a, cmap = 'gray')\n",
    "ax2.imshow(b, cmap = 'gray')\n",
    "ax3.imshow(c, cmap = 'gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute SNR\n",
    "We can compute the SNR by looking at the ratio of the mean to the standard deviation in a region that is supposed to be constant\n",
    "\n",
    "$$ SNR = \\frac{\\mu_{img}}{\\sigma_{img}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify the region\n",
    "x1 = 75\n",
    "x2 = 200\n",
    "y1 = 875\n",
    "y2 = 1000\n",
    "\n",
    "# extract a sub image\n",
    "subA1=a[x1:x2,y1:y2];\n",
    "snrA1=np.mean(subA1)/np.std(subA1) # compute the snr\n",
    "print(\"SNR for A_1 is {}\".format(snrA1))\n",
    "plt.imshow(subA1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## __Tasks__\n",
    "1. Find a second region in a\n",
    "1. Repeat the procedure with images b and c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d=np.mean(imread('data/testpattern.png'),2)\n",
    "plt.imshow(d, cmap= 'gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import randn\n",
    "\n",
    "def show_noisy_images(scale_100, scale_10, scale_5, scale_2, scale_1):\n",
    "    fig, ((ax1,ax2,ax3),(ax4,ax5,ax6)) = plt.subplots(2, 3, figsize=(8.5,5))\n",
    "    ax1.imshow(d)\n",
    "    ax1.set_title('Original')\n",
    "\n",
    "    d_snr100=d+scale_100*randn(*d.shape);\n",
    "    ax2.imshow(d_snr100)\n",
    "    ax2.set_title('SNR 100')\n",
    "\n",
    "\n",
    "    d_snr10=d+scale_10*randn(*d.shape);\n",
    "    ax3.imshow(d_snr10)\n",
    "    ax3.set_title('SNR 10')\n",
    "\n",
    "    scale = 100 \n",
    "    d_snr5=d+scale_5*randn(*d.shape);\n",
    "    ax4.imshow(d_snr5)\n",
    "    ax4.set_title('SNR 5')\n",
    "\n",
    "    scale = 1000 \n",
    "    d_snr2=d+scale_2*randn(*d.shape);\n",
    "    ax5.imshow(d_snr100)\n",
    "    ax5.set_title('SNR 2')\n",
    "\n",
    "    scale = 5000 \n",
    "    d_snr1=d+scale_1*randn(*d.shape);\n",
    "    ax6.imshow(d_snr1)\n",
    "    ax6.set_title('SNR 1')\n",
    "    return idict({1: d_snr1, 2: d_snr2, 5: d_snr5, 10: d_snr10, 100: d_snr100})\n",
    "noisy_images = interactive(show_noisy_images, scale_100 = (0.0,1.0), scale_10 = (0.0,10.0), scale_5 = (0.0,50.0), scale_2 = (0.0,100.0), scale_1 = (0.0,200.0))\n",
    "display(noisy_images)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter the images\n",
    "You can read about the standard filters in scipy by looking at the documentation in http://docs.scipy.org/doc/scipy-0.14.0/reference/ndimage.html#module-scipy.ndimage.filters alternatively more (different) filters are available using OpenCV for the more advanced students\n",
    "\n",
    "## Uniform Filters\n",
    "The specific documentation on the filter is as below\n",
    "http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.filters.uniform_filter.html#scipy.ndimage.filters.uniform_filter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage.filters import uniform_filter\n",
    "# Size of the filter window\n",
    "N=[3,5,7];\n",
    "# Images\n",
    "noisy_image_list = noisy_images.result\n",
    "for i,filter_size in enumerate(N):\n",
    "    fig, all_axes = plt.subplots(5, 2, figsize=(6,22))\n",
    "    for ((snr,img),(ax1,ax2)) in zip(noisy_image_list.items(),all_axes):\n",
    "        ax1.imshow(img, cmap='gray')\n",
    "        ax1.set_title(\"Raw, SNR:{}\".format(snr))\n",
    "        ax2.imshow(uniform_filter(img,filter_size), cmap='gray')\n",
    "        ax2.set_title(\"N:{}, SNR:{}\".format(filter_size,snr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Median Filter\n",
    "The median filter function can be found here, complete the same exercise as before using this instead\n",
    "http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.filters.median_filter.html#scipy.ndimage.filters.median_filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# insert your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 - Diffusion filter\n",
    "Tune the filter parameters for the three experiment images "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def anisodiff(img,niter=1,kappa=50,gamma=0.1,step=(1.,1.),sigma=0.0,option=1,ploton=False):\n",
    "    \"\"\"\n",
    "    Anisotropic diffusion.\n",
    "\n",
    "    Usage:\n",
    "    imgout = anisodiff(im, niter, kappa, gamma, option)\n",
    "\n",
    "    Arguments:\n",
    "            img    - input image\n",
    "            niter  - number of iterations\n",
    "            kappa  - conduction coefficient 20-100 ?\n",
    "            gamma  - max value of .25 for stability\n",
    "            step   - tuple, the distance between adjacent pixels in (y,x)\n",
    "            option - 1 Perona Malik diffusion equation No 1\n",
    "                     2 Perona Malik diffusion equation No 2\n",
    "            ploton - if True, the image will be plotted on every iteration\n",
    "\n",
    "    Returns:\n",
    "            imgout   - diffused image.\n",
    "\n",
    "    kappa controls conduction as a function of gradient.  If kappa is low\n",
    "    small intensity gradients are able to block conduction and hence diffusion\n",
    "    across step edges.  A large value reduces the influence of intensity\n",
    "    gradients on conduction.\n",
    "\n",
    "    gamma controls speed of diffusion (you usually want it at a maximum of\n",
    "    0.25)\n",
    "\n",
    "    step is used to scale the gradients in case the spacing between adjacent\n",
    "    pixels differs in the x and y axes\n",
    "\n",
    "    Diffusion equation 1 favours high contrast edges over low contrast ones.\n",
    "    Diffusion equation 2 favours wide regions over smaller ones.\n",
    "\n",
    "    Reference: \n",
    "    P. Perona and J. Malik. \n",
    "    Scale-space and edge detection using ansotropic diffusion.\n",
    "    IEEE Transactions on Pattern Analysis and Machine Intelligence, \n",
    "    12(7):629-639, July 1990.\n",
    "\n",
    "    Original MATLAB code by Peter Kovesi  \n",
    "    School of Computer Science & Software Engineering\n",
    "    The University of Western Australia\n",
    "    pk @ csse uwa edu au\n",
    "    <http://www.csse.uwa.edu.au>\n",
    "\n",
    "    Translated to Python and optimised by Alistair Muldal\n",
    "    Department of Pharmacology\n",
    "    University of Oxford\n",
    "    <alistair.muldal@pharm.ox.ac.uk>\n",
    "\n",
    "    June 2000  original version.       \n",
    "    March 2002 corrected diffusion eqn No 2.\n",
    "    July 2012 translated to Python\n",
    "\n",
    "    February 2020 Revised for Python 3.6 / A. Kaestner\n",
    "    \"\"\"\n",
    "\n",
    "    # ...you could always diffuse each color channel independently if you\n",
    "    # really want\n",
    "    if img.ndim == 3:\n",
    "        warnings.warn(\"Only grayscale images allowed, converting to 2D matrix\")\n",
    "        img = img.mean(2)\n",
    "\n",
    "    # initialize output array\n",
    "    img = img.astype('float32')\n",
    "    imgout = img.copy()\n",
    "\n",
    "    # initialize some internal variables\n",
    "    deltaS = np.zeros_like(imgout)\n",
    "    deltaE = deltaS.copy()\n",
    "    NS = deltaS.copy()\n",
    "    EW = deltaS.copy()\n",
    "    gS = np.ones_like(imgout)\n",
    "    gE = gS.copy()\n",
    "\n",
    "    # create the plot figure, if requested\n",
    "    if ploton:\n",
    "        import matplotlib.pyplot as plt\n",
    "        from time import sleep\n",
    "\n",
    "        plt.figure(figsize=(20,5.5),num=\"Anisotropic diffusion\")\n",
    "        plt.subplot(1,3,1)\n",
    "        plt.imshow(img,interpolation='nearest')\n",
    "        plt.title('Original')\n",
    "        plt.colorbar()\n",
    "\n",
    "        \n",
    "    for ii in np.arange(0,niter):\n",
    "        smoothimgout = imgout\n",
    "        \n",
    "        if sigma != 0 :\n",
    "            smoothimgout = imgout         ###### Introduce gradient smoothing here \n",
    "            \n",
    "        # calculate the diffs\n",
    "        deltaS[:-1,: ] = np.diff(smoothimgout,axis=0)\n",
    "        deltaE[: ,:-1] = np.diff(smoothimgout,axis=1)\n",
    "        \n",
    "\n",
    "        # conduction gradients (only need to compute one per dim!)\n",
    "        if option == 1:\n",
    "            gS = np.exp(-(deltaS/kappa)**2.)/step[0]\n",
    "            gE = np.exp(-(deltaE/kappa)**2.)/step[1]\n",
    "        elif option == 2:\n",
    "            gS = 1./(1.+(deltaS/kappa)**2.)/step[0]\n",
    "            gE = 1./(1.+(deltaE/kappa)**2.)/step[1]\n",
    "\n",
    "        # update matrices\n",
    "        E = gE*deltaE\n",
    "        S = gS*deltaS\n",
    "\n",
    "        # subtract a copy that has been shifted 'North/West' by one\n",
    "        # pixel. don't as questions. just do it. trust me.\n",
    "        NS[:] = S\n",
    "        EW[:] = E\n",
    "        NS[1:,:] -= S[:-1,:]\n",
    "        EW[:,1:] -= E[:,:-1]\n",
    "\n",
    "        # update the image\n",
    "        imgout += gamma*(NS+EW)\n",
    "\n",
    "    if ploton:\n",
    "        iterstring = \"Iteration %i\" %(ii+1)\n",
    "        plt.subplot(1,3,2)\n",
    "        plt.imshow(imgout)\n",
    "        plt.title(iterstring)\n",
    "        \n",
    "        plt.subplot(1,3,3)\n",
    "        plt.imshow(img-imgout)\n",
    "        plt.title('Difference before - after')\n",
    "\n",
    "\n",
    "    return imgout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3.1 Find filter parameters\n",
    "In this exercise you need to find the parameters to suppress the noise in the image. __Hint:__ compute the histogram of the gradient image, this will guide you. A very important parameter is $\\kappa$ that acts a sensitivity threshold on the gradient contribution. The solution time  i.e. $niter \\times step$ also plays a role in the filter strength. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "niter=10 \n",
    "kappa=50  \n",
    "gamma=0.1 \n",
    "step=(0.25, 0.25) \n",
    "option=2  # select weighing equation\n",
    "ploton=True\n",
    "res=anisodiff(img=a, niter=niter, kappa=kappa, gamma=gamma, step=step, option=option,ploton=ploton)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Part 3.2 Introduce some regularization (more advanced)\n",
    "The original Perona Malik filter is very sensitive to noise. A way improve this is to smooth the image before computing the gradient. \n",
    "\n",
    "- Locate the line ```###### Introduce gradient smoothing here``` in the filter function and add a line to filter the image that is used to compute the gradient using a Gauss filter.\n",
    "- Use filtered image for the gradient but keep ```imgout``` to update the filter iteration.\n",
    "- Use the parameter ```sigma``` to the function definition to control the width of the Gauss kernel.\n",
    "- Compare the preformance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "niter=10 \n",
    "kappa=50  \n",
    "gamma=0.1 \n",
    "sigma=1.0\n",
    "step=(0.25, 0.25) \n",
    "option=2  # select weighing equation\n",
    "ploton=True\n",
    "res=anisodiff(img=a, niter=niter, kappa=kappa, gamma=gamma, step=step, option=option,sigma=sigma,ploton=ploton)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}