{"cells":[
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport sklearn\nimport glob, os"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from subprocess import check_output\nprint(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\nsmjpegs = [f for f in glob.glob(\"../input/train_sm/*.jpeg\")]\nprint(smjpegs[:9])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "set175 = [smj for smj in smjpegs if \"set175\" in smj]\nprint(set175)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import math\nfrom PIL import Image\nfor f in sorted(smjpegs):\n  w,h=Image.open(f).size\n  print (f,\" \",w,\" \",h)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "np.min(first)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "print(first)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "pixel_matrix = np.reshape(first, (dims[0] * dims[1], dims[2]))\nprint(np.shape(pixel_matrix))"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "Scatter plots are a go to to look for clusters and separatbility in the data, but these are busy and don't reveal density well, so we\nswitch to using 2d histograms instead. The data between bands is really correlated, typical with\nvisible imagery and why most satellite image analysts prefer to at least have near infrared values."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "#plt.scatter(pixel_matrix[:,0], pixel_matrix[:,1])\n_ = plt.hist2d(pixel_matrix[:,1], pixel_matrix[:,2], bins=(50,50))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fifth = plt.imread('../input/train_sm/set175_5.jpeg')\ndims = np.shape(fifth)\npixel_matrix5 = np.reshape(fifth, (dims[0] * dims[1], dims[2]))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "_ = plt.hist2d(pixel_matrix5[:,1], pixel_matrix5[:,2], bins=(50,50))"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "We can look at variations between the scenes now and see that there's a significant\namount of difference, probably due to sensor angle and illumination variation. Raw band\ndifferences will need to be scaled or thresholded for any traditional approach."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "_ = plt.hist2d(pixel_matrix[:,2], pixel_matrix5[:,2], bins=(50,50))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(first)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(fifth)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "Without coregistering portions of the image, the naive red band subtraction for change indication\nbasically just shows the location shift between images."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(first[:,:,2] - fifth[:,:,1])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "second = plt.imread('../input/train_sm/set175_2.jpeg')\nplt.imshow(first[:,:,2] - second[:,:,2])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(second)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Initial impressions\n\nImages aren't registered, so an image registration process between images with common overlap would probably be the first step in a traditional approach.\nUsing a localizer in a deep learning context would probably be the newfangled way to tackle this.\n\nImage content and differences will be dominated by topographic and built variations\ndue to sensor orientation, resolution differences between scenes, and some registration accuracy will be impossible to factor out as\nthe image hasn't been orthorectified and some anciliary data would be required for it\nto be done, e.g. georeferenceing against a previously orthorectified image.\n\nSo this is basically a basic computer vision task that deep learning will be a good fit for. The usual preprocessing steps\nand data expectations you'd see in remote sensing aren't fulfilled by this dataset."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# simple k means clustering\nfrom sklearn import cluster\n\nkmeans = cluster.KMeans(5)\nclustered = kmeans.fit_predict(pixel_matrix)\n\ndims = np.shape(first)\nclustered_img = np.reshape(clustered, (dims[0], dims[1]))\nplt.imshow(clustered_img)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(first)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "ind0, ind1, ind2, ind3 = [np.where(clustered == x)[0] for x in [0, 1, 2, 3]]"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "This code doesn't run on the server.\n\n```python\nfrom mpl_toolkits.mplot3d import Axes3D\n\nfig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\n\nplot_vals = [('r', 'o', ind0),\n             ('b', '^', ind1),\n             ('g', '8', ind2),\n             ('m', '*', ind3)]\n\nfor c, m, ind in plot_vals:\n    xs = pixel_matrix[ind, 0]\n    ys = pixel_matrix[ind, 1]\n    zs = pixel_matrix[ind, 2]\n    ax.scatter(xs, ys, zs, c=c, marker=m)\n\nax.set_xlabel('Blue channel')\nax.set_ylabel('green channel')\nax.set_zlabel('Red channel')\n```"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# quick look at color value histograms for pixel matrix from first image\nimport seaborn as sns\nsns.distplot(pixel_matrix[:,0], bins=12)\nsns.distplot(pixel_matrix[:,1], bins=12)\nsns.distplot(pixel_matrix[:,2], bins=12)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# even subsampling is throwing memory error for me, :p\n#length = np.shape(pixel_matrix)[0]\n#rand_ind = np.random.choice(length, size=50000)\n#sns.pairplot(pixel_matrix[rand_ind,:])"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Day 2\n\nWe'll start by considering the entire sequence of a different image set this time and look at strategies\nfor matching features across scenes."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "set79 = [smj for smj in smjpegs if \"set79\" in smj]\nprint(set79)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img79_1, img79_2, img79_3, img79_4, img79_5 = \\\n  [plt.imread(\"../input/train_sm/set79_\" + str(n) + \".jpeg\") for n in range(1, 6)]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img_list = (img79_1, img79_2, img79_3, img79_4, img79_5)\n\nplt.figure(figsize=(8,10))\nplt.imshow(img_list[0])\nplt.show()"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "Tracking dimensions across image transforms is annoying, so we'll make a class to do that.\nAlso I'm going to use this brightness normalization transform and visualize the image that\nway, good test scenario for class."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "class MSImage():\n    \"\"\"Lightweight wrapper for handling image to matrix transforms. No setters,\n    main point of class is to remember image dimensions despite transforms.\"\"\"\n    \n    def __init__(self, img):\n        \"\"\"Assume color channel interleave that holds true for this set.\"\"\"\n        self.img = img\n        self.dims = np.shape(img)\n        self.mat = np.reshape(img, (self.dims[0] * self.dims[1], self.dims[2]))\n\n    @property\n    def matrix(self):\n        return self.mat\n        \n    @property\n    def image(self):\n        return self.img\n    \n    def to_flat_img(self, derived):\n        \"\"\"\"Use dims property to reshape a derived matrix back into image form when\n        derived image would only have one band.\"\"\"\n        return np.reshape(derived, (self.dims[0], self.dims[1]))\n    \n    def to_matched_img(self, derived):\n        \"\"\"\"Use dims property to reshape a derived matrix back into image form.\"\"\"\n        return np.reshape(derived, (self.dims[0], self.dims[1], self.dims[2]))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "msi79_1 = MSImage(img79_1)\nprint(np.shape(msi79_1.matrix))\nprint(np.shape(msi79_1.img))"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Brightness Normalization\n\nBrightness Normalization is preprocessing strategy you can apply prior to using strategies\nto identify materials in a scene, if you want your matching algorithm\nto be robust across variations in illumination. See [Wu's paper](https://pantherfile.uwm.edu/cswu/www/my%20publications/2004_RSE.pdf)."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "def bnormalize(mat):\n    \"\"\"much faster brightness normalization, since it's all vectorized\"\"\"\n    bnorm = np.zeros_like(mat, dtype=np.float32)\n    maxes = np.max(mat, axis=1)\n    bnorm = mat / np.vstack((maxes, maxes, maxes)).T\n    return bnorm"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "bnorm = bnormalize(msi79_1.matrix)\nbnorm_img = msi79_1.to_matched_img(bnorm)\nplt.figure(figsize=(8,10))\nplt.imshow(bnorm_img)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "msi79_2 = MSImage(img79_2)\nbnorm79_2 = bnormalize(msi79_2.matrix)\nbnorm79_2_img = msi79_2.to_matched_img(bnorm79_2)\nplt.figure(figsize=(8,10))\nplt.imshow(bnorm79_2_img)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "msinorm79_1 = MSImage(bnorm_img)\nmsinorm79_2 = MSImage(bnorm79_2_img)\n\n_ = plt.hist2d(msinorm79_1.matrix[:,2], msinorm79_2.matrix[:,2], bins=(50,50))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "_ = plt.hist2d(msinorm79_1.matrix[:,1], msinorm79_2.matrix[:,1], bins=(50,50))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "_ = plt.hist2d(msinorm79_1.matrix[:,0], msinorm79_2.matrix[:,0], bins=(50,50))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import seaborn as sns\nsns.distplot(msinorm79_1.matrix[:,0], bins=12)\nsns.distplot(msinorm79_1.matrix[:,1], bins=12)\nsns.distplot(msinorm79_1.matrix[:,2], bins=12)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(8,10))\nplt.imshow(img79_1)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "np.max(img79_1[:,:,0])"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Using thresholds with brightness normalization\n\nOk, so what am I even doing here? Well, my goal is to try and figure out simple threshold selection\nmethods for getting high albedo targets out of a scene so I could then theoretically track them\nbetween scenes. For example, a simple blob/aggregation to centroid (in coordinates or in subsampled\nimage bins) would give me a means to look at plausible structural similarities in distributions\nbetween scenes, then use that to anchor a comparison of things that change.\n\nThe brightness normalization step is helpful because thresholds that aren't anchored by a\npreprocessing step end up being arbitrary and can't generalize between scenes even in the same\nimage set, whereas thresholds following brightness normalization tend to pull out materils that stand\nout from the background more reliably. See the following demonstration:"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(img79_1[:,:,0] > 230)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(img79_2[:,:,0] > 230)\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "print(np.min(bnorm79_2_img[:,:,0]))\nprint(np.max(bnorm79_2_img[:,:,0]))\nprint(np.mean(bnorm79_2_img[:,:,0]))\nprint(np.std(bnorm79_2_img[:,:,0]))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm79_2_img[:,:,0] > 0.98)\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm_img[:,:,0] > 0.98)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow((bnorm79_2_img[:,:,0] > 0.9999) & \\\n           (bnorm79_2_img[:,:,1] < 0.9999) & \\\n           (bnorm79_2_img[:,:,2] < 0.9999))\nplt.subplot(122)\nplt.imshow(img79_2)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(121)\nplt.imshow(bnorm_img[:,:,0] > 0.995)\nplt.subplot(122)\nplt.imshow(img79_1)\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,6))\nplt.subplot(121)\nplt.plot(bnorm_img[2000, 1000, :])\nplt.subplot(122)\nplt.plot(img79_1[2000, 1000, :])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from scipy import spatial\n\npixel = msi79_1.matrix[2000 * 1000, :]\nnp.shape(pixel)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Something's borked here\n\nThink I'm gonna have to verify cosine similarity behavior for scipy here.\n\n```python\ndef spectral_angle_mapper(pixel):\n    return lambda p2: spatial.distance.cosine(pixel, p2)\n\nmatch_pixel = np.apply_along_axis(spectral_angle_mapper(pixel), 1, msi79_1.matrix)\n\nplt.figure(figsize=(10,6))\nplt.imshow(msi79_1.to_flat_img(match_pixel < 0.0000001))\n\ndef summary(mat):\n    print(\"Max: \", np.max(mat),\n          \"Min: \", np.min(mat),\n          \"Std: \", np.std(mat),\n          \"Mean: \", np.mean(mat))\n    \nsummary(match_pixel)\n```"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Rudimentary Transforms, Edge Detection, Texture"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "set144 = [MSImage(plt.imread(smj)) for smj in smjpegs if \"set144\" in smj]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(set144[0].image)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import skimage\nfrom skimage.feature import greycomatrix, greycoprops\nfrom skimage.filters import sobel"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Sobel Edge Detection\n\nA Sobel filter is one means of getting a basic edge magnitude/gradient image. Can be useful to\nthreshold and find prominent linear features, etc. Several other similar filters in skimage.filters\nare also good edge detectors: `roberts`, `scharr`, etc. and you can control direction, i.e. use\nan anisotropic version."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# a sobel filter is a basic way to get an edge magnitude/gradient image\nfig = plt.figure(figsize=(8, 8))\nplt.imshow(sobel(set144[0].image[:750,:750,2]))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from skimage.filters import sobel_h\n\n# can also apply sobel only across one direction.\nfig = plt.figure(figsize=(8, 8))\nplt.imshow(sobel_h(set144[0].image[:750,:750,2]), cmap='BuGn')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from sklearn.decomposition import PCA\n\npca = PCA(3)\npca.fit(set144[0].matrix)\nset144_0_pca = pca.transform(set144[0].matrix)\nset144_0_pca_img = set144[0].to_matched_img(set144_0_pca)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.imshow(set144_0_pca_img[:,:,0], cmap='BuGn')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.imshow(set144_0_pca_img[:,:,1], cmap='BuGn')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.imshow(set144_0_pca_img[:,:,2], cmap='BuGn')"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# GLCM Textures\n\nProcessing time can be pretty brutal so we subset the image. We'll create texture images so\nwe can characterize each pixel by the texture of its neighborhood.\n\nGLCM is inherently anisotropic but can be averaged so as to be rotation invariant. For more on GLCM, see [the tutorial](http://www.fp.ucalgary.ca/mhallbey/tutorial.htm).\n\nA good article on use in remote sensing is [here](http://ieeexplore.ieee.org/xpl/login.jsp?tp=&arnumber=4660321&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%3Farnumber%3D4660321):\n\nPesaresi, M., Gerhardinger, A., & Kayitakire, F. (2008). A robust built-up area presence index by anisotropic rotation-invariant textural measure. Selected Topics in Applied Earth Observations and Remote Sensing, IEEE Journal of, 1(3), 180-192."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "sub = set144[0].image[:150,:150,2]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "def glcm_image(img, measure=\"dissimilarity\"):\n    \"\"\"TODO: allow different window sizes by parameterizing 3, 4. Also should\n    parameterize direction vector [1] [0]\"\"\"\n    texture = np.zeros_like(sub)\n\n    # quadratic looping in python w/o vectorized routine, yuck!\n    for i in range(img.shape[0] ):  \n        for j in range(sub.shape[1] ):  \n          \n            # don't calculate at edges\n            if (i < 3) or \\\n               (i > (img.shape[0])) or \\\n               (j < 3) or \\\n               (j > (img.shape[0] - 4)):          \n                continue  \n        \n            # calculate glcm matrix for 7 x 7 window, use dissimilarity (can swap in\n            # contrast, etc.)\n            glcm_window = img[i-3: i+4, j-3 : j+4]  \n            glcm = greycomatrix(glcm_window, [1], [0],  symmetric = True, normed = True )   \n            texture[i,j] = greycoprops(glcm, measure)  \n    return texture"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "dissimilarity = glcm_image(sub, \"dissimilarity\")"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.subplot(1,2,1)\nplt.imshow(dissimilarity, cmap=\"bone\")\nplt.subplot(1,2,2)\nplt.imshow(sub, cmap=\"bone\")"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# HSV Transform\n\nSince this contest is about time series ordering, I think it's possible there may be useful\ninformation in a transform to HSV color space. HSV is useful for identifying shadows and illumination, as well\nas giving us a means to identify similar objects that are distinct by color between scenes (hue), \nthough there's no guarantee the hue will be stable."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from skimage import color\n\nhsv = color.rgb2hsv(set144[0].image)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.subplot(2,2,1)\nplt.imshow(set144[0].image, cmap=\"bone\")\nplt.subplot(2,2,2)\nplt.imshow(hsv[:,:,0], cmap=\"bone\")\nplt.subplot(2,2,3)\nplt.imshow(hsv[:,:,1], cmap='bone')\nplt.subplot(2,2,4)\nplt.imshow(hsv[:,:,2], cmap='bone')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.subplot(2,2,1)\nplt.imshow(set144[0].image[:200,:200,:])\nplt.subplot(2,2,2)\nplt.imshow(hsv[:200,:200,0], cmap=\"PuBuGn\")\nplt.subplot(2,2,3)\nplt.imshow(hsv[:200,:200,1], cmap='bone')\nplt.subplot(2,2,4)\nplt.imshow(hsv[:200,:200,2], cmap='bone')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 6))\nplt.imshow(hsv[200:500,200:500,0], cmap='bone')"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "hsvmsi = MSImage(hsv)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Shadow Detection\n\nWe can apply a threshold to the V band now to find dark areas that are probably thresholds. Let's\nlook at the distribution of all values then work interactively to find a good filter value."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import seaborn as sns\nsns.distplot(hsvmsi.matrix[:,0], bins=12)\nsns.distplot(hsvmsi.matrix[:,1], bins=12)\nsns.distplot(hsvmsi.matrix[:,2], bins=12)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(hsvmsi.image[:,:,2] < 0.4, cmap=\"plasma\")"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nplt.subplot(1,2,1)\nplt.imshow(set144[0].image[:250,:250,:])\nplt.subplot(1,2,2)\nplt.imshow(hsvmsi.image[:250,:250,2] < 0.4, cmap=\"plasma\")"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "fig = plt.figure(figsize=(8, 8))\nimg2 = plt.imshow(set144[0].image[:250,:250,:], interpolation='nearest')\nimg3 = plt.imshow(hsvmsi.image[:250,:250,2] < 0.4, cmap='binary_r', alpha=0.4)\nplt.show()"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "Could we glean something useful about sun position from shadow orientation if we could accurately\nreference the image?"
 }
],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}}, "nbformat": 4, "nbformat_minor": 0}