Working with images in python/matplotlib

Working with Images.ipynb

{
 "metadata": {
  "name": "Working with Images"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Working with Images"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "attempt at Notebook transcript of Py4Science meeting, during interactive demos\n",
      "\n",
      "Includes working with image data, and automatic display of PIL Image objects."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imgpath = matplotlib.cbook.get_sample_data('stinkbug.png', False)\n",
      "print \"loading image from: %s\" % imgpath\n",
      "img = plt.imread(imgpath)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Information about the image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "type(img)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img.ndim"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img.dtype"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's view the image:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(img)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "I want to pull out just one color, how would I throw away the other two?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imgg = img[:,:,0]\n",
      "print \"shape:\", imgg.shape\n",
      "imshow(imgg)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Why does it have those colors?\n",
      "\n",
      "Because imshow(2d array) shows a colormap."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(imgg)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we have a color image that appears grey, and a grey image that has color!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Images can have 1 (bw), 3 (rgb), or 4 (rgba) colors, but 2 is invalid:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(img[:,:,:2])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can set the color map with `implot_.set_cmap(name)`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "implot = imshow(imgg)\n",
      "implot.set_cmap('hot')\n",
      "display(implot.get_figure())\n",
      "implot.set_cmap('gray')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "matplotlib has a `plt.hist` function for doing a histogram"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(imgg);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But that doesn't look right for a luminosity histogram.\n",
      "If we *flatten* the image, the our luminosity will look more sensible, \n",
      "because the data will be 1D, and use just 256 bins for the histogram."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.hist(imgg.flat, 256, range=(0,1));"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now let's have a look at dynamic range"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "matplotlib has a very useful little module called `matlotlib.cm` for all the color maps.  \n",
      "So you can tab-complete on them, and see all the (many) options.\n",
      "\n",
      "You can also set the cmap by name with the `cmap` kwarg (we called `set_cmap()` previously)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import cm\n",
      "imshow(imgg, cmap=cm.hot)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The colorbar reveals that we don't actually have any data in the top 15% of the range (1.0 is white).\n",
      "\n",
      "We can clip the color bar at any point, for instance, let's cut off at 70%:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p2 = imshow(imgg, cmap=cm.hot)\n",
      "colorbar()\n",
      "p2.set_clim(0, 0.7)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now the highlights stand out more (they are white, because they are off scale).\n",
      "\n",
      "Setting the color limits are a common operation, \n",
      "often setting the min and max for each plot to a *global* min and max, across all your data.\n",
      "That way, all your plots are on the same scale, and are more easily comparable."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Now the PIL\n",
      "\n",
      "The [Python Imaging Library](http://pythonware.com/products/pil/)\n",
      "provides more traditional *image* manipulation.  It imports as `Image`.\n",
      "\n",
      "The PIL has very good support for a vast array of image formats, and provides\n",
      "standard 'filter-style' image manipulation operations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "import Image"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "we open an image with `Image.open(fname)`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img = Image.open(imgpath)\n",
      "img"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can register PIL Images for automatic display:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from io import BytesIO\n",
      "\n",
      "def display_image(img, format='PNG'):\n",
      "    fp = BytesIO()\n",
      "    img.save(fp, format=format)\n",
      "    return fp.getvalue()\n",
      "\n",
      "ip = get_ipython()\n",
      "png_formatter = ip.display_formatter.formatters['image/png']\n",
      "png_formatter.for_type_by_name('Image', 'Image', display_image)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can do common image operations, like resize:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "rsize = img.resize((img.size[0]/10, img.size[1]/10))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "we can load the new small image into an array:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "rsize_arr = np.asarray(rsize)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print img.size\n",
      "print rsize_arr.shape\n",
      "print rsize_arr.dtype"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can use `np.asarray(img)` because PIL Image objects provide something called the Python buffer interface,\n",
      "which lets various different containers to talk to each other."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(rsize_arr, interpolation='nearest')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "nearest-neighbor interpolation just draws big square pixels, since the output is larger than the input.\n",
      "\n",
      "But if we use bicubic interpolation, it will look smooth (and blurry):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(rsize_arr, interpolation='bicubic')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To get an Image back from an array, you can use `Image.fromarray(arr)`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Image.fromarray(rsize_arr)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Moon Landing and denoising\n",
      "\n",
      "We have an image with lots of high frequency noise:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imgpath = matplotlib.cbook.get_sample_data('moonlanding.png', False)\n",
      "print \"loading image from: %s\" % imgpath\n",
      "img = imread(imgpath)\n",
      "imshow(img, cmap=cm.Greys_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So let's try a super-simple band-pass, simply chopping off the high-frequency elements.\n",
      "It might be blurry, but we might be able to recover some low-frequency information.\n",
      "\n",
      "For this, we need to take [Fast Fourier Transform](http://en.wikipedia.org/wiki/Fft) of the image:\n",
      "\n",
      "$X_k =  \\sum_{n=0}^{N-1} x_n e^{-{i 2\\pi k \\frac{n}{N}}}\n",
      "\\qquad\n",
      "k = 0,\\dots,N-1.$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First thing, we need to understand how numpy organizes the frequencies in the FFT.\n",
      "\n",
      "In a 2d FFT, numpy puts the high frequency data in the *center* of the image array,\n",
      "so the low-frequency data is focused at the corners of the array."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "F = np.fft.fft2(img)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "F.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "F.dtype"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "it's complex!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(abs(F))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "weird, and kind of useless.  Let's put a colorbar on it."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(abs(F))\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Wow, dynamic range is an issue, because we have a few pixels dominating everything).\n",
      "\n",
      "Let's try a log to smooth it out:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(log(abs(F)))\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That's better, but let's try a cutoff, like we did before, but this time\n",
      "cut off everything but the very bottom:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax = imshow(abs(F))\n",
      "ax.set_clim(None, 1500)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Cut off *even lower*:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax = imshow(abs(F))\n",
      "ax.set_clim(None, 300)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "let's look at the histogram"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(abs(F).flat, 256);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "see, there's nothing anywhere but the first bin.\n",
      "So what if we only look at the histogram of the most recent"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(abs(F).flat, 256, range=(0,200));"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The high frequency information is in the *center* of the FFT array.\n",
      "If we want to eliminate high frequency noise, we can simply discard (set to zero)\n",
      "the high frequency region.\n",
      "\n",
      "So how do we do that?"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's copy our FFT, so we don't destroy the array, and define the fraction we want to keep:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "Fc = F.copy()\n",
      "keep = 0.1 # 10%"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Get our number of rows and number of columns:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "nrow, ncol = img.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now zero out everything but the *outer* band, defined by `keep`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "# zero out the rows\n",
      "Fc[keep*nrow : (1-keep)*nrow, :] = 0\n",
      "# zero out the columns\n",
      "Fc[:, keep*ncol : (1-keep)*ncol] = 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's look at what we have done:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax = imshow(abs(Fc))\n",
      "ax.set_clim(0,200)\n",
      "colorbar()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's look at the resulting image, with almost all of the frequency information thrown away.\n",
      "\n",
      "(just take the real part, because we know we have a real signal, but will get some imaginary part,\n",
      "due to "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img_new = np.fft.ifft2(Fc).real\n",
      "imshow(img_new, cmap=cm.Greys_r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not bad for fingerpainting with Fourier Transforms."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}