{"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": "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": ""
 },
 {
  "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\nprint(\"Image \" + str(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": "I initially defined a @np.vectorize function for this, but the loop runs faster for some\nreason."
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "bnorm = np.zeros_like(msi79_1.matrix, dtype=np.float32)\nfor x in range(7219900):\n    bnorm[x,:] = msi79_1.matrix[x,:] / float(np.max(msi79_1.matrix[x,:]))"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "bnorm_img = msi79_1.to_matched_img(bnorm)\n\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)\n\ndef bnormalize(mat):\n    bnorm = np.zeros_like(mat, dtype=np.float32)\n    for x in range(np.shape(mat)[0]):\n        bnorm[x,:] = mat[x,:] / float(np.max(mat[x,:]))\n    return bnorm\n\nbnorm79_2 = bnormalize(msi79_2.matrix)\nbnorm79_2_img = msi79_2.to_matched_img(bnorm79_2)"
 },
 {
  "cell_type": "code",
  "execution_count": null,
  "metadata": {
   "collapsed": false
  },
  "outputs": [],
  "source": "plt.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": ""
 }
],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}}, "nbformat": 4, "nbformat_minor": 0}