{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "a373440b-e2b5-4e97-a85f-cad130bd1d8f"
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import sklearn\n",
"import glob, os"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "565430d3-8ae0-41e7-8d70-29865ac0bd2e"
},
"source": [
"# Read in files\n",
"\n",
"This 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": {
"_cell_guid": "8946eb4b-0f02-4873-9262-b4a0bb73a0fa"
},
"outputs": [],
"source": [
"from subprocess import check_output\n",
"print(check_output([\"ls\", \"../input\"]).decode(\"utf8\"))\n",
"smjpegs = [f for f in glob.glob(\"../input/train_sm/*.jpeg\")]\n",
"print(smjpegs[:9])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "18c5a6a8-e8fb-4484-993c-d0010ae663f5"
},
"outputs": [],
"source": [
"set175 = [smj for smj in smjpegs if \"set175\" in smj]\n",
"print(set175)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "89334ce0-fe0b-43a7-be5e-34bf6cb12fdb"
},
"source": [
"# Basic exploration\n",
"\n",
"Just look at image dimensions, confirm it's 3 band (RGB), byte scaled (0-255)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "503c0083-c8da-4e7f-8f35-b133b9ff49b0"
},
"outputs": [],
"source": [
"first = plt.imread('../input/train_sm/set175_1.jpeg')\n",
"dims = np.shape(first)\n",
"print(dims)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "3c71599d-94c7-485b-8af5-bb83feba036f"
},
"outputs": [],
"source": [
"np.min(first), np.max(first)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "8b04fce9-658f-48eb-b76d-84d6a946ba24"
},
"source": [
"For any image specific classification, clustering, etc. transforms we'll want to \n",
"collapse spatial dimensions so that we have a matrix of pixels by color channels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "c23dfd15-42f3-469f-af9b-d922fc090653"
},
"outputs": [],
"source": [
"pixel_matrix = np.reshape(first, (dims[0] * dims[1], dims[2]))\n",
"print(np.shape(pixel_matrix))"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "02fbb025-2467-4b63-8b0a-b6ff63a485e5"
},
"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\n",
"switch to using 2d histograms instead. The data between bands is really correlated, typical with\n",
"visible imagery and why most satellite image analysts prefer to at least have near infrared values."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "76af2f59-4ef3-4659-b8fd-bedd68ecb7d3"
},
"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": {
"_cell_guid": "2d7b7e4a-7646-4f4f-aaea-9be4bbb8025c"
},
"outputs": [],
"source": [
"fifth = plt.imread('../input/train_sm/set175_5.jpeg')\n",
"dims = np.shape(fifth)\n",
"pixel_matrix5 = np.reshape(fifth, (dims[0] * dims[1], dims[2]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e4e91466-81fa-4470-85d4-a0774e128f48"
},
"outputs": [],
"source": [
"_ = plt.hist2d(pixel_matrix5[:,1], pixel_matrix5[:,2], bins=(50,50))"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "360f6da1-0fb6-41c1-9240-3d31ee5a4d63"
},
"source": [
"We can look at variations between the scenes now and see that there's a significant\n",
"amount of difference, probably due to sensor angle and illumination variation. Raw band\n",
"differences will need to be scaled or thresholded for any traditional approach."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e7c1e14d-fb66-4145-ba0a-06f168ce754f"
},
"outputs": [],
"source": [
"_ = plt.hist2d(pixel_matrix[:,2], pixel_matrix5[:,2], bins=(50,50))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "39f4e335-a1d5-4462-9ba3-a416854066a8"
},
"outputs": [],
"source": [
"plt.imshow(first)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "23e0d226-28ca-449b-bcec-38198b0e0ffa"
},
"outputs": [],
"source": [
"plt.imshow(fifth)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "c4194a14-3b71-4e68-af18-2a81f497761e"
},
"source": [
"Without coregistering portions of the image, the naive red band subtraction for change indication\n",
"basically just shows the location shift between images."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "21cfae7f-b49f-43b3-a62f-48364379c176"
},
"outputs": [],
"source": [
"plt.imshow(first[:,:,2] - fifth[:,:,1])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "a9be7930-35e1-4d15-a292-10324a7c7966"
},
"outputs": [],
"source": [
"second = plt.imread('../input/train_sm/set175_2.jpeg')\n",
"plt.imshow(first[:,:,2] - second[:,:,2])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "71a18cc6-f596-4ed8-9fc8-2bf2cc6c3d76"
},
"outputs": [],
"source": [
"plt.imshow(second)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "cf88bba9-8e4f-4c03-b584-a6f4e007454a"
},
"source": [
"# Initial impressions\n",
"\n",
"Images aren't registered, so an image registration process between images with common overlap would probably be the first step in a traditional approach.\n",
"Using a localizer in a deep learning context would probably be the newfangled way to tackle this.\n",
"\n",
"Image content and differences will be dominated by topographic and built variations\n",
"due to sensor orientation, resolution differences between scenes, and some registration accuracy will be impossible to factor out as\n",
"the image hasn't been orthorectified and some anciliary data would be required for it\n",
"to be done, e.g. georeferenceing against a previously orthorectified image.\n",
"\n",
"So this is basically a basic computer vision task that deep learning will be a good fit for. The usual preprocessing steps\n",
"and data expectations you'd see in remote sensing aren't fulfilled by this dataset."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "b68f3943-d96b-47d6-af52-295de21d22da"
},
"outputs": [],
"source": [
"# simple k means clustering\n",
"from sklearn import cluster\n",
"\n",
"kmeans = cluster.KMeans(5)\n",
"clustered = kmeans.fit_predict(pixel_matrix)\n",
"\n",
"dims = np.shape(first)\n",
"clustered_img = np.reshape(clustered, (dims[0], dims[1]))\n",
"plt.imshow(clustered_img)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "d885beba-5b77-4d2a-8953-de327dfb204a"
},
"outputs": [],
"source": [
"plt.imshow(first)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "43c0250c-b7f3-4528-9a93-37a983040f59"
},
"outputs": [],
"source": [
"ind0, ind1, ind2, ind3 = [np.where(clustered == x)[0] for x in [0, 1, 2, 3]]"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "551a4dfc-20a0-40fd-90de-2882d1659071"
},
"source": [
"This code doesn't run on the server.\n",
"\n",
"```python\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"\n",
"plot_vals = [('r', 'o', ind0),\n",
" ('b', '^', ind1),\n",
" ('g', '8', ind2),\n",
" ('m', '*', ind3)]\n",
"\n",
"for 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",
"\n",
"ax.set_xlabel('Blue channel')\n",
"ax.set_ylabel('green channel')\n",
"ax.set_zlabel('Red channel')\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "0edcda33-0862-46ba-97f9-3babaeb95ed6"
},
"outputs": [],
"source": [
"# quick look at color value histograms for pixel matrix from first image\n",
"import seaborn as sns\n",
"sns.distplot(pixel_matrix[:,0], bins=12)\n",
"sns.distplot(pixel_matrix[:,1], bins=12)\n",
"sns.distplot(pixel_matrix[:,2], bins=12)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "6cdc7d36-9e4a-45dd-85b0-9e3b0797fe20"
},
"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": {
"_cell_guid": "5726c6ab-5f09-4dab-9592-3a1183396f04"
},
"source": [
"# Day 2\n",
"\n",
"We'll start by considering the entire sequence of a different image set this time and look at strategies\n",
"for matching features across scenes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "b425bdfd-91f6-491e-b5ed-50628effb3d5"
},
"outputs": [],
"source": [
"set79 = [smj for smj in smjpegs if \"set79\" in smj]\n",
"print(set79)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "06a5a947-54a5-42e9-8f9b-abaf1e5b8894"
},
"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": {
"_cell_guid": "d4e80ef3-b20d-45ab-8c98-39efb3b4595e"
},
"outputs": [],
"source": [
"img_list = (img79_1, img79_2, img79_3, img79_4, img79_5)\n",
"\n",
"plt.figure(figsize=(8,10))\n",
"plt.imshow(img_list[0])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "39192a54-b42e-4eda-b9dd-062a8574476b"
},
"source": [
"Tracking dimensions across image transforms is annoying, so we'll make a class to do that.\n",
"Also I'm going to use this brightness normalization transform and visualize the image that\n",
"way, good test scenario for class."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "bd278077-037d-4ccb-a862-8979d59d08ed"
},
"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": {
"_cell_guid": "873cc990-36eb-4bf1-8fdf-7d093fbafee7"
},
"outputs": [],
"source": [
"msi79_1 = MSImage(img79_1)\n",
"print(np.shape(msi79_1.matrix))\n",
"print(np.shape(msi79_1.img))"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "98e8078c-6054-458c-9f1b-543978a8f460"
},
"source": [
"# Brightness Normalization\n",
"\n",
"Brightness Normalization is preprocessing strategy you can apply prior to using strategies\n",
"to identify materials in a scene, if you want your matching algorithm\n",
"to 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": {
"_cell_guid": "d1e21cc4-f57c-4440-a991-328d1cea055d"
},
"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": {
"_cell_guid": "a0461571-4d3c-48d8-bf90-efe6c84ad63b"
},
"outputs": [],
"source": [
"bnorm = bnormalize(msi79_1.matrix)\n",
"bnorm_img = msi79_1.to_matched_img(bnorm)\n",
"plt.figure(figsize=(8,10))\n",
"plt.imshow(bnorm_img)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "9bea4215-28f4-4adf-8b72-de0e45944913"
},
"outputs": [],
"source": [
"msi79_2 = MSImage(img79_2)\n",
"bnorm79_2 = bnormalize(msi79_2.matrix)\n",
"bnorm79_2_img = msi79_2.to_matched_img(bnorm79_2)\n",
"plt.figure(figsize=(8,10))\n",
"plt.imshow(bnorm79_2_img)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "b9bb8d22-3e4c-404e-8398-4dd55efffd37"
},
"outputs": [],
"source": [
"msinorm79_1 = MSImage(bnorm_img)\n",
"msinorm79_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": {
"_cell_guid": "cf4e849b-3c38-4f70-8a25-f1a4feb7ff43"
},
"outputs": [],
"source": [
"_ = plt.hist2d(msinorm79_1.matrix[:,1], msinorm79_2.matrix[:,1], bins=(50,50))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "10ea0358-e917-4e3e-96dc-6d7b53334111"
},
"outputs": [],
"source": [
"_ = plt.hist2d(msinorm79_1.matrix[:,0], msinorm79_2.matrix[:,0], bins=(50,50))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e4ba466d-7b75-4154-92e5-c3404e2dab28"
},
"outputs": [],
"source": [
"import seaborn as sns\n",
"sns.distplot(msinorm79_1.matrix[:,0], bins=12)\n",
"sns.distplot(msinorm79_1.matrix[:,1], bins=12)\n",
"sns.distplot(msinorm79_1.matrix[:,2], bins=12)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "a8357cc4-a69d-4765-89ef-8bb0b4852f26"
},
"outputs": [],
"source": [
"plt.figure(figsize=(8,10))\n",
"plt.imshow(img79_1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "632d9902-dc30-4b46-b5fd-3a315573d738"
},
"outputs": [],
"source": [
"np.max(img79_1[:,:,0])"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "ce163ec2-7d77-4b64-b3ef-46f9c091a29d"
},
"source": [
"# Using thresholds with brightness normalization\n",
"\n",
"Ok, so what am I even doing here? Well, my goal is to try and figure out simple threshold selection\n",
"methods for getting high albedo targets out of a scene so I could then theoretically track them\n",
"between scenes. For example, a simple blob/aggregation to centroid (in coordinates or in subsampled\n",
"image bins) would give me a means to look at plausible structural similarities in distributions\n",
"between scenes, then use that to anchor a comparison of things that change.\n",
"\n",
"The brightness normalization step is helpful because thresholds that aren't anchored by a\n",
"preprocessing step end up being arbitrary and can't generalize between scenes even in the same\n",
"image set, whereas thresholds following brightness normalization tend to pull out materils that stand\n",
"out from the background more reliably. See the following demonstration:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "cbac949a-6815-40af-b675-2967598aab03"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow(img79_1[:,:,0] > 230)\n",
"plt.subplot(122)\n",
"plt.imshow(img79_1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "711f6e12-42db-4972-af87-065e8d30b304"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow(img79_2[:,:,0] > 230)\n",
"plt.subplot(122)\n",
"plt.imshow(img79_2)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "6b2d5512-bc1f-4870-a4e8-38041b868ed6"
},
"outputs": [],
"source": [
"print(np.min(bnorm79_2_img[:,:,0]))\n",
"print(np.max(bnorm79_2_img[:,:,0]))\n",
"print(np.mean(bnorm79_2_img[:,:,0]))\n",
"print(np.std(bnorm79_2_img[:,:,0]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "5d14565f-daf2-47ae-9830-ab2570ce57aa"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow(bnorm79_2_img[:,:,0] > 0.98)\n",
"plt.subplot(122)\n",
"plt.imshow(img79_2)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "9c92a2fa-bb58-4407-82fd-6fc10b7cb6ae"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow(bnorm_img[:,:,0] > 0.98)\n",
"plt.subplot(122)\n",
"plt.imshow(img79_1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "65bd07d2-e6c4-4409-a50e-6c906936b76e"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow((bnorm79_2_img[:,:,0] > 0.9999) & \\\n",
" (bnorm79_2_img[:,:,1] < 0.9999) & \\\n",
" (bnorm79_2_img[:,:,2] < 0.9999))\n",
"plt.subplot(122)\n",
"plt.imshow(img79_2)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "36abf15c-5b0f-48ce-a77b-5852c53b973c"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(121)\n",
"plt.imshow(bnorm_img[:,:,0] > 0.995)\n",
"plt.subplot(122)\n",
"plt.imshow(img79_1)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "4ba15127-8c37-43a2-af6c-a65f0005c931"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,6))\n",
"plt.subplot(121)\n",
"plt.plot(bnorm_img[2000, 1000, :])\n",
"plt.subplot(122)\n",
"plt.plot(img79_1[2000, 1000, :])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "0b864974-539f-4efd-a11e-8063e8803fb3"
},
"outputs": [],
"source": [
"from scipy import spatial\n",
"\n",
"pixel = msi79_1.matrix[2000 * 1000, :]\n",
"np.shape(pixel)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "0c95ed52-9018-4585-8479-93e6b34c4208"
},
"source": [
"# Something's borked here\n",
"\n",
"Think I'm gonna have to verify cosine similarity behavior for scipy here.\n",
"\n",
"```python\n",
"def spectral_angle_mapper(pixel):\n",
" return lambda p2: spatial.distance.cosine(pixel, p2)\n",
"\n",
"match_pixel = np.apply_along_axis(spectral_angle_mapper(pixel), 1, msi79_1.matrix)\n",
"\n",
"plt.figure(figsize=(10,6))\n",
"plt.imshow(msi79_1.to_flat_img(match_pixel < 0.0000001))\n",
"\n",
"def summary(mat):\n",
" print(\"Max: \", np.max(mat),\n",
" \"Min: \", np.min(mat),\n",
" \"Std: \", np.std(mat),\n",
" \"Mean: \", np.mean(mat))\n",
" \n",
"summary(match_pixel)\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "18ccdfb5-0937-4abf-b344-b370dce391d4"
},
"source": [
"# Rudimentary Transforms, Edge Detection, Texture"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "05470414-7685-43f8-9e91-2ce478f4fce3"
},
"outputs": [],
"source": [
"set144 = [MSImage(plt.imread(smj)) for smj in smjpegs if \"set144\" in smj]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "43f32ef1-6748-43ea-a3ed-8ee92cebd996"
},
"outputs": [],
"source": [
"plt.imshow(set144[0].image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "4016c910-da89-4973-8601-5198ac5b242d"
},
"outputs": [],
"source": [
"import skimage\n",
"from skimage.feature import greycomatrix, greycoprops\n",
"from skimage.filters import sobel"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "0e1ebeee-6567-4246-ba19-fe0945f223c4"
},
"source": [
"# Sobel Edge Detection\n",
"\n",
"A Sobel filter is one means of getting a basic edge magnitude/gradient image. Can be useful to\n",
"threshold and find prominent linear features, etc. Several other similar filters in skimage.filters\n",
"are also good edge detectors: `roberts`, `scharr`, etc. and you can control direction, i.e. use\n",
"an anisotropic version."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "67e36e21-ba5b-43e9-9cbb-a99f0fc14d75"
},
"outputs": [],
"source": [
"# a sobel filter is a basic way to get an edge magnitude/gradient image\n",
"fig = plt.figure(figsize=(8, 8))\n",
"plt.imshow(sobel(set144[0].image[:750,:750,2]))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "99a2ce33-1228-4c91-b634-bac2eef99c56"
},
"outputs": [],
"source": [
"from skimage.filters import sobel_h\n",
"\n",
"# can also apply sobel only across one direction.\n",
"fig = plt.figure(figsize=(8, 8))\n",
"plt.imshow(sobel_h(set144[0].image[:750,:750,2]), cmap='BuGn')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "49af3363-56e4-425c-af4c-f9d8e6c489cd"
},
"outputs": [],
"source": [
"from sklearn.decomposition import PCA\n",
"\n",
"pca = PCA(3)\n",
"pca.fit(set144[0].matrix)\n",
"set144_0_pca = pca.transform(set144[0].matrix)\n",
"set144_0_pca_img = set144[0].to_matched_img(set144_0_pca)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "a4e61335-6594-4a8f-8411-529f883ac1ac"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.imshow(set144_0_pca_img[:,:,0], cmap='BuGn')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e3d99ae2-7512-4d32-9c0a-a26d1562f586"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.imshow(set144_0_pca_img[:,:,1], cmap='BuGn')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "2a41fe2d-0d25-4395-a66e-64c9fa3e29fc"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.imshow(set144_0_pca_img[:,:,2], cmap='BuGn')"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "1e2db1e2-474a-4df1-9c00-92938b2b1c2e"
},
"source": [
"# GLCM Textures\n",
"\n",
"Processing time can be pretty brutal so we subset the image. We'll create texture images so\n",
"we can characterize each pixel by the texture of its neighborhood.\n",
"\n",
"GLCM 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",
"\n",
"A 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",
"\n",
"Pesaresi, 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": {
"_cell_guid": "e7cfc63e-1d26-4cfd-84c6-f93cc1ff07f6"
},
"outputs": [],
"source": [
"sub = set144[0].image[:150,:150,2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "b9525d5a-6e7c-42af-8b1e-960a6af6e43b"
},
"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": {
"_cell_guid": "c1de324f-f12f-4c63-a8cd-b5c21b220ae3"
},
"outputs": [],
"source": [
"dissimilarity = glcm_image(sub, \"dissimilarity\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "26b2e488-2226-4772-a598-85cd27541221"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.subplot(1,2,1)\n",
"plt.imshow(dissimilarity, cmap=\"bone\")\n",
"plt.subplot(1,2,2)\n",
"plt.imshow(sub, cmap=\"bone\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "8b8fac31-86a0-4694-809a-facf0c9ba11c"
},
"source": [
"# HSV Transform\n",
"\n",
"Since this contest is about time series ordering, I think it's possible there may be useful\n",
"information in a transform to HSV color space. HSV is useful for identifying shadows and illumination, as well\n",
"as giving us a means to identify similar objects that are distinct by color between scenes (hue), \n",
"though there's no guarantee the hue will be stable."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "45ab5ebf-fe02-4cdc-889e-6a7a2978fe83"
},
"outputs": [],
"source": [
"from skimage import color\n",
"\n",
"hsv = color.rgb2hsv(set144[0].image)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "4c74a90f-ca3a-4819-bae4-68764f17d69a"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.subplot(2,2,1)\n",
"plt.imshow(set144[0].image, cmap=\"bone\")\n",
"plt.subplot(2,2,2)\n",
"plt.imshow(hsv[:,:,0], cmap=\"bone\")\n",
"plt.subplot(2,2,3)\n",
"plt.imshow(hsv[:,:,1], cmap='bone')\n",
"plt.subplot(2,2,4)\n",
"plt.imshow(hsv[:,:,2], cmap='bone')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "1e6f1ee1-3dea-47ef-9092-7e6f0d26b8d0"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.subplot(2,2,1)\n",
"plt.imshow(set144[0].image[:200,:200,:])\n",
"plt.subplot(2,2,2)\n",
"plt.imshow(hsv[:200,:200,0], cmap=\"PuBuGn\")\n",
"plt.subplot(2,2,3)\n",
"plt.imshow(hsv[:200,:200,1], cmap='bone')\n",
"plt.subplot(2,2,4)\n",
"plt.imshow(hsv[:200,:200,2], cmap='bone')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "009dcfd7-d380-4c45-90ba-cb3d2da7cb95"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 6))\n",
"plt.imshow(hsv[200:500,200:500,0], cmap='bone')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "053bb5b8-cad1-4aeb-b4d4-e6bce2727abd"
},
"outputs": [],
"source": [
"hsvmsi = MSImage(hsv)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "f2ad8c07-961a-4715-9e2f-c803f0a0285e"
},
"source": [
"# Shadow Detection\n",
"\n",
"We can apply a threshold to the V band now to find dark areas that are probably thresholds. Let's\n",
"look at the distribution of all values then work interactively to find a good filter value."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "38dbfc30-0451-4f52-8b3b-2b5c2e840d94"
},
"outputs": [],
"source": [
"import seaborn as sns\n",
"sns.distplot(hsvmsi.matrix[:,0], bins=12)\n",
"sns.distplot(hsvmsi.matrix[:,1], bins=12)\n",
"sns.distplot(hsvmsi.matrix[:,2], bins=12)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "d564ff67-5706-4875-bbf9-22c9f379ad80"
},
"outputs": [],
"source": [
"plt.imshow(hsvmsi.image[:,:,2] < 0.4, cmap=\"plasma\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "90b69fb2-2131-4c80-b966-ba67f59bdfb4"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"plt.subplot(1,2,1)\n",
"plt.imshow(set144[0].image[:250,:250,:])\n",
"plt.subplot(1,2,2)\n",
"plt.imshow(hsvmsi.image[:250,:250,2] < 0.4, cmap=\"plasma\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "ec864f2d-1487-4abe-8b93-96ad81e21901"
},
"outputs": [],
"source": [
"fig = plt.figure(figsize=(8, 8))\n",
"img2 = plt.imshow(set144[0].image[:250,:250,:], interpolation='nearest')\n",
"img3 = plt.imshow(hsvmsi.image[:250,:250,2] < 0.4, cmap='binary_r', alpha=0.4)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "87c6cb06-84d8-407f-9187-587f4736c77e"
},
"source": [
"Could we glean something useful about sun position from shadow orientation if we could accurately\n",
"reference the image?"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "974be7a5-2424-4664-83e6-93f367ae46f3"
},
"source": [
"# Image Registration\n",
"\n",
"This is an earlier form the library found [here](https://github.com/matejak/imreg_dft).\n",
"\n",
"BSD family license, reproduced below with copyright so I can utilize similar functions here where\n",
"import isn't available.\n",
"\n",
"This version can be found [here](http://www.lfd.uci.edu/~gohlke/code/imreg.py.html)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "de05542e-7376-40be-a42d-433f6da252ae"
},
"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",
"\n",
"Implements an FFT-based technique for translation, rotation and scale-invariant\n",
"image registration [1].\n",
"\n",
":Author:\n",
" `Christoph Gohlke `_\n",
"\n",
":Organization:\n",
" Laboratory for Fluorescence Dynamics, University of California, Irvine\n",
"\n",
":Version: 2013.01.18\n",
"\n",
"Requirements\n",
"------------\n",
"* `CPython 2.7 or 3.3 `_\n",
"* `Numpy 1.7 `_\n",
"* `Scipy 0.12 `_\n",
"* `Matplotlib 1.2 `_ (optional for plotting)\n",
"\n",
"Notes\n",
"-----\n",
"The API and algorithms are not stable yet and are expected to change between\n",
"revisions.\n",
"\n",
"References\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",
"\n",
"Examples\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",
"\n",
"from __future__ import division, print_function\n",
"\n",
"import math\n",
"\n",
"import numpy\n",
"from numpy.fft import fft2, ifft2, fftshift\n",
"\n",
"try:\n",
" import scipy.ndimage.interpolation as ndii\n",
"except ImportError:\n",
" import ndimage.interpolation as ndii\n",
"\n",
"__version__ = '2013.01.18'\n",
"__docformat__ = 'restructuredtext en'\n",
"__all__ = ['translation', 'similarity']\n",
"\n",
"\n",
"def 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",
"\n",
"def 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",
"\n",
"def 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",
"\n",
"def 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",
"\n",
"def 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",
"\n",
"def 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",
"\n",
"def 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": {
"_cell_guid": "e8f3c2c2-f748-4abe-a966-40496dd1814a"
},
"source": [
"We read in two files from the original set and compare them."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "f956fe2b-2267-4c37-a17e-be5396f6b50d"
},
"outputs": [],
"source": [
"img1 = plt.imread(set175[1])\n",
"img2 = plt.imread(set175[3])\n",
"plt.figure(figsize=(8,10))\n",
"plt.subplot(121)\n",
"plt.imshow(img1)\n",
"plt.subplot(122)\n",
"plt.imshow(img2)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "44866fde-8c95-40cc-8a72-ff4ab59d4d7b"
},
"source": [
"Now transform `img2` red band to align with `img1` red band."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "d8d6e9c1-66d7-4c2c-b992-8ee5acfd85d8"
},
"outputs": [],
"source": [
"img3, scale, angle, (t0, t1) = similarity(img1[:,:,2], img2[:,:,2])"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "4f328193-f623-4e69-80d0-075dfa1dfb7c"
},
"source": [
"# Viewing registration\n",
"\n",
"We 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": {
"_cell_guid": "4ed3a7a8-f7df-4162-aa49-a93e03869512"
},
"outputs": [],
"source": [
"imshow(img1[:,:,2], img2[:,:,2], img3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "12b4d780-4750-4222-a4dd-86d24dfadd24"
},
"outputs": [],
"source": [
"# just diff\n",
"plt.imshow(img3 - img1[:,:,2])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "24836f0c-92a7-40e1-a58a-845646d91b86"
},
"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.\n",
"plt.imshow(img3)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "2199af01-62bc-49ce-adff-0188f40c156d"
},
"source": [
"# NN with Downscaled Red Band Image\n",
"\n",
"I'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",
"\n",
"Ideally 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\n",
"siamese-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",
"\n",
"There are lots of details to be worked out here, still -- mainly trying to select overlapping image patches."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "8004d39e-ead9-4183-aeed-0479d79a18e0"
},
"outputs": [],
"source": [
"from skimage.transform import downscale_local_mean\n",
"\n",
"img1 = plt.imread(set175[1])[:,:,2]\n",
"img2 = plt.imread(set175[3])[:,:,2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "b95fa028-17c1-45b8-991e-f0fd7d3db7b0"
},
"outputs": [],
"source": [
"np.shape(img1), np.shape(img2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "4749673b-97a5-45e8-add0-14aca6b2b8a8"
},
"outputs": [],
"source": [
"img1ds = downscale_local_mean(img1, (10, 10))\n",
"img2ds = downscale_local_mean(img2, (10, 10))\n",
"plt.figure(figsize=(8,10))\n",
"plt.subplot(121)\n",
"plt.imshow(img1ds[:225,:300])\n",
"plt.subplot(122)\n",
"plt.imshow(img2ds[:225,:300])"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "596b03a5-d201-4b44-99e4-72e3686807a2"
},
"source": [
"For prototyping my goal is to get a bunch of downscaled images so that I can feed them in to test\n",
"plausibility of basic NN architecture decisions.\n",
"\n",
"Dims 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": {
"_cell_guid": "2617280a-fbd7-4bf6-b02e-d5f743862674"
},
"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": {
"_cell_guid": "77a30430-b5bd-439d-a1ef-478bbbad9a8d"
},
"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": {
"_cell_guid": "fbc9428a-8d3c-45ba-959d-e0fa5b043d1a"
},
"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": {
"_cell_guid": "8d48964c-0c2b-4af0-8de7-930ab0c7eacd"
},
"outputs": [],
"source": [
"set79 = read_ds_set(79)\n",
"set285 = read_ds_set(285)\n",
"set35 = read_ds_set(35)\n",
"set175 = read_ds_set(175)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "7836d280-50bd-4448-a849-168a344c4ebb"
},
"outputs": [],
"source": [
"plt.imshow(set79[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "9983cea1-16c7-4661-a737-22059f010d46"
},
"outputs": [],
"source": [
"plt.imshow(set285[2])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "50da6f39-1294-4ce4-921b-5463039444f6"
},
"outputs": [],
"source": [
"plt.imshow(set35[3])"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "5631c88c-57b5-40e4-9627-57da18036bd6"
},
"source": [
"# pairing function\n",
"\n",
"We 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": {
"_cell_guid": "553e8c7d-dd0c-4c7e-80b5-6a74ddc5d995"
},
"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",
" \n",
"paired = get_pairs([set35, set79, set285, set175])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "279d67fc-bd24-443c-899d-e788fdec3f94"
},
"outputs": [],
"source": [
"# Let's just start with a dorky example.\n",
"rev_pairs = [(imgb, imga) for imga, imgb in paired]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "7f2d9237-5298-483b-85fe-c24220479bf5"
},
"outputs": [],
"source": [
"img_a, img_b = rev_pairs[0]\n",
"concat_img = np.vstack((img_a, img_b))\n",
"print(np.shape(concat_img))\n",
"plt.imshow(concat_img)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "7fa9e0f0-b9c1-406e-aa88-0366635b34e8"
},
"outputs": [],
"source": [
"import random\n",
"\n",
"def concatter(imgpairs):\n",
" for a, b in imgpairs:\n",
" yield np.vstack((a, b))\n",
"\n",
"concats = [cimg for cimg in concatter(paired + rev_pairs)]\n",
"random.shuffle(concats)"
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "2b7c5d69-6b17-46c3-8bdc-5247927ce88c"
},
"source": [
"# Oops!\n",
"\n",
"Shuffled before supplying labels elsewhere about whether it's forward or backward comparison.\n",
"\n",
"Going to have to leave this for tonight anyways."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e319f204-0237-4655-a98b-8f466fdf4d65"
},
"outputs": [],
"source": [
"plt.figure(figsize=(10,15))\n",
"plt.subplot(321)\n",
"plt.imshow(concats[0])\n",
"plt.subplot(322)\n",
"plt.imshow(concats[1])\n",
"plt.subplot(323)\n",
"plt.imshow(concats[2])\n",
"plt.subplot(324)\n",
"plt.imshow(concats[3])\n",
"plt.subplot(325)\n",
"plt.imshow(concats[4])\n",
"plt.subplot(326)\n",
"plt.imshow(concats[5])\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"_cell_guid": "e16066be-01e8-4889-be3c-921f55cb2a57"
},
"outputs": [],
"source": [
"# Need to simplify or build up next block from convnet template."
]
},
{
"cell_type": "markdown",
"metadata": {
"_cell_guid": "e7348ce5-b862-42f4-8bb5-6cf30ac4a24c"
},
"source": [
"```python\n",
"from keras.models import Sequential\n",
"from keras.layers.core import Dense, Dropout, Activation, Flatten\n",
"from keras.layers.convolutional import Convolution2D, MaxPooling2D\n",
"\n",
"# Need to add dropouts\n",
"model = Sequential()\n",
"\n",
"# conv layer 1\n",
"model.add(Convolution2D(50,1,2,2))\n",
"model.add(Activation('relu'))\n",
"\n",
"# conv layer 2\n",
"model.add(Convolution2D(50, 32, 2, 2))\n",
"model.add(Activation('relu')) \n",
"model.add(MaxPooling2D(poolsize=(2,2)))\n",
"\n",
"# conv layer 3\n",
"model.add(Convolution2D(50, 50, 2, 2))\n",
"model.add(Activation('relu'))\n",
"model.add(MaxPooling2D(poolsize=(2,2)))\n",
"\n",
"# feed to fully connected \n",
"model.add(Flatten())\n",
"\n",
"# first fully connected\n",
"model.add(Dense(1000, 128, init='glorot_uniform'))\n",
"model.add(Activation('relu'))\n",
"model.add(Dropout(0.25))\n",
"\n",
"# next fully connected\n",
"model.add(Dense(128, 64, init='glorot_uniform'))\n",
"model.add(Activation('relu'))\n",
"model.add(Dropout(0.25))\n",
"\n",
"# last fully connected which outputs comparison result\n",
"model.add(Dense(64, 1, init='glorot_uniform'))\n",
"model.add(Activation('sigmoid'))\n",
"\n",
"# compile model\n",
"# model.compile(loss='binary_crossentropy', optimizer=\"rmsprop\")\n",
"```"
]
}
],
"metadata": {
"_change_revision": 0,
"_is_fork": false,
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}