{"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": "markdown",
  "metadata": {},
  "source": "# Read in files\n\nThis is pretty routine stuff.\n\n* We get a list of jpeg files, reading them in as needed with `matplotlib.pyplot.imread`."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": ""
 },
 {
  "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": "markdown",
  "metadata": {},
  "source": "# Basic exploration\n\nJust look at image dimensions, confirm it's 3 band (RGB), byte scaled (0-255)."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": ""
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "first = plt.imread('../input/train_sm/set175_1.jpeg')\ndims = np.shape(first)\nprint(dims)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "np.min(first), np.max(first)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "For any image specific classification, clustering, etc. transforms we'll want to \ncollapse spatial dimensions so that we have a matrix of pixels by color channels."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": ""
 },
 {
  "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": ""
 },
 {
  "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": ""
 },
 {
  "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": ""
 },
 {
  "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": ""
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": ""
 },
 {
  "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": ""
 },
 {
  "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?"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Image Registration\n\nThis is an earlier form the library found [here](https://github.com/matejak/imreg_dft).\n\nBSD family license, reproduced below with copyright so I can utilize similar functions here where\nimport isn't available.\n\nThis version can be found [here](http://www.lfd.uci.edu/~gohlke/code/imreg.py.html)."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# -*- coding: utf-8 -*-\n# imreg.py\n\n# Copyright (c) 2011-2014, Christoph Gohlke\n# Copyright (c) 2011-2014, The Regents of the University of California\n# Produced at the Laboratory for Fluorescence Dynamics\n# All rights reserved.\n#\n# Redistribution and use in source and binary forms, with or without\n# modification, are permitted provided that the following conditions are met:\n#\n# * Redistributions of source code must retain the above copyright\n#   notice, this list of conditions and the following disclaimer.\n# * Redistributions in binary form must reproduce the above copyright\n#   notice, this list of conditions and the following disclaimer in the\n#   documentation and/or other materials provided with the distribution.\n# * Neither the name of the copyright holders nor the names of any\n#   contributors may be used to endorse or promote products derived\n#   from this software without specific prior written permission.\n#\n# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\"\n# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE\n# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE\n# ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE\n# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR\n# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF\n# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS\n# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN\n# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)\n# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE\n# POSSIBILITY OF SUCH DAMAGE.\n\n\"\"\"FFT based image registration.\n\nImplements an FFT-based technique for translation, rotation and scale-invariant\nimage registration [1].\n\n:Author:\n  `Christoph Gohlke <http://www.lfd.uci.edu/~gohlke/>`_\n\n:Organization:\n  Laboratory for Fluorescence Dynamics, University of California, Irvine\n\n:Version: 2013.01.18\n\nRequirements\n------------\n* `CPython 2.7 or 3.3 <http://www.python.org>`_\n* `Numpy 1.7 <http://www.numpy.org>`_\n* `Scipy 0.12 <http://www.scipy.org>`_\n* `Matplotlib 1.2 <http://www.matplotlib.org>`_  (optional for plotting)\n\nNotes\n-----\nThe API and algorithms are not stable yet and are expected to change between\nrevisions.\n\nReferences\n----------\n(1) An FFT-based technique for translation, rotation and scale-invariant\n    image registration. BS Reddy, BN Chatterji.\n    IEEE Transactions on Image Processing, 5, 1266-1271, 1996\n(2) An IDL/ENVI implementation of the FFT-based algorithm for automatic\n    image registration. H Xiea, N Hicksa, GR Kellera, H Huangb, V Kreinovich.\n    Computers & Geosciences, 29, 1045-1055, 2003.\n(3) Image Registration Using Adaptive Polar Transform. R Matungka, YF Zheng,\n    RL Ewing. IEEE Transactions on Image Processing, 18(10), 2009.\n\nExamples\n--------\n>>> im0 = imread('t400')\n>>> im1 = imread('Tr19s1.3')\n>>> im2, scale, angle, (t0, t1) = similarity(im0, im1)\n>>> imshow(im0, im1, im2)\n\n>>> im0 = imread('t350380ori')\n>>> im1 = imread('t350380shf')\n>>> t0, t1 = translation(im0, im1)\n\n\"\"\"\n\nfrom __future__ import division, print_function\n\nimport math\n\nimport numpy\nfrom numpy.fft import fft2, ifft2, fftshift\n\ntry:\n    import scipy.ndimage.interpolation as ndii\nexcept ImportError:\n    import ndimage.interpolation as ndii\n\n__version__ = '2013.01.18'\n__docformat__ = 'restructuredtext en'\n__all__ = ['translation', 'similarity']\n\n\ndef translation(im0, im1):\n    \"\"\"Return translation vector to register images.\"\"\"\n    shape = im0.shape\n    f0 = fft2(im0)\n    f1 = fft2(im1)\n    ir = abs(ifft2((f0 * f1.conjugate()) / (abs(f0) * abs(f1))))\n    t0, t1 = numpy.unravel_index(numpy.argmax(ir), shape)\n    if t0 > shape[0] // 2:\n        t0 -= shape[0]\n    if t1 > shape[1] // 2:\n        t1 -= shape[1]\n    return [t0, t1]\n\n\ndef similarity(im0, im1):\n    \"\"\"Return similarity transformed image im1 and transformation parameters.\n\n    Transformation parameters are: isotropic scale factor, rotation angle (in\n    degrees), and translation vector.\n\n    A similarity transformation is an affine transformation with isotropic\n    scale and without shear.\n\n    Limitations:\n    Image shapes must be equal and square.\n    All image areas must have same scale, rotation, and shift.\n    Scale change must be less than 1.8.\n    No subpixel precision.\n\n    \"\"\"\n    if im0.shape != im1.shape:\n        raise ValueError(\"Images must have same shapes.\")\n    elif len(im0.shape) != 2:\n        raise ValueError(\"Images must be 2 dimensional.\")\n\n    f0 = fftshift(abs(fft2(im0)))\n    f1 = fftshift(abs(fft2(im1)))\n\n    h = highpass(f0.shape)\n    f0 *= h\n    f1 *= h\n    del h\n\n    f0, log_base = logpolar(f0)\n    f1, log_base = logpolar(f1)\n\n    f0 = fft2(f0)\n    f1 = fft2(f1)\n    r0 = abs(f0) * abs(f1)\n    ir = abs(ifft2((f0 * f1.conjugate()) / r0))\n    i0, i1 = numpy.unravel_index(numpy.argmax(ir), ir.shape)\n    angle = 180.0 * i0 / ir.shape[0]\n    scale = log_base ** i1\n\n    if scale > 1.8:\n        ir = abs(ifft2((f1 * f0.conjugate()) / r0))\n        i0, i1 = numpy.unravel_index(numpy.argmax(ir), ir.shape)\n        angle = -180.0 * i0 / ir.shape[0]\n        scale = 1.0 / (log_base ** i1)\n        if scale > 1.8:\n            raise ValueError(\"Images are not compatible. Scale change > 1.8\")\n\n    if angle < -90.0:\n        angle += 180.0\n    elif angle > 90.0:\n        angle -= 180.0\n\n    im2 = ndii.zoom(im1, 1.0/scale)\n    im2 = ndii.rotate(im2, angle)\n\n    if im2.shape < im0.shape:\n        t = numpy.zeros_like(im0)\n        t[:im2.shape[0], :im2.shape[1]] = im2\n        im2 = t\n    elif im2.shape > im0.shape:\n        im2 = im2[:im0.shape[0], :im0.shape[1]]\n\n    f0 = fft2(im0)\n    f1 = fft2(im2)\n    ir = abs(ifft2((f0 * f1.conjugate()) / (abs(f0) * abs(f1))))\n    t0, t1 = numpy.unravel_index(numpy.argmax(ir), ir.shape)\n\n    if t0 > f0.shape[0] // 2:\n        t0 -= f0.shape[0]\n    if t1 > f0.shape[1] // 2:\n        t1 -= f0.shape[1]\n\n    im2 = ndii.shift(im2, [t0, t1])\n\n    # correct parameters for ndimage's internal processing\n    if angle > 0.0:\n        d = int((int(im1.shape[1] / scale) * math.sin(math.radians(angle))))\n        t0, t1 = t1, d+t0\n    elif angle < 0.0:\n        d = int((int(im1.shape[0] / scale) * math.sin(math.radians(angle))))\n        t0, t1 = d+t1, d+t0\n    scale = (im1.shape[1] - 1) / (int(im1.shape[1] / scale) - 1)\n\n    return im2, scale, angle, [-t0, -t1]\n\n\ndef similarity_matrix(scale, angle, vector):\n    \"\"\"Return homogeneous transformation matrix from similarity parameters.\n\n    Transformation parameters are: isotropic scale factor, rotation angle (in\n    degrees), and translation vector (of size 2).\n\n    The order of transformations is: scale, rotate, translate.\n\n    \"\"\"\n    S = numpy.diag([scale, scale, 1.0])\n    R = numpy.identity(3)\n    angle = math.radians(angle)\n    R[0, 0] = math.cos(angle)\n    R[1, 1] = math.cos(angle)\n    R[0, 1] = -math.sin(angle)\n    R[1, 0] = math.sin(angle)\n    T = numpy.identity(3)\n    T[:2, 2] = vector\n    return numpy.dot(T, numpy.dot(R, S))\n\n\ndef logpolar(image, angles=None, radii=None):\n    \"\"\"Return log-polar transformed image and log base.\"\"\"\n    shape = image.shape\n    center = shape[0] / 2, shape[1] / 2\n    if angles is None:\n        angles = shape[0]\n    if radii is None:\n        radii = shape[1]\n    theta = numpy.empty((angles, radii), dtype=numpy.float64)\n    theta.T[:] = -numpy.linspace(0, numpy.pi, angles, endpoint=False)\n    #d = radii\n    d = numpy.hypot(shape[0]-center[0], shape[1]-center[1])\n    log_base = 10.0 ** (math.log10(d) / (radii))\n    radius = numpy.empty_like(theta)\n    radius[:] = numpy.power(log_base, numpy.arange(radii,\n                                                   dtype=numpy.float64)) - 1.0\n    x = radius * numpy.sin(theta) + center[0]\n    y = radius * numpy.cos(theta) + center[1]\n    output = numpy.empty_like(x)\n    ndii.map_coordinates(image, [x, y], output=output)\n    return output, log_base\n\n\ndef highpass(shape):\n    \"\"\"Return highpass filter to be multiplied with fourier transform.\"\"\"\n    x = numpy.outer(\n        numpy.cos(numpy.linspace(-math.pi/2., math.pi/2., shape[0])),\n        numpy.cos(numpy.linspace(-math.pi/2., math.pi/2., shape[1])))\n    return (1.0 - x) * (2.0 - x)\n\n\ndef imread(fname, norm=True):\n    \"\"\"Return image data from img&hdr uint8 files.\"\"\"\n    with open(fname+'.hdr', 'r') as fh:\n        hdr = fh.readlines()\n    img = numpy.fromfile(fname+'.img', numpy.uint8, -1)\n    img.shape = int(hdr[4].split()[-1]), int(hdr[3].split()[-1])\n    if norm:\n        img = img.astype(numpy.float64)\n        img /= 255.0\n    return img\n\n\ndef imshow(im0, im1, im2, im3=None, cmap=None, **kwargs):\n    \"\"\"Plot images using matplotlib.\"\"\"\n    from matplotlib import pyplot\n    if cmap is None:\n        cmap = 'coolwarm'\n    if im3 is None:\n        im3 = abs(im2 - im0)\n    pyplot.subplot(221)\n    pyplot.imshow(im0, cmap, **kwargs)\n    pyplot.subplot(222)\n    pyplot.imshow(im1, cmap, **kwargs)\n    pyplot.subplot(223)\n    pyplot.imshow(im3, cmap, **kwargs)\n    pyplot.subplot(224)\n    pyplot.imshow(im2, cmap, **kwargs)\n    pyplot.show()"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "We read in two files from the original set and compare them."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img1 = plt.imread(set175[1])\nimg2 = plt.imread(set175[3])\nplt.figure(figsize=(8,10))\nplt.subplot(121)\nplt.imshow(img1)\nplt.subplot(122)\nplt.imshow(img2)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "Now transform `img2` red band to align with `img1` red band."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img3, scale, angle, (t0, t1) = similarity(img1[:,:,2], img2[:,:,2])"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Viewing registration\n\nWe should be looking at:\n    \n* UL: template for transformation\n* UR: image that was transformed\n* LL: diff after transformation\n* LR: transformed image"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "imshow(img1[:,:,2], img2[:,:,2], img3)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# just diff\nplt.imshow(img3 - img1[:,:,2])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# well, was working in a standalone private notebook, need to figure out what's off.\n# it's fairly obviously offset/translation/shift is wrong somehow.\nplt.imshow(img3)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# NN with Downscaled Red Band Image\n\nI'm exploring the best way to feed images into a siamese net like the one described [here](http://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Zagoruyko_Learning_to_Compare_2015_CVPR_paper.pdf) and [here](http://www.cs.utoronto.ca/~gkoch/files/msc-thesis.pdf).\n\nIdeally you'd want to build the DNN to function this way:\n\n* Feed in two images as in diagram in paper [1]:\n\"A central-surround two-stream network that uses a\nsiamese-type architecture to process each stream\"\n* Train the DNN to learn a simple comparison function: 1 if image is before, 0 if false (no simultaneously captured images so no equality).\n* Once you have this comparator, then you can apply a sorting alrogithm to sort each dataset.\n\nThere are lots of details to be worked out here, still -- mainly trying to select overlapping image patches."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "from skimage.transform import downscale_local_mean\n\nimg1 = plt.imread(set175[1])[:,:,2]\nimg2 = plt.imread(set175[3])[:,:,2]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "np.shape(img1), np.shape(img2)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img1ds = downscale_local_mean(img1, (10, 10))\nimg2ds = downscale_local_mean(img2, (10, 10))\nplt.figure(figsize=(8,10))\nplt.subplot(121)\nplt.imshow(img1ds[:225,:300])\nplt.subplot(122)\nplt.imshow(img2ds[:225,:300])"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "For prototyping my goal is to get a bunch of downscaled images so that I can feed them in to test\nplausibility of basic NN architecture decisions.\n\nDims are subsampled and hard-coded for now so we don't have wrong dimensions showing up here and there."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "def read_and_downscale(f):\n    img = plt.imread(f)\n    ds = downscale_local_mean(img[:,:,2], (10, 10))\n    return ds[:225, :300]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "set79 = [read_and_downscale(\"../input/train_sm/set79_\" + str(n) + \".jpeg\") for n in range(1, 6)]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "def read_ds_set(setn):\n    match_pre = \"../input/train_sm/set\" + str(setn) + \"_\"\n    return [read_and_downscale(match_pre + str(n) + \".jpeg\") for n in range(1, 6)]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "set79 = read_ds_set(79)\nset285 = read_ds_set(285)\nset35 = read_ds_set(35)\nset175 = read_ds_set(175)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(set79[0])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(set285[2])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.imshow(set35[3])"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# pairing function\n\nWe want to traverse multiple lists of images and construct all naive before/after pairs.\n\n* Note there are actually more before/after pairs in this space, but this blows up quickly, at\n  least for little Kaggle notebooks."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "def get_pairs(imgls):\n    pairl = []\n    for imgl in imgls:\n        pairl += [(a,b) for a,b in zip(imgl[:-1], imgl[1:])]\n    return pairl\n            \npaired = get_pairs([set35, set79, set285, set175])"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# Let's just start with a dorky example.\nrev_pairs = [(imgb, imga) for imga, imgb in paired]"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "img_a, img_b = rev_pairs[0]\nconcat_img = np.vstack((img_a, img_b))\nprint(np.shape(concat_img))\nplt.imshow(concat_img)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "import random\n\ndef concatter(imgpairs):\n    for a, b in imgpairs:\n        yield np.vstack((a, b))\n\nconcats = [cimg for cimg in concatter(paired + rev_pairs)]\nrandom.shuffle(concats)"
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "# Oops!\n\nShuffled before supplying labels elsewhere about whether it's forward or backward comparison.\n\nGoing to have to leave this for tonight anyways."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.figure(figsize=(10,15))\nplt.subplot(321)\nplt.imshow(concats[0])\nplt.subplot(322)\nplt.imshow(concats[1])\nplt.subplot(323)\nplt.imshow(concats[2])\nplt.subplot(324)\nplt.imshow(concats[3])\nplt.subplot(325)\nplt.imshow(concats[4])\nplt.subplot(326)\nplt.imshow(concats[5])\nplt.show()"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "# Need to simplify or build up next block from convnet template."
 },
 {
  "cell_type": "markdown",
  "metadata": {},
  "source": "```python\nfrom keras.models import Sequential\nfrom keras.layers.core import Dense, Dropout, Activation, Flatten\nfrom keras.layers.convolutional import Convolution2D, MaxPooling2D\n\n# Need to add dropouts\nmodel = Sequential()\n\n# conv layer 1\nmodel.add(Convolution2D(50,1,2,2))\nmodel.add(Activation('relu'))\n\n# conv layer 2\nmodel.add(Convolution2D(50, 32, 2, 2))\nmodel.add(Activation('relu')) \nmodel.add(MaxPooling2D(poolsize=(2,2)))\n\n# conv layer 3\nmodel.add(Convolution2D(50, 50, 2, 2))\nmodel.add(Activation('relu'))\nmodel.add(MaxPooling2D(poolsize=(2,2)))\n\n# feed to fully connected \nmodel.add(Flatten())\n\n# first fully connected\nmodel.add(Dense(1000, 128, init='glorot_uniform'))\nmodel.add(Activation('relu'))\nmodel.add(Dropout(0.25))\n\n# next fully connected\nmodel.add(Dense(128, 64, init='glorot_uniform'))\nmodel.add(Activation('relu'))\nmodel.add(Dropout(0.25))\n\n# last fully connected which outputs comparison result\nmodel.add(Dense(64, 1, init='glorot_uniform'))\nmodel.add(Activation('sigmoid'))\n\n# compile model\n# model.compile(loss='binary_crossentropy', optimizer=\"rmsprop\")\n```"
 }
],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}}, "nbformat": 4, "nbformat_minor": 0}