scipy-2014/1_image_filters.ipynb

{
 "metadata": {
  "celltoolbar": "Slideshow",
  "name": "",
  "signature": "sha256:bcccf9281ea745f8eb64bccf56f7c8e3fbcca1dd0a1148a88dbdfe5c1347c9c6"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Image filtering"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Image filtering theory"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Filtering is one of the most basic and common image operations in image processing. You can filter an image to remove noise or to enhance features; the filtered image could be the desired result or just a preprocessing step. Regardless, filtering is an important topic to understand."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Local filtering"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Let's get started by setting the colormap to grayscale and turning off pixel interpolation, as discussed in the previous section."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "import numpy as np\n",
      "\n",
      "plt.rcParams['image.cmap'] = 'gray'\n",
      "plt.rcParams['image.interpolation'] = 'none'"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The \"local\" in local filtering simply means that a pixel is adjusted by values in some surrounding neighborhood. These surrounding elements are identified or weighted based on a \"footprint\", \"structuring element\", or \"kernel\".\n",
      "\n",
      "Let's go to back to basics and look at a 1D step-signal"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "step_signal = np.zeros(100)\n",
      "step_signal[50:] = 10\n",
      "plt.plot(step_signal)\n",
      "plt.margins(0.1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Now add some (very) artificial noise to this signal:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "noisy_signal = np.copy(step_signal)\n",
      "noisy_signal[::2] += 1\n",
      "plt.plot(noisy_signal)\n",
      "plt.margins(0.1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The simplest way to recover something that looks like the original image is to take the average between neighboring pixels:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Take the mean of neighboring pixels\n",
      "smooth_signal = (noisy_signal[:-1] + noisy_signal[1:]) / 2.0\n",
      "plt.plot(smooth_signal)\n",
      "plt.margins(0.1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This same concept, nearest-neighbor averages, can be expressed as a convolution with an averaging kernel:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# The neighboring\n",
      "mean_kernel = np.array([1, 1]) / 2.0\n",
      "smooth_signal = np.convolve(noisy_signal, mean_kernel, mode='valid')\n",
      "plt.plot(smooth_signal)\n",
      "plt.margins(0.1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Local filtering of images"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Now let's apply this to a real image. Let's start with an incredibly simple image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "\n",
      "bright_square = np.zeros((7, 7), dtype=float)\n",
      "bright_square[2:5, 2:5] = 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This gives the values below:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print bright_square"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "and looks like a white square centered on a black square:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.imshow(bright_square);"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "The mean filter"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "For our first example of a filter, consider the following filtering array, which we'll call a \"mean kernel\". For each pixel, a kernel defines which neighboring pixels to consider when filtering, and how much to weight those pixels."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mean_kernel = 1.0/9.0 * np.ones((3, 3))\n",
      "\n",
      "print mean_kernel"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Now, let's take our mean kernel and apply it to every pixel of the image."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Applying a (linear) filter essentially means:\n",
      "* Center a kernel on a pixel\n",
      "* Multiply the pixels *under* that kernel by the values *in* the kernel\n",
      "* Sum all the those results\n",
      "* Replace the center pixel with the summed result"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This process is known as convolution. Let's actually walk through an example of this process below:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import skdemo\n",
      "skdemo.mean_filter_interactive_demo(bright_square)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Let's take a look at the numerical result:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.ndimage import convolve\n",
      "\n",
      "%precision 2\n",
      "print convolve(bright_square, mean_kernel)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The meaning of \"mean kernel\" should be clear now: Each pixel was replaced with the mean value within the 3x3 neighborhood of that pixel. When the kernel was over `n` bright pixels, the pixel in the kernel's center was changed to n/9 (= n * 0.111). When no bright pixels were under the kernel, the result was 0."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This filter is a simple smoothing filter and produces two important results:\n",
      "1. The intensity of the bright pixel decreased.\n",
      "2. The intensity of the region near the bright pixel increased."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Slight aside:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print mean_kernel.sum()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Note that all the values of the kernel sum to 1. Why might that be important? Sure, a definition of a mean requires this property, but why might this be a favorable property for many image filters?"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Downsampled image"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Let's consider a real image now. It'll be easier to see some of the filtering we're doing if we downsample the image a bit. We can slice into the image using the \"step\" argument to sub-sample it:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from skimage import data\n",
      "\n",
      "image = data.camera()\n",
      "pixelated = image[::10, ::10]\n",
      "skdemo.imshow_all(image, pixelated)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Here we use a step of 10, giving us every tenth column and every tenth row of the original image. You can see the highly pixelated result on the right."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Mean filter on a real image"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Now we can apply the filter to this downsampled image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "filtered = convolve(pixelated, mean_kernel)\n",
      "skdemo.imshow_all(pixelated, filtered)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Comparing the filtered image to the pixelated image, we can see that this filtered result is smoother: Sharp edges (which are just borders between dark and bright pixels) are smoothed because dark pixels reduce the intensity of neighboring pixels and bright pixels do the opposite."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Essential filters"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "If you read through the last section, you're already familiar with the essential concepts of image filtering. But, of course, you don't have to create custom filter kernels for all of your filtering needs.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Gaussian filter"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The classic image filter is the Gaussian filter. This is similar to the mean filter, in that it tends to smooth images. The Gaussian filter, however, doesn't weight all values in the neighborhood equally. Instead, pixels closer to the center are weighted more than those farther away."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from skimage import filter\n",
      "\n",
      "sigma = 1\n",
      "smooth = filter.gaussian_filter(bright_square, sigma)\n",
      "skdemo.imshow_all(bright_square, smooth)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "For the Gaussian filter, `sigma`, the standard deviation, defines the size of the neighborhood.\n",
      "\n",
      "For a real image, we get the following:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from skimage import img_as_float\n",
      "# The Gaussian filter returns a float image, regardless of input.\n",
      "# Cast to float so the images have comparable intensity ranges.\n",
      "pixelated_float = img_as_float(pixelated)\n",
      "smooth = filter.gaussian_filter(pixelated_float, 1)\n",
      "skdemo.imshow_all(pixelated_float, smooth)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This doesn't look drastically different than the mean filter, but the Gaussian filter is typically preferred because of the distance-dependent weighting. For a more detailed image and a larger filter, you can see artifacts in the mean filter since it doesn't take distance into account:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "size = 20\n",
      "structuring_element = np.ones((3*size, 3*size))\n",
      "smooth_mean = filter.rank.mean(image, structuring_element)\n",
      "smooth_gaussian = filter.gaussian_filter(image, size)\n",
      "titles = ['mean', 'gaussian']\n",
      "skdemo.imshow_all(smooth_mean, smooth_gaussian, titles=titles)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The size of the structuring element used for the mean filter and the size (standard deviation) of the Gaussian filter are tweaked to produce an approximately equal amount of smoothing in the two results."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Basic edge filtering"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "For images, edges are basically boundaries between light and dark values. The detection of edges can be useful on its own, or it can be used as preliminary step in other algorithms (which we'll see later)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "A 1D difference filter"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Let's first consider a 1D edge, which is just the boundary in a step function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "step = np.zeros(40)\n",
      "step[20:] = 1"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(step, 'k')\n",
      "plt.margins(0.1) "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "To \"detect\" the edge, we can just take a difference of neighboring values. Here, we'll use convolution or cross-correlation to do just that:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "edge = np.convolve(step, np.array([1, -1]), mode='valid')"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.plot(step, 'k:', edge, 'r')\n",
      "plt.margins(0.1)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Whenever neighboring values are equal, the filter response is 0. Right at the boundary of a step, we're subtracting a large value from a small value and and get a spike in the response. This spike \"identifies\" our edge."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Difference filters in 2D"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "For images, you can think of an edge as points where the gradient is large in one direction. We can approximate gradients with difference filters. There are many ways to compute intensity differences between neighboring pixels (by weighting neighbors differently). At its simplest, you can just subtract one neighbor from the other."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "horizontal_edges = pixelated[1:, :] - pixelated[:-1, :]\n",
      "vertical_edges = pixelated[:, 1:] - pixelated[:, :-1]"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This is exactly like the convolution filter we used in 1D, but using slicing operations instead of a filtering kernel. We're avoiding convolution here only because you'll be developing the filtering kernel in an exercise soon enough."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "skdemo.imshow_all(pixelated, horizontal_edges, vertical_edges)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "That's obviously not what we were hoping for: It all looks like noise. What's wrong here?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "In addition to the more obvious issues above, this operation has two additional issues, which can be seen below:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "horizontal_edges = bright_square[1:, :-1] - bright_square[:-1, :-1]\n",
      "vertical_edges = bright_square[:-1, 1:] - bright_square[:-1, :-1]\n",
      "skdemo.imshow_all(bright_square, horizontal_edges, vertical_edges)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print bright_square.shape, horizontal_edges.shape"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Note here that:\n",
      "1. The shape of the image isn't preserved\n",
      "2. The operation skews edges to one corner of the image.\n",
      "\n",
      "This difference operation gives the gradient *in-between* pixels, but we typically want the gradient at the same pixels as the original image."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "<span style=\"color:cornflowerblue\">Exercise:</span>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Create a simple difference filter to **find the horizontal or vertical edges** of an image. Try to ensure that the filtering operation doesn't shift the edge position preferentially. (Don't use slicing to produce the difference image; use convolution.)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This should get you started:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Replace the kernels below with your difference filter\n",
      "# `ones` is used just for demonstration and your kernel should be larger than (1, 1)\n",
      "horizontal_edge_kernel = np.ones((1, 1))\n",
      "vertical_edge_kernel = np.ones((1, 1))\n",
      "\n",
      "# As discussed earlier, you may want to replace pixelated with a different image.\n",
      "image = pixelated\n",
      "horizontal_edges = convolve(image, horizontal_edge_kernel)\n",
      "vertical_edges = convolve(image, vertical_edge_kernel)\n",
      "skdemo.imshow_all(horizontal_edges, vertical_edges)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Sobel edge filter"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The Sobel filter, the most commonly used edge filter, should look pretty similar to what you developed above. Take a look at the vertical and horizontal components of the Sobel kernel to see how they differ from your earlier implementation:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "* http://scikit-image.org/docs/dev/api/skimage.filter.html#vsobel\n",
      "* http://scikit-image.org/docs/dev/api/skimage.filter.html#hsobel"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The standard Sobel filter gives the gradient magnitude. This is similar to what we saw above, except that horizontal and vertical components are combined such that the direction of the gradient is ignored."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "skdemo.imshow_all(bright_square, filter.sobel(bright_square))"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Notice that the size of the output matches the input, and the edges aren't preferentially shifted to a corner of the image. Furthermore, the weights used in the Sobel filter produce diagonal edges with reponses that are comparable to horizontal or vertical edges.\n",
      "\n",
      "Like any derivative, noise can have a strong impact on the result:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pixelated_edges = filter.sobel(pixelated)\n",
      "skdemo.imshow_all(pixelated, pixelated_edges)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Smoothing is often used as a preprocessing step in preparation for feature detection and image-enhancement operations because sharp features can distort results."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "edges = filter.sobel(smooth)\n",
      "titles = ['edges before smoothing', 'edges after smoothing']\n",
      "# Scale smoothed edges up so they're of comparable brightness.\n",
      "skdemo.imshow_all(pixelated_edges, edges*1.8, titles=titles)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Notice how the edges look more continuous in the smoothed image."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "<span style=\"color:cornflowerblue\">Exercise:</span>"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Using the filter kernels from the `vsobel` and `hsobel` documentation, **find the direction of the maximum gradient** in an image. Here's some code to get you started"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Kernel copied from `vsobel` docstring.\n",
      "# Vertical edge-reponse is the *horizontal* gradient.\n",
      "dx_kernel = np.array([\n",
      "    [1,   0,  -1],\n",
      "    [2,   0,  -2],\n",
      "    [1,   0,  -1],\n",
      "])\n",
      "# Rotate array by 90 degrees to get y-gradient kernel\n",
      "dy_kernel = np.rot90(dx_kernel)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "And here are some test images"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "image_45 = np.tril(-np.ones([7, 7]))\n",
      "image_135 = np.rot90(image_45)\n",
      "skdemo.imshow_all(image_45, image_135)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Denoising filters"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This is a bit arbitrary, but here, we distinguish smoothing filters from denoising filters. We'll label denoising filters as those that are edge preserving.\n",
      "\n",
      "As you can see from our earlier examples, mean and Gaussian filters smooth an image rather uniformly, including the edges of objects in an image. When denoising, however, you typically want to preserve features and just remove noise. The distinction between noise and features can, of course, be highly situation-dependent and subjective."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Median Filter"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "The median filter is the classic edge-preserving filter. As the name implies, this filter takes a set of pixels (i.e. the pixels within a kernel or \"structuring element\") and returns the median value within that neighborhood. Because regions near a sharp edge will have many dark values and many light values (but few values in between) the median at an edge will most likely be either light or dark, rather than some value in between. In that way, we don't end up with edges that are smoothed."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from skimage.morphology import disk\n",
      "selem = disk(2)  # \"selem\" is often the name used for \"structuring element\"\n",
      "median = filter.rank.median(pixelated, selem)\n",
      "titles = ['image', 'gaussian', 'median']\n",
      "skdemo.imshow_all(pixelated, smooth, median, titles=titles)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "This difference is more noticeable with a more detailed image."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "selem = disk(10)\n",
      "coins = data.coins()\n",
      "mean_coin = filter.rank.mean(coins, selem)\n",
      "median_coin = filter.rank.median(coins, selem)\n",
      "titles = ['image', 'mean', 'median']\n",
      "skdemo.imshow_all(coins, mean_coin, median_coin, titles=titles)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "Notice how the edges of coins are preserved after using the median filter."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Further reading"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "notes"
      }
     },
     "source": [
      "`scikit-image` also provides more sophisticated denoising filters:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# If this import fails for you, you're probably using an old version\n",
      "# Instead, use: from skimage.filter import denoise_tv_bregman\n",
      "from skimage.restoration import denoise_tv_bregman\n",
      "denoised = denoise_tv_bregman(image, 4)\n",
      "titles = ['image', 'median', 'denoised']\n",
      "skdemo.imshow_all(image, median, denoised, titles=titles)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "* [Denoising examples](http://scikit-image.org/docs/dev/auto_examples/plot_denoise.html)\n",
      "* [Rank filters example](http://scikit-image.org/docs/dev/auto_examples/applications/plot_rank_filters.html)"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}