{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive Image Processing with Numba and Bokeh\n",
    "\n",
    "This demo shows off how interactive image processing can be done in the notebook, using [Numba](http://numba.pydata.org) for numerics, [Bokeh](http://bokeh.pydata.org) for plotting, and Ipython interactors for widgets. The demo runs entirely inside the Ipython notebook, with no Bokeh server required.\n",
    "\n",
    "Numba must be installed in order to run this demo. To run, click on, `Cell->Run All` in the top menu, then scroll down to individual examples and play around with their controls. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "from bokeh.plotting import figure, show, output_notebook\n",
    "from bokeh.models import GlyphRenderer, LinearColorMapper\n",
    "from bokeh.io import push_notebook\n",
    "from numba import jit, njit\n",
    "\n",
    "from ipywidgets import interact\n",
    "import numpy as np\n",
    "import scipy.misc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Blur\n",
    "\n",
    "This first section demonstrates performing a simple Gaussian blur on an image. It presents the image, as well as a slider that controls how much blur is applied. Numba is used to compile the python blur kernel, which is invoked when the user modifies the slider. \n",
    "\n",
    "*Note:* This simple example does not handle the edge case, so the edge of the image will remain unblurred as the slider is increased. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# smaller image\n",
    "img_blur = (scipy.misc.ascent()[::-1,:]/255.0)[:250, :250].copy(order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "palette = ['#%02x%02x%02x' %(i,i,i) for i in range(256)]\n",
    "width, height = img_blur.shape\n",
    "p_blur = figure(x_range=(0, width), y_range=(0, height))\n",
    "r_blur = p_blur.image(image=[img_blur], x=[0], y=[0], dw=[width], dh=[height], palette=palette, name='blur')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def blur(outimg, img, amt):\n",
    "    iw, ih = img.shape\n",
    "    for i in range(amt, iw-amt):\n",
    "        for j in range(amt, ih-amt):\n",
    "            px = 0.\n",
    "            for w in range(-amt//2, amt//2):\n",
    "                for h in range(-amt//2, amt//2):\n",
    "                    px += img[i+w, j+h]\n",
    "            outimg[i, j]= px/(amt*amt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def update(i=0):\n",
    "    level = 2*i + 1\n",
    "    \n",
    "    out = img_blur.copy()\n",
    "    \n",
    "    ts = timer()\n",
    "    blur(out, img_blur, level)\n",
    "    te = timer()\n",
    "    print('blur takes:', te - ts)\n",
    "    \n",
    "    renderer = p_blur.select(dict(name=\"blur\", type=GlyphRenderer))\n",
    "    r_blur.data_source.data['image'] = [out]\n",
    "    push_notebook(handle=t_blur)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_blur = show(p_blur)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(update, i=(0, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3x3 Image Kernels\n",
    "\n",
    "Many image processing filters can be expressed as 3x3 matrices. This more sophisticated example demonstrates how numba can be used to compile kernels for arbitrary 3x3 kernels, and then provides serveral predefined kernels for the user to experiment with. \n",
    "\n",
    "The UI presents the image to process (along with a dropdown to select a different image) as well as a dropdown that lets the user select which kernel to apply. Additioanlly there are sliders the permit adjustment to the bias and scale of the final greyscale image. \n",
    "\n",
    "*Note:* Right now, adjusting the scale and bias are not as efficient as possible, because the update function always also applies the kernel (even if it has not changed). A better implementation might have a class that keeps track of the current kernal and output image so that bias and scale can be applied by themselves. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def getitem(img, x, y):\n",
    "    w, h = img.shape\n",
    "    if x >= w:\n",
    "        x = w - 1 - (x - w)\n",
    "    if y >= h:\n",
    "        y = h - 1 - (y - h)\n",
    "    return img[x, y]\n",
    " \n",
    "def filter_factory(kernel):\n",
    "    ksum = np.sum(kernel)\n",
    "    if ksum == 0:\n",
    "        ksum = 1\n",
    "    k9 = kernel / ksum\n",
    " \n",
    "    @jit\n",
    "    def kernel_apply(img, out, x, y):\n",
    "        tmp = 0\n",
    "        for i in range(3):\n",
    "            for j in range(3):\n",
    "                tmp += img[x+i-1, y+j-1] * k9[i, j]\n",
    "        out[x, y] = tmp\n",
    " \n",
    "    @jit\n",
    "    def kernel_apply_edge(img, out, x, y):\n",
    "        tmp = 0\n",
    "        for i in range(3):\n",
    "            for j in range(3):\n",
    "                tmp += getitem(img, x+i-1, y+j-1) * k9[i, j]\n",
    "        out[x, y] = tmp\n",
    " \n",
    "    @jit\n",
    "    def kernel_k9(img, out):\n",
    "        # Loop through all internals\n",
    "        for x in range(1, img.shape[0] -1):\n",
    "            for y in range(1, img.shape[1] -1):\n",
    "                kernel_apply(img, out, x, y)\n",
    " \n",
    "        # Loop through all the edges\n",
    "        for x in range(img.shape[0]):\n",
    "            kernel_apply_edge(img, out, x, 0)\n",
    "            kernel_apply_edge(img, out, x, img.shape[1] - 1)\n",
    " \n",
    "        for y in range(img.shape[1]):\n",
    "            kernel_apply_edge(img, out, 0, y)\n",
    "            kernel_apply_edge(img, out, img.shape[0] - 1, y)\n",
    " \n",
    "    return kernel_k9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "average = np.array([\n",
    "    [1, 1, 1],\n",
    "    [1, 1, 1],\n",
    "    [1, 1, 1],\n",
    "], dtype=np.float32)\n",
    "\n",
    "sharpen = np.array([\n",
    "    [-1, -1, -1],\n",
    "    [-1, 12, -1],\n",
    "    [-1, -1, -1],\n",
    "], dtype=np.float32)\n",
    "\n",
    "edge = np.array([\n",
    "    [ 0, -1,  0],\n",
    "    [-1,  4, -1],\n",
    "    [ 0, -1,  0],\n",
    "], dtype=np.float32)\n",
    "\n",
    "edge_h = np.array([\n",
    "    [ 0,  0,  0],\n",
    "    [-1,  2, -1],\n",
    "    [ 0,  0,  0],\n",
    "], dtype=np.float32)\n",
    "\n",
    "edge_v = np.array([\n",
    "    [0, -1, 0],\n",
    "    [0,  2, 0],\n",
    "    [0, -1, 0],\n",
    "], dtype=np.float32)\n",
    "\n",
    "gradient_h = np.array([\n",
    "    [-1, -1, -1],\n",
    "    [ 0,  0,  0],\n",
    "    [ 1,  1,  1],\n",
    "], dtype=np.float32)\n",
    "\n",
    "gradient_v = np.array([\n",
    "    [-1, 0, 1],\n",
    "    [-1, 0, 1],\n",
    "    [-1, 0, 1],\n",
    "], dtype=np.float32)\n",
    "\n",
    "sobol_h = np.array([\n",
    "    [ 1,  2,  1],\n",
    "    [ 0,  0,  0],\n",
    "    [-1, -2, -1],\n",
    "], dtype=np.float32)\n",
    "\n",
    "sobol_v = np.array([\n",
    "    [-1, 0, 1],\n",
    "    [-2, 0, 2],\n",
    "    [-1, 0, 1],\n",
    "], dtype=np.float32)\n",
    " \n",
    "emboss = np.array([    \n",
    "    [-2, -1, 0],\n",
    "    [-1,  1, 1],\n",
    "    [ 0,  1, 2],\n",
    "], dtype=np.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kernels = {\n",
    "    \"average\"               : filter_factory(average),\n",
    "    \"sharpen\"               : filter_factory(sharpen),\n",
    "    \"edge (both)\"           : filter_factory(edge),\n",
    "    \"edge (horizontal)\"     : filter_factory(edge_h),\n",
    "    \"edge (vertical)\"       : filter_factory(edge_v),\n",
    "    \"gradient (horizontal)\" : filter_factory(gradient_h),\n",
    "    \"gradient (vertical)\"   : filter_factory(gradient_v),\n",
    "    \"sobol (horizontal)\"    : filter_factory(sobol_h),\n",
    "    \"sobol (vertical)\"      : filter_factory(sobol_v),\n",
    "    \"emboss\"                : filter_factory(emboss),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "images = {\n",
    "    \"ascent\" : np.copy(scipy.misc.ascent().astype(np.float32)[::-1, :]),\n",
    "    \"face\"   : np.copy(scipy.misc.face(gray=True).astype(np.float32)[::-1, :]),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "palette = ['#%02x%02x%02x' %(i,i,i) for i in range(256)]\n",
    "cm = LinearColorMapper(palette=palette, low=0, high=256)\n",
    "width, height = images['ascent'].shape\n",
    "p_kernel = figure(x_range=(0, width), y_range=(0, height))\n",
    "r_kernel = p_kernel.image(image=[images['ascent']], x=[0], y=[0], dw=[width], dh=[height], color_mapper=cm, name=\"kernel\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def update(image=\"ascent\", kernel_name=\"none\", scale=100, bias=0):\n",
    "    global _last_kname\n",
    "    global _last_out\n",
    "    \n",
    "    img_kernel = images.get(image)\n",
    "\n",
    "    kernel = kernels.get(kernel_name, None)\n",
    "    if kernel == None:\n",
    "        out = np.copy(img_kernel)\n",
    "\n",
    "    else:\n",
    "        out = np.zeros_like(img_kernel)\n",
    "\n",
    "        ts = timer()\n",
    "        kernel(img_kernel, out)\n",
    "        te = timer()\n",
    "        print('kernel takes:', te - ts)\n",
    "\n",
    "    out *= scale / 100.0\n",
    "    out += bias\n",
    "\n",
    "    r_kernel.data_source.data['image'] = [out]\n",
    "    push_notebook(handle=t_kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_kernel = show(p_kernel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "knames = [\"none\"] + sorted(kernels.keys())\n",
    "interact(update, image=[\"ascent\" ,\"face\"], kernel_name=knames, scale=(10, 100, 10), bias=(0, 255))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wavelet Decomposition\n",
    "\n",
    "This last example demostrates a Haar wavelet decomposition using a Numba-compiled function. Play around with the slider to see differnet levels of decomposition of the image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def wavelet_decomposition(img, tmp):\n",
    "    \"\"\"\n",
    "    Perform inplace wavelet decomposition on `img` with `tmp` as\n",
    "    a temporarily buffer.\n",
    "\n",
    "    This is a very simple wavelet for demonstration\n",
    "    \"\"\"\n",
    "    w, h = img.shape\n",
    "    halfwidth, halfheight = w//2, h//2\n",
    " \n",
    "    lefthalf, righthalf = tmp[:halfwidth, :], tmp[halfwidth:, :]\n",
    " \n",
    "    # Along first dimension\n",
    "    for x in range(halfwidth):\n",
    "        for y in range(h):\n",
    "            lefthalf[x, y] = (img[2 * x, y] + img[2 * x + 1, y]) / 2\n",
    "            righthalf[x, y] = img[2 * x, y] - img[2 * x + 1, y]\n",
    " \n",
    "    # Swap buffer\n",
    "    img, tmp = tmp, img\n",
    "    tophalf, bottomhalf = tmp[:, :halfheight], tmp[:, halfheight:]\n",
    " \n",
    "    # Along second dimension\n",
    "    for y in range(halfheight):\n",
    "        for x in range(w):\n",
    "            tophalf[x, y] = (img[x, 2 * y] + img[x, 2 * y + 1]) / 2\n",
    "            bottomhalf[x, y] = img[x, 2 * y] - img[x, 2 * y + 1]\n",
    " \n",
    "    return halfwidth, halfheight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img_wavelet = np.copy(scipy.misc.face(gray=True)[::-1, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "palette = ['#%02x%02x%02x' %(i,i,i) for i in range(256)]\n",
    "width, height = img_wavelet.shape\n",
    "p_wavelet = figure(x_range=(0, width), y_range=(0, height))\n",
    "r_wavelet = p_wavelet.image(image=[img_wavelet], x=[0], y=[0], dw=[width], dh=[height], palette=palette, name=\"wavelet\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def update(level=0):\n",
    "\n",
    "    out = np.copy(img_wavelet)\n",
    "    tmp = np.zeros_like(img_wavelet)\n",
    "\n",
    "    ts = timer()\n",
    "    hw, hh = img_wavelet.shape\n",
    "    while level > 0 and hw > 1 and hh > 1:\n",
    "        hw, hh = wavelet_decomposition(out[:hw, :hh], tmp[:hw, :hh])\n",
    "        level -= 1\n",
    "    te = timer()\n",
    "    print('wavelet takes:', te - ts)\n",
    "\n",
    "    r_wavelet.data_source.data['image'] = [out]\n",
    "    push_notebook(handle=t_wavelet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "t_wavelet = show(p_wavelet)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(update, level=(0, 7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}