dobler_111714_part2.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:fbe519096d56c59d1cf874578304b6093e89c2a18b00ddb991c9a4178707656f"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(always be aware of your imports and <b><u><i>preserve namespaces</i></u></b>!!!)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import os\n",
      "import sys\n",
      "import time\n",
      "import numpy as np\n",
      "import scipy.ndimage as nd\n",
      "import matplotlib.pyplot as plt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---\n",
      "\n",
      "## Advanced (and relational) indexing: Part 2\n",
      "\n",
      "In addition to indexing based on position within an image, we can index on pixel <i>values</i> which is useful for selecting subcomponents of an image that are not conected compoenents.  As an example, let's create an image consisting of a blue circle and line:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rind = np.arange(20000).reshape(200,100) / 100\n",
      "cind = np.arange(20000).reshape(200,100) % 100\n",
      "two_ln = ((rind>130) & (rind<150)) | \\\n",
      "    (np.sqrt((rind-50)**2+(cind-50)**2)<20)\n",
      "two_ln_c = np.dstack([64*two_ln, 128*two_ln, 192*two_ln]\n",
      "                     ).astype(np.uint8)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Plotting it:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig0, ax0 = plt.subplots(num=0,figsize=[5,10])\n",
      "fig0.subplots_adjust(0,0,1,1)\n",
      "ax0.grid('off')\n",
      "im0 = ax0.imshow(two_ln_c)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now identify the pixels in 2D which have the specified color and set those pixels to brown,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "two_ln_c.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "(200, 100, 3)"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "index = (two_ln_c[:,:,0]==64) & \\\n",
      "    (two_ln_c[:,:,1]==128) & \\\n",
      "    (two_ln_c[:,:,2]==192)\n",
      "two_ln_c[index] = [192,128,64]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "and replot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "im0.set_data(two_ln_c)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.close(0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's try it on an actual image."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dpath = 'images'\n",
      "fname = 'dalmatian.jpg'\n",
      "infile = os.path.join(dpath,fname)\n",
      "img_dog = nd.imread(infile)\n",
      "\n",
      "nrow, ncol = img_dog.shape[:2]\n",
      "xsize = 8.\n",
      "ysize = xsize*float(nrow)/float(ncol)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig0, ax0 = plt.subplots(num=0,figsize=[xsize,ysize])\n",
      "fig0.subplots_adjust(0,0,1,1)\n",
      "ax0.grid('off')\n",
      "im0 = ax0.imshow(img_dog)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### TODO: Turn the background black and turn the white fur red."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ind_bkg = (img_dog[:,:,0]==255) & (img_dog[:,:,1]==255) & \\\n",
      "    (img_dog[:,:,2]==255)\n",
      "im0.set_data(ind_bkg)\n",
      "im0.set_clim([0,1])\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "img_dog[ind_bkg] = [0,0,0]\n",
      "im0.set_data(img_dog)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ind_fur = (img_dog[:,:,0]>150) & (img_dog[:,:,1]>150) & \\\n",
      "    (img_dog[:,:,2]>150)\n",
      "img_dog[ind_fur] = [255,128,128]\n",
      "im0.set_data(img_dog)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.close(0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---\n",
      "\n",
      "## Example: Invasive Species"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Freshkills is a 2,200 acre former landfill on Staten Island that is in the process of being converted into a park (Freshkills park).  There is significant contamination by an invasive species called phragmites.\n",
      "\n",
      "#### CUSP PhD student Nick Johnson flew a balloon over Freshkills Park to image the extent of the phragmites invasion."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dpath = 'images'\n",
      "dfile = 'phrag-annotated.png'\n",
      "infile = os.path.join(dpath,dfile)\n",
      "\n",
      "img = nd.imread(infile)[:,:,:3] # ignore the alpha channel\n",
      "\n",
      "nrow, ncol = img.shape[:2]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "xsize = 7.\n",
      "ysize = xsize*float(nrow)/float(ncol)\n",
      "\n",
      "fig0, ax0 = plt.subplots(num=0,figsize=[xsize,ysize])\n",
      "fig0.subplots_adjust(0,0,1,1); ax0.axis('off')\n",
      "im0 = ax0.imshow(img)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What separates the phragmites from the healthy vegetation is the <i>color</i>.  But we already know how to index on color, so let's get an idea of what color we're looking for.  If we zoom in on a region dominated by phragmites:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "prow = [1460,1605]\n",
      "pcol = [4015,4354]\n",
      "\n",
      "ax0.add_patch(plt.Rectangle((pcol[0],prow[0]),\n",
      "                           pcol[1]-pcol[0],prow[1]-prow[0],\n",
      "              facecolor='none',edgecolor='orange',lw=2))\n",
      "ax0.set_xlim(pcol[0]-100,pcol[1]+100)\n",
      "ax0.set_ylim(prow[0]-100,prow[1]+100)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can see that it's a bit bluer than a region with healthy vegetation:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hrow = [2650,2905]\n",
      "hcol = [570,915]\n",
      "\n",
      "ax0.add_patch(plt.Rectangle((hcol[0],hrow[0]),\n",
      "                           hcol[1]-hcol[0],hrow[1]-hrow[0],\n",
      "              facecolor='none',edgecolor='orange',lw=2))\n",
      "ax0.set_xlim(hcol[0]-100,hcol[1]+100)\n",
      "ax0.set_ylim(hrow[0]-100,hrow[1]+100)\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ax0.set_xlim(0,ncol); ax0.set_ylim(nrow,0); fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To start understanding what kind of indexing we need to do, let's make some histograms."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "stamp_ph = img[prow[0]:prow[1],pcol[0]:pcol[1]]\n",
      "stamp_he = img[hrow[0]:hrow[1],hcol[0]:hcol[1]]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig1, ax1 = plt.subplots(2,2,num=1,figsize=[2*xsize,xsize])\n",
      "\n",
      "[ax1[0,0].hist(stamp_ph[:,:,i].flatten(),bins=128,range=[0,255], \n",
      "               normed=True) for i in range(3)]\n",
      "[ax1[1,0].hist(stamp_he[:,:,i].flatten(),bins=128,range=[0,255], \n",
      "               normed=True) for i in range(3)]\n",
      "[i.set_yticklabels('') for j in ax1 for i in j]\n",
      "\n",
      "ax1[0,0].set_ylabel('phragmites')\n",
      "ax1[1,0].set_ylabel('healthy')\n",
      "fig1.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So the first thing we can do is cut on the ratio of green to blue since the healthy vegetation has a ratio >1.5 for (more or less) all pixels in our patch."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hind = (1.0*img[:,:,1])/(1.0*img[:,:,2]) < 1.8"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "-c:1: RuntimeWarning: divide by zero encountered in divide\n",
        "-c:1: RuntimeWarning: invalid value encountered in divide\n",
        "-c:1: RuntimeWarning: invalid value encountered in less\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hind = (1.0*img[:,:,1])/(1.0*img[:,:,2] + (img[:,:,2]==0)) < 1.5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.close(2)\n",
      "fig2, ax2 = plt.subplots(num=2,figsize=[xsize,ysize])\n",
      "fig2.subplots_adjust(0,0,1,1); ax2.axis('off')\n",
      "im2 = ax2.imshow(hind)\n",
      "fig2.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Looks like we're getting phragmites, dirt, and road.  Let's look at the histograms of the dirt and road:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "drow = [2450,2880]\n",
      "dcol = [1905,2235]\n",
      "rrow = [1665,1790]\n",
      "rcol = [3600,3725]\n",
      "\n",
      "stamp_di = img[drow[0]:drow[1],dcol[0]:dcol[1]]\n",
      "stamp_rd = img[rrow[0]:rrow[1],rcol[0]:rcol[1]]\n",
      "\n",
      "for reg in [[drow,dcol],[rrow,rcol]]:\n",
      "    ax0.add_patch(plt.Rectangle((reg[1][0],reg[0][0]),\n",
      "                               reg[1][1]-reg[1][0],reg[0][1]-reg[0][0],\n",
      "                  facecolor='none',edgecolor='orange',lw=2))\n",
      "fig0.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "[ax1[0,1].hist(stamp_di[:,:,i].flatten(),bins=128,range=[0,255], \n",
      "               normed=True) for i in range(3)]\n",
      "[ax1[1,1].hist(stamp_rd[:,:,i].flatten(),bins=128,range=[0,255], \n",
      "               normed=True) for i in range(3)]\n",
      "\n",
      "ax1[0,1].set_ylabel('dirt')\n",
      "ax1[1,1].set_ylabel('road')\n",
      "\n",
      "fig1.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "It looks like \"dirt\" has a high red-to-blue ratio compared to phragmites, so we update our definition of phragmites colors:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hind = ((1.0*img[:,:,1])/(1.0*img[:,:,2] + (img[:,:,2]==0)) < 1.5) & \\\n",
      "    ((1.0*img[:,:,0])/(1.0*img[:,:,2] + (img[:,:,2]==0)) < 1.2)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "im2.set_data(hind)\n",
      "fig2.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Lastly, the road mostly has values $>200$ while the phragmites are $<200$ (also let's use the fact that phragmites have red values which are $<150$ and the lines are black)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hind = ((1.0*img[:,:,1])/(1.0*img[:,:,2] + (img[:,:,2]==0)) < 1.5) & \\\n",
      "    ((1.0*img[:,:,0])/(1.0*img[:,:,2] + (img[:,:,2]==0)) < 1.2) & \\\n",
      "    (img.max(2)<200) & \\\n",
      "    (img[:,:,0]<150) & \\\n",
      "    (img.max(2)>50)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "im2.set_data(hind)\n",
      "fig2.canvas.draw()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}