svank/Python Earth Mover's Distance.ipynb

## Python Earth Mover's Distance.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trying to understand the Python wrapper for OpenCV 3's EMD function, whose C++ documentation is [here](https://docs.opencv.org/3.4.3/d6/dc7/group__imgproc__hist.html#ga902b8e60cc7075c8947345489221e0e0). (Fun fact, OpenCV's Python bindings are [automatically generated](https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_bindings/py_bindings_basics/py_bindings_basics.html), so Python documentation isn't guaranteed.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EMD function expects \"signatures\" rather than images/matrices. For an N-dimensional matrix with a total of M elements, the signature is an M x (N+1) array. Each of the M rows corresponds to a single pixel/element in the original image/matrix. Within each row, the first value is the pixel/element value that occurs in the source image/matrix, and the remaining N values are the coordinates along each dimension of that pixel/element. In short, it’s a list of pixel values and their corresponding coordinates.\n",
"\n",
"We can generate this signature by iterating through the values in our source image and filling in the signature’s rows one-by-one. The order we go through the source image pixels doesn’t matter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def img_to_sig(arr):\n",
    "    \"\"\"Convert a 2D array to a signature for cv2.EMD\"\"\"\n",
    "    \n",
    "    # cv2.EMD requires single-precision, floating-point input\n",
    "    sig = np.empty((arr.size, 3), dtype=np.float32)\n",
    "    count = 0\n",
    "    for i in range(arr.shape[0]):\n",
    "        for j in range(arr.shape[1]):\n",
    "            sig[count] = np.array([arr[i,j], i, j])\n",
    "            count += 1\n",
    "    return sig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2. 0. 0.]\n",
      " [0. 0. 1.]\n",
      " [0. 0. 2.]\n",
      " [2. 1. 0.]\n",
      " [0. 1. 1.]\n",
      " [0. 1. 2.]\n",
      " [0. 2. 0.]\n",
      " [0. 2. 1.]\n",
      " [2. 2. 2.]]\n",
      "[[0. 0. 0.]\n",
      " [1. 0. 1.]\n",
      " [1. 0. 2.]\n",
      " [2. 1. 0.]\n",
      " [0. 1. 1.]\n",
      " [0. 1. 2.]\n",
      " [0. 2. 0.]\n",
      " [2. 2. 1.]\n",
      " [0. 2. 2.]]\n"
     ]
    }
   ],
   "source": [
    "arr1 = np.array([[2, 0, 0],\n",
    "                 [2, 0, 0],\n",
    "                 [0, 0, 2]])\n",
    "\n",
    "arr2 = np.array([[0, 1, 1],\n",
    "                 [2, 0, 0],\n",
    "                 [0, 2, 0]])\n",
    "\n",
    "sig1 = img_to_sig(arr1)\n",
    "sig2 = img_to_sig(arr2)\n",
    "\n",
    "print(sig1)\n",
    "print(sig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of `cv2.EMD` has three values. The first is the distance between the two signatures. This appears to be normalized in some way---adding non-moving elements will reduce the distance, and doubling all pixel values doesn't affect the distance. I'm not sure what the second element---a `None`---is. The third value is the \"flow matrix\", telling you what was moved where. Per the documentation, if the input arrays (before being converted to signatures) have `size1` and `size2` elements each, the flow is a \"`𝚜𝚒𝚣𝚎𝟷×𝚜𝚒𝚣𝚎𝟸` flow matrix: `𝚏𝚕𝚘𝚠_i,j` is a flow from i-th point of signature1 to j-th point of signature2 .\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8333333134651184\n",
      "None\n",
      "[[0. 1. 1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 2. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 2. 0.]]\n"
     ]
    }
   ],
   "source": [
    "dist, _, flow = cv2.EMD(sig1, sig2, cv2.DIST_L2)\n",
    "\n",
    "print(dist)\n",
    "print(_)\n",
    "print(flow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So that 1 in row 1, column 2 is saying that one unit of \"earth\" was moved from the coordinates in row 1 of the first signature, or (0, 0), to the coordinates in row 2 of the second signature, or (0, 1). Note that unmoved earth shows up in this flow matrix---along the diagonal in this case, since the coordinates are in the same order in the two signatures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to visualize this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_flow(sig1, sig2, flow, arrow_width_scale=3):\n",
    "    \"\"\"Plots the flow computed by cv2.EMD\n",
    "    \n",
    "    The source images are retrieved from the signatures and\n",
    "    plotted in a combined image, with the first image in the\n",
    "    red channel and the second in the green. Arrows are\n",
    "    overplotted to show moved earth, with arrow thickness\n",
    "    indicating the amount of moved earth.\"\"\"\n",
    "    \n",
    "    img1 = sig_to_img(sig1)\n",
    "    img2 = sig_to_img(sig2)\n",
    "    combined = np.dstack((img1, img2, 0*img2))\n",
    "    # RGB values should be between 0 and 1\n",
    "    combined /= combined.max()\n",
    "    print('Red channel is \"before\"; green channel is \"after\"; yellow means \"unchanged\"')\n",
    "    plt.imshow(combined)\n",
    "\n",
    "    flows = np.transpose(np.nonzero(flow))\n",
    "    for src, dest in flows:\n",
    "        # Skip the pixel value in the first element, grab the\n",
    "        # coordinates. It'll be useful later to transpose x/y.\n",
    "        start = sig1[src, 1:][::-1]\n",
    "        end = sig2[dest, 1:][::-1]\n",
    "        if np.all(start == end):\n",
    "            # Unmoved earth shows up as a \"flow\" from a pixel\n",
    "            # to that same exact pixel---don't plot mini arrows\n",
    "            # for those pixels\n",
    "            continue\n",
    "        \n",
    "        # Add a random shift to arrow positions to reduce overlap.\n",
    "        shift = np.random.random(1) * .3 - .15\n",
    "        start = start + shift\n",
    "        end = end + shift\n",
    "        \n",
    "        mag = flow[src, dest] * arrow_width_scale\n",
    "        plt.quiver(*start, *(end - start), angles='xy',\n",
    "                   scale_units='xy', scale=1, color='white',\n",
    "                   edgecolor='black', linewidth=mag/3,\n",
    "                   width=mag, units='dots',\n",
    "                   headlength=5,\n",
    "                   headwidth=3,\n",
    "                   headaxislength=4.5)\n",
    "    \n",
    "    plt.title(\"Earth moved from img1 to img2\")\n",
    "    \n",
    "def sig_to_img(sig):\n",
    "    \"\"\"Convert a signature back to a 2D image\"\"\"\n",
    "    intsig = sig.astype(int)\n",
    "    img = np.empty((intsig[:, 1].max()+1, intsig[:, 2].max()+1), dtype=float)\n",
    "    for i in range(sig.shape[0]):\n",
    "        img[intsig[i, 1], intsig[i, 2]] = sig[i, 0]\n",
    "    return img"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAADHCAYAAADxqlPLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFcZJREFUeJzt3Xu4HVV5x/HvzySA3NFEDLkQeIwCUm8cwVttWsQCgvGxVIPKrWLUlgoWqxQFAa1SSrFSUAyCXEQuosVgsYgPUEABOaHcQgqmKOY0EUJCQsI98PaPWYfs7Oxz3bMvZ9bv8zz7OXv2rL3Wmj3vvGdmzew9igjMzCwvL+t0B8zMrP2c/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENO/iWS9FFJP+9g++dIOqGkuqZLWitpXJq+UdKRZdSd6vuZpMPKqs+aJ+kkSd9PzzdY/x3oS6nbkqSFkmal5y8tZ0l1Hy/pu2XV1y6VTf6S3iXpV5JWS1op6ZeS3prmHS7plibrnyEpJI3vfy0iLomI9zbb9wHa+52kpyWtkbQqLdunJL20DiPiUxHxlWHW9Z7BykTE7yNiy4h4oYS+b7SxRcR+EXFhs3XnriYu1tY8zmq23jLXfz1JF0h6LsXyGkn3Sfq6pG1q2h/WtpTq+upQ5SLi9RFxY5NdR9IsSX11dX8tIkrbMWqXSiZ/SVsDPwX+DXgFMAU4GXi2pPrHD12qJQ6MiK2AHYFTgS8A55XdSAeXz0bnwJSo+x9HdbpDw3BaiuVJwBHA24BfStqizEYcy4OIiMo9gB5g1QDzdgWeAV4A1vaXA94H/DfwBLAEOKnmPTOAAD4O/B64Kf2NVMda4O3A4cAtNe8L4FPAb4DHgbMBpXnjgH8BHgN+CxyVyo8foN+/A95T99qewIvA7mn6AuCr6flEin+Aq4CVwM0U/+wvTu95OvX78wMsX/9r41N9NwJfB34NrAZ+ArwizZsF9DXqL7Av8BzwfGrv7pr6jkzPXwZ8CXgYeBS4CNim7rM/LPXtMeCLnY6xbnk0iouaeYcDtwCnp/j7LbBfzfydgP8C1gDXAWcB36/73GvX/1eAX6byPwcm1tR1aFp/K4AThujXS3Fa89pWwDLgqNq+p+cCvpFiYzVwD7A7MDfF1XMptq6u+Uy+kMo9C4yv7Q9wEnAlcHlaljuBN9Ztt6+p7y+wBcV28yLrt/sdUn3fryn/fmAhxbZ3I7Br3fr6XOrb6tSHzToRO5Xc8wceBF6QdKGk/SRt1z8jIhZRJORbo9hL2jbNepIigLel+EfwaUkfqKv3Tyj+efw58O702rapnlsH6MsBwFuBNwIfSu8F+ASwH/Am4C1AfVtDiohfA33AHzeYfWyaNwnYHji+eEscQpFE+/cWTxtg+Ro5FPgrioBfB5w5jD7+J/A14PLU3hsbFDs8Pf4U2BnYkiIR1XoX8Dpgb+BESbsO1bYBsBfwAMXOwGnAeZKU5v0AWJDmfYXiH+xgPkKxl/4qYBOKJIak3YBvAR8FJgPbUBxtD1tE9P8DahTL76XY3l5LsX1+GFgREfOASyiOIraMiANr3nMwxXa8bUSsa1DnbOCHFCMDPwCukjRhiD4+SbHNLo31R1lLa8tIei1wKXAMxbZ3DXC1pE1qin2IYqdoJ+ANFLHfdpVM/hHxBEWyCOBcYLmk+ZK2H+Q9N0bEvRHxYkTcQ7EC/6Su2EkR8WREPD2C7pwaEasi4vfADRTJHooA+GZE9EXE4xTDOKOxlCKA6z1PsSHuGBHPR8TNkXY9BjHU8l0cEfeljeAE4EMlnRD8KHBGRDwUEWuBfwDm1B2ynxwRT0fE3cDdFP9MrXBVOg/U//hEzbyHI+LcKMbuL6SIie0lTafYKTkhIp6NiJuAq4do53sR8WCKjytYH8sHUex13xIRzwEnUmx7IzVYLG8F7EJx5LwoIpYNUdeZEbFkkFheEBFXRsTzwBnAZhRDT836MPAfEXFdqvt04OXAO+r6tjQiVlJ85m9qUE/LVTL5Q7GHHxGHR8RUikPEHYB/Hai8pL0k3SBpuaTVFEcHE+uKLRlFV/5Q8/wpir1aUn9q6xtN3VDsYa1s8Po/A4uBn0t6SNJxw6hrqD7Uzn8YmMDGn9Fo7JDqq617PMURS7+BPkeDD0TEtjWPc2vmvfS5RcRT6emWFJ/54+kfeb/addDIsGI5tbNihMsAA8RyRFxPcSR4NvCIpHnpvN5ghh3LEfEixVHyDiPrbkMbxHKqewkbHgl1RSxXNvnXioj/oRi3273/pQbFfgDMB6ZFxDbAORRjjRtUNcDz0VgGTK2ZnjbSCtLVS1MoxnU3EBFrIuLYiNgZOBD4O0l7988eoMqhlqm2j9Mp9sgeoxgy27ymX+MoDnmHW+9SipPYtXWvAx4Z4n02esuA7epOsE5voq6XYlnSy4FXjqQCSVtSnCO6udH8iDgzIvYAXk8x/PP3/bMGqHLYsZyumJtKEYdQJOTNa8q+egT1bhDLaYhtGvB/Q7yv7SqZ/CXtIulYSVPT9DSKMcDbUpFHgKl143BbASsj4hlJe1KMbw5mOcWJn51H2c0rgKMlTZG0LcUJqmGRtLWkA4DLKE403dugzAGSXpOC7wmKE9z9l+09Msp+f0zSbpI2B04BrkzDCQ8Cm0l6Xxo3/RKwac37HgFm1F6WWudS4LOSdkpJoP8cQaOxWitBRDwM9AInS9pE0rsodhJG40rgQEnvSNvUyWy849SQpE0l7QFcRXFS+nsNyrw1HZlPoNjR6L9gA0Yfy3tI+mAaWjyG4sRwf364C/iIpHGS9mXD4d9HgFfWXpZa5wrgfZL2Tv09NtX9q1H0saUqmfwpzuDvBdwu6UmKlXofxYoAuJ7ibPwfJD2WXvtr4BRJayjGLK8YrIF0aPuPFJenrZI00vHCcymumLiH4iqjayj2dge7rvrq1L8lwBcpxiqPGKDsTOAXFFck3Ap8K9Zf5/x14Eup358bQZ8vpjiC+gPFGOlnACJiNcXn912KPZwnKQ6j+/0w/V0h6c4G9Z6f6r6J4oqUZ4C/HUG/cnd13XX+/z7M932EYjtZCXyZ4iqrEYuIhRTr6zKKo4A1FFfmDHZp9edTLK9M7S4A3lE3DNVva4rt5XHWX1F0epp3HrBbiuWrRtDtn1CMzz8OHAJ8MI3RAxxN8Y9wFcX5qJfqTaMIlwIPpTY3GCqKiAeAj1FcZv5YqufAdC6kq/RfdmgdJmk/4JyI2HHIwmZdLB29rQJmRsRvO90fa6yqe/5dT9LLJe0vabykKRR7XsPdYzPrKpIOlLR5OodwOnAvxTXt1qWaSv6SXiHpOkm/SX+3G6DcC5LuSo/5zbRZIaIYG32cYthnEcVwk3UBx/aIzaY42bmUYshxzjAuLbYOamrYR9JpFCdJT02XEm4XERuduJS0NiJ8aZ6NGY5tq7pmk/8DwKyIWCZpMnBjRLyuQTlvIDamOLat6pod89++/5t26e+rBii3maReSbdp459MMOtGjm2rtCF/8U7SL9jwSw79vjiCdqZHxFJJOwPXS7o3Iv63QVtzKX6siS222GKPXXbZZQRNdK8FCxZ0ugvW2PMU31Go19LYBvYYeVe70+TJkzvdBauzatUqnnrqqSG/ZzFk8o+IAX/3XdIjkibXHBo/OkAdS9PfhyTdCLwZ2GgDST/UNA+gp6cnent7h+remCAN6/su1n73RERPoxmtjG1JlTkR+slPfrLTXbA63/nOd4ZVrtlhn/ms/yXAwyi+OLEBSdtJ2jQ9nwi8E7i/yXbNWs2xbZXWbPI/FdhH0m+AfdI0knq0/rZmuwK9ku6m+FXLUyPCG4h1O8e2VVpTd7mJiBUUv69e/3ovcGR6/ivgj5ppx6zdHNtWdf6Gr5lZhpz8zcwy5ORvZpYhJ38zsww5+ZuZZcjJ38wsQ07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENO/mZmGXLyNzPLkJO/mVmGSkn+kvaV9ICkxZKOazB/U0mXp/m3S5pRRrtmrebYtqpqOvlLGgecDewH7AYcLGm3umIfBx6PiNcA3wD+qdl2zVrNsW1VVsae/57A4oh4KCKeAy4DZteVmQ1cmJ5fCewtSSW0bdZKjm2rrDKS/xRgSc10X3qtYZmIWAesBl5ZX5GkuZJ6JfUuX768hK6ZNaUlsd2ivpqNSBnJv9FeToyiDBExLyJ6IqJn0qRJJXTNrCktie1SembWpDKSfx8wrWZ6KrB0oDKSxgPbACtLaNuslRzbVlllJP87gJmSdpK0CTAHmF9XZj5wWHp+EHB9RGy0d2TWZRzbVlnjm60gItZJOgq4FhgHnB8RCyWdAvRGxHzgPOBiSYsp9ormNNuuWas5tq3Kmk7+ABFxDXBN3Wsn1jx/BvjLMtoyayfHtlWVv+FrZpYhJ38zsww5+ZuZZcjJ38wsQ07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENO/mZmGXLyNzPLkJO/mVmGnPzNzDJUSvKXtK+kByQtlnRcg/mHS1ou6a70OLKMds1azbFtVdX0nbwkjQPOBvahuJn1HZLmR8T9dUUvj4ijmm3PrF0c21ZlZez57wksjoiHIuI54DJgdgn1mnWaY9sqq4x7+E4BltRM9wF7NSj3F5LeDTwIfDYiltQXkDQXmAswffr0Erpm1pSWxfbDDz/cgu62n6ROd8FGqYw9/0ZrP+qmrwZmRMQbgF8AFzaqKCLmRURPRPRMmjSphK6ZNcWxbZVVRvLvA6bVTE8FltYWiIgVEfFsmjwX2KOEds1azbFtlVVG8r8DmClpJ0mbAHOA+bUFJE2umXw/sKiEds1azbFtldX0mH9ErJN0FHAtMA44PyIWSjoF6I2I+cBnJL0fWAesBA5vtl2zVnNsW5Upon4Iszv09PREb29vp7tRCp8U61oLIqKn3Y06tq3VImLIFeNv+JqZZcjJ38wsQ07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENO/mZmGXLyNzPLkJO/mVmGnPzNzDLk5G9mliEnfzOzDJWS/CWdL+lRSfcNMF+SzpS0WNI9kt5SRrtmreS4tiora8//AmDfQebvB8xMj7nAt0tq16yVLsBxbRVVSvKPiJso7l86kNnARVG4Ddi27sbXZl3HcW1V1q4x/ynAkprpvvSa2VjmuLYxq13Jv9HNhDe6c7ykuZJ6JfUuX768Dd0ya8qw4hoc29Z92pX8+4BpNdNTgaX1hSJiXkT0RETPpEmT2tQ1s1EbVlyDY9u6T7uS/3zg0HR1xNuA1RGxrE1tm7WK49rGrPFlVCLpUmAWMFFSH/BlYAJARJwDXAPsDywGngKOKKNds1ZyXFuVlZL8I+LgIeYH8DdltGXWLo5rqzJ/w9fMLENO/mZmGXLyNzPLkJO/mVmGnPzNzDLk5G9mliEnfzOzDDn5m5llyMnfzCxDTv5mZhly8jczy5CTv5lZhpz8zcwy5ORvZpYhJ38zsww5+ZuZZcjJ38wsQ6Ukf0nnS3pU0n0DzJ8labWku9LjxDLaNWslx7VVWSm3cQQuAM4CLhqkzM0RcUBJ7Zm1wwU4rq2iStnzj4ibgJVl1GXWLRzXVmVl7fkPx9sl3Q0sBT4XEQvrC0iaC8ytmW5j92w4inuWV0NJ8TVkXKe2KhnbjoexS2WtPEkzgJ9GxO4N5m0NvBgRayXtD3wzImYOUV91oqpCKraxL4iIniHKzKDEuE7vq8yHWLF46HQXShMRQy5MW672iYgnImJten4NMEHSxHa0bdYqjmsby9qS/CW9WunfqqQ9U7sr2tG2Was4rm0sK2XMX9KlwCxgoqQ+4MvABICIOAc4CPi0pHXA08CcqNLxolWS49qqrLQx/7JVaVy0Sro1XkZjOGP+LWq3Mh9ixeKh010oTdeM+ZuZWXdx8jczy5CTv5lZhpz8zcwy5ORvZpYhJ38zsww5+ZuZZcjJ38wsQ07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENNJ39J0yTdIGmRpIWSjm5QRpLOlLRY0j2S3tJsu2at5ti2KivjNo7rgGMj4k5JWwELJF0XEffXlNkPmJkeewHfTn/Nuplj2yqr6T3/iFgWEXem52uARcCUumKzgYuicBuwraTJzbZt1kqObauyUsf8Jc0A3gzcXjdrCrCkZrqPjTciJM2V1Cupt8x+mTXLsW1VU8awDwCStgR+BBwTEU/Uz27wlo3u/BwR84B5qb7q3BnaxjTHtlVRKXv+kiZQbByXRMSPGxTpA6bVTE8FlpbRtlkrObatqsq42kfAecCiiDhjgGLzgUPTlRFvA1ZHxLJm2zZrJce2VVkZwz7vBA4B7pV0V3rteGA6QEScA1wD7A8sBp4CjiihXbNWc2xbZSmiO4cfPS7anbo1XkZD0oKI6OlAu5X5ECsWD53uQmkiYsiF8Td8zcwy5ORvZpYhJ38zsww5+ZuZZcjJ38wsQ07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXIyd/MLENO/mZmGXLyNzPLkJO/mVmGyriN4zRJN0haJGmhpKMblJklabWku9LjxGbbNWs1x7ZVWRm3cVwHHBsRd0raClgg6bqIuL+u3M0RcUAJ7Zm1i2PbKqvpPf+IWBYRd6bna4BFwJRm6zXrNMe2VVmpY/6SZgBvBm5vMPvtku6W9DNJry+zXbNWc2xb1ZQx7AOApC2BHwHHRMQTdbPvBHaMiLWS9geuAmY2qGMuMDdNrgUeKKt/g5gIPNaGdtqh5cvSpptct2ud7DicQmM0ttvyGVYsHtqhHcsyvLiOiKZbkjQB+ClwbUScMYzyvwN6IqLjK1RSb0T0dLofZajKsnTTcozV2O6mz7BZXpbWKONqHwHnAYsG2jgkvTqVQ9Keqd0VzbZt1kqObauyMoZ93gkcAtwr6a702vHAdICIOAc4CPi0pHXA08CcKOOQw6y1HNtWWU0n/4i4BRh04C8izgLOaratFpnX6Q6UqCrL0hXLMcZjuys+w5J4WVqglDF/MzMbW/zzDmZmGco2+UvaV9IDkhZLOq7T/RktSedLelTSfZ3uS7OG83MKNjTHdnfp1rjOcthH0jjgQWAfoA+4Azi4wdf2u56kd1NcN35RROze6f40Q9JkYHLtzykAHxiL66VTHNvdp1vjOtc9/z2BxRHxUEQ8B1wGzO5wn0YlIm4CVna6H2XwzymUwrHdZbo1rnNN/lOAJTXTfXTByrD1hvg5BRuYY7uLdVNc55r8G12+l9/4V5ca4ucUbHCO7S7VbXGda/LvA6bVTE8FlnaoL1Yj/ZzCj4BLIuLHne7PGOTY7kLdGNe5Jv87gJmSdpK0CTAHmN/hPmVvOD+nYENybHeZbo3rLJN/RKwDjgKupTj5ckVELOxsr0ZH0qXArcDrJPVJ+nin+9SE/p9T+LOaO2Pt3+lOjSWO7a7UlXGd5aWeZma5y3LP38wsd07+ZmYZcvI3M8uQk7+ZWYac/M3MMuTkb2aWISd/M7MMOfmbmWXo/wG8A7DePaSMyAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe9b5b81c18>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red channel is \"before\"; green channel is \"after\"; yellow means \"unchanged\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEICAYAAABBKnGGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAGeFJREFUeJzt3XmUXGWd//H3h6wEkbBpFkIACYgKcwwRcVzIEJiBQIBx/M0PRAUBm6AgztFRjnhMwijgnBkHGR0zMMOPZTAgeyL7FmBkkU6EhEVCAkjahDUQEhKC3fn+/rg3SdGp7up+6tbSnc8r556uqvvUfb51q+pT9z51b0oRgZlZb23V6ALMrG9yeJhZEoeHmSVxeJhZEoeHmSVxeJhZEodHL0jaTVJIGtjoWqoh6QVJh3Qxb2tJcyStlHRNA2r7vqT/qne/RZN0vKQ7Gl1HLfWr8MjfFGslrS6Zfl7l8sq+yfqxLwAfBHaMiP9T784j4tyIOKUWy5b095IelLRG0twKbSdKakvtKyKujIi/Tr1/dyT9o6QnJK2S9Lykf6xFP5X06U/QLkyJiLuqWYCkgRHRXlRBfcxYYFFXj7+Pr5sVwAXAh4GDG1xLNQR8BVgAfAi4Q9LSiLiqrlVERL+ZgBeAQ7qY9yHgHuB14DXgSmB4p/t+L39C1gGzgPXAWmA18F1gNyCAE4AX8+Wc3U09lwL/AdyaL+O3wAiyF/AbwB+Aj5e03weYC7wJPAkcld9+IPASMKCk7d8CC/LLWwFnAUvyx/drYIeStl8G/pjPO7ur9QTMAN4F/pzXezJwYl73v5G9+X6U9/eDfJmvAJcD2+XL2LCOvgoszR/nVOAT+bp9E/h5N+tsOvA/KcsCBgD/mj8vzwOn5/cf2KmPU4C53dSwTf68r8/Xw2pgFDAkf+6W5dMFwJAulnEi8L8l1wP4OvAssAr4J7LX5EPAW/lzNrik/XeB5Xk/p+T337OLvi4E/r3u77dGv+ELfTDdh8eewKH5C2Bn4H7ggk73fQwYA2xdbnklL+aLga2BvyALmn266PPS/IW8PzCULLyeJ/vUGJC/Ee/N2w4CFgPfBwaTfTKuAvbO5y8BDi1Z9jXAWfnlbwEPA7vkj+8/gVn5vI/kL/7P5fN+CrR3s56mk795S94E7cAZZFuqWwMn5bXuAbwPuB64otM6mpk/5r8G3gFuBD4AjCYLnIMq9d/bZZEFy1P5etgeuIuE8MjbTATaOt12Tr6eP0D2GnoQ+Kcu7n8im4fHbOD9wEfz183d+TrcLq/7hLztYWQfFh8FhgFX0EV4kG2F/B6YWvf3W6Pf8IU+mOzNvprsE2nD9LUu2h4D/L7TfU8qs7xy4bFLyW2/A47too9LgYtLrp8BPF1yfV/gzfzyZ/MXzFYl82cB0/PLPwIuyS9vC7wNjM2vPw1MKrnfSLKth4HAD4GrSuZtQ7Z10ZvweLFTm7uBr5dc37ukvw3raHTJ/NeB/1ty/TrgW5X67+2yyML51JJ5h1BseCwBJpdc/xvghS7ufyKbh8enS67PA75Xcv1fyT/MgEuA80rm7UnX4TEDeJwutoBqOfXHMY9josyYh6QPkG3efZbszbcV2WZwqaU97OOlkstryD59u/JyyeW1Za5vuO8oYGlErC+Z/0eyT1eAXwEPSjoN+DwwPyL+mM8bC9wgqfS+HWQDn6MoeVwR8bak17upt5zO62VUXltpnQPz/jbo6ePuiV6tw5J5PX0+e6rc4x7Vi/tXehwjSvppLZlX9nFIOp1sK/azEbGuF3UUol9921LBeWTpvV9EvB/4EtkmX6nOpxjX85TjZcAYSaXPya7AnwAi4imyF+vhwBfJwmSDpcDhETG8ZBoaEX8i228es6GhpGHAjr2srfN6WEYWWKV1tvPeN0MjLCfbZdlgTFcNe6Dcc1/ucS+roo+uVHwckk4iG+eaFBHJ3wpVY0sKj23Jd2kkjQZ68vXWy2T7pPXwCNmuyHclDZI0EZgClI6g/wr4Jtn4RekxGDOBH0saCyBpZ0lH5/OuBY6U9BlJg8n226t93mcB/yBpd0nvA84Fro7Gfwvza+BMSaMlDScbAN9I0gBJQ8m2kraSNFTSoC6W9TKwo6TtSm6bBfwgX787ke0S/k/xD4NfA1+VtE8e9j8snSnpeLJ1fmhEPFeD/nukP4bHnE7HedyQ3z4DGA+sBG4mG+Sr5DyyF8ubkr5To3oBiIh3gaPItixeI/uW5isR8YeSZrPI9sXviYjXSm7/Gdlg3B2SVpEN6n0yX+6TwDfIgmc52a5atZ9Ul5AN4t1PNgD8Dtl4TqNdDNxB9k3M74FbyLaIOvL5XybbPfgl2e7r2vw+m8nX+yzgufz5H0U27tSaL38hMD+/rVARcSvZLva9ZAPTD+WzNuya/Ihs6/HRktf5zKLrqET5oItZvyPpcGBmRIyt2LiJSdoHeIJsULTRW3cb9cctD9tC5YfWT5Y0MN81nQbcUOl+zUjS30oaLGl74CfAnGYKDqgyPCTtIOlOSc/mf7fvol2HpMfyaXY1fZp1Q2S7p2+Q7bY8Tafxgj7kVOBVsq+HO4DTGlvO5qrabZH0z8CKiDhf0lnA9hHxvTLtVkdEb76aM7MmV214PANMjIjlkkaSHXizd5l2Dg+zfqba8HgzIoaXXH8jIjbbdZHUTnbodztwfkTc2MXyWoAWgG1g/w8nV9b/zRvZ6AqsX1jOaxGxc8pdKx5hKukuNh35VursXvSza0Qsk7QHcI+khRGxpHOjiLgIuAhgghStnRvYRmppdAXWL8x4zxGzvVIxPCKiy//PQtLLkkaW7La80sUyluV/n8v/H4WPkw0EmVkfVe1XtbPJTk8n/3tT5waStpc0JL+8E/BpsjMIzawPqzY8zgcOlfQs2enu5wNImlDyX8ntA7RKepzsiLnz8/M0zKwPq+qs2oh4HZhU5vZWstOeiYgHyU49N7N+xEeYmlkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4mFkSh4eZJXF4bCleBR4n+yltswJU9T+JWR/yGnAjDN1mKINHDmbV2FXEXgE7k/3OmlkvOTy2FB+GbUZvwxX/fgXDhg3j2huu5YYbb+Cd9nf4855/5t0PvQtj8SvCeqyqH32qpXr+bsvcfOpLZhyUcKel8LHhH2PhvIUARAQLFy7kptk3cdX1V7Fk0RIG7TmI1buthnHANoWWbM1oBvMiYkLKXR0eZFvt06ZNq1NvjTVixAimTp1adt7LL7/MLbfcwqzrZvHA3Ae8e7MlcHhUR2SfwrbJW2+9xfTp07nwwgvp6OiA44C9Gl2VFa6K8PAebm769OmNLqEuKm153Hzzzfzqul/x2/t+y+CRg1n/V+uzXZikXzO1/sxbHsCFwIo69VWUosY8FixYsGnM49klDBo3iLd3ezsLjGGFlmzNyLstWx71dogmYJv/zr5t2Xrrrbnmhmu4ac5NvNPR6duWAbWo1pqWd1usoj/A2396m+O+dBxDRg1h1W6riGM8EGrpHB5bip2AY2DduHWsG7au0dVYP1DI4emSDpP0jKTFks4qM3+IpKvz+Y9I2q2Ifq0Xdgb+Ao9jWGGqDg9JA4BfAIcDHwGOk/SRTs1OBt6IiD2BfwN+Um2/ZtZYRWx5HAAsjojnIuJd4Crg6E5tjgYuyy9fC0yS5D1t6/vmA21Ac37vUFNFjHmMBpaWXG8DPtlVm4hol7QS2JHsdC2zPmtY6zDWr1wPg0B7ibW7r4UPAYMbXVntFbHlUW4LonMO96QNkloktUpqfbWAwsxqbc2n1rDXR/fiiXlPcN7x5/GJZZ9g8AWD2faabeFRYGWjK6ydIrY82oAxJdd3AZZ10aZN0kBgO8oclxURFwEXQXacRwG11cxcGnwyXUM7t43Ww8L5C3nxxRc588wzOfPMM1m5ciW33347V113FXf8vzvYarutWLP7Gjr27Mi2wfvJDnvVB4nlYbAImAT8iSxvvxgRT5a0+Qawb0RMlXQs8PmI+PvultvsB4ltSSfTWWXHH38848aN2+z2jo4OHnroIa6/8XquufEaXnvttebavWn0EaaSJgMXkB2feElE/FjSOUBrRMyWNBS4Avg42RbHsRHxXHfL7Avh0axH51rzuvPOOzn99NNZtGgRA8YMoOOkjsYW1OjwqIW+EB7e8rANutryaG9vf8+Wx4oVK4hxwTt7vAN70Ke3PHyEaaKfAStmzGhY/0knxlnx1oMeFAcddNDG8Fi5ciW33XYbV19/NXfcdgdbDc/HPA7ugFH0mzEPh0eibza4/xkTG1yAZRbCvuP3Zdddd+WCCy7gymuvZMH8BQzZYwirxq6Ck8i+HuiHHB5mVRj20DAWrVzEvhP23bQ78jl4d/C7jS6t5hweZlVYM2ENfJB+tTvSUw4Ps2qMb3QBjeMffTKzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJA4PM0vi8DCzJIWEh6TDJD0jabGks8rMP1HSq5Iey6dTiujXzBqn6h99kjQA+AVwKNAGPCppdkQ81anp1RFxerX9mVlzKGLL4wBgcUQ8FxHvAlcBRxewXDNrYkX83ORoYGnJ9Tbgk2Xa/Z2kzwGLgH+IiKWdG0hqAVoAdt0V+GMB1fVXW9jvolrzKWLLo9zLODpdnwPsFhH7AXcBl5VbUERcFBETImLCzjsXUJmZ1UwR4dEGjCm5vguwrLRBRLweEevyqxcD+xfQr5k1UBHh8SgwTtLukgYDxwKzSxtIGlly9Sjg6QL6NbMGqnrMIyLaJZ0O3A4MAC6JiCclnQO0RsRs4JuSjgLagRXAidX2a2aNpYjOwxPNYcIERWtro6toXvKAqRVjXkRMSLmjjzA1syQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsyQODzNL4vAwsySFhIekSyS9IumJLuZL0oWSFktaIGl8Ef2aWeMUteVxKXBYN/MPB8blUwvwy4L6NbMGKSQ8IuJ+YEU3TY4GLo/Mw8BwSSOL6NvMGqNeYx6jgaUl19vy295DUoukVkmtr75ap8rMLEm9wkNlbovNboi4KCImRMSEnXeuQ1Vmlqxe4dEGjCm5vguwrE59m1kN1Cs8ZgNfyb91ORBYGRHL69S3mdXAwCIWImkWMBHYSVIbMA0YBBARM4FbgMnAYmAN8NUi+jWzxikkPCLiuArzA/hGEX2ZWXPwEaZmlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlsThYWZJHB5mlqSQ8JB0iaRXJD3RxfyJklZKeiyfflhEv2bWOIX80DVwKfBz4PJu2jwQEUcW1J+ZNVghWx4RcT+woohlmVnfUNSWR098StLjwDLgOxHxZOcGklqAlk3X61hdXxONLqD5hV8/FVWziuoVHvOBsRGxWtJk4EZgXOdGEXERcBGAJL89zJpYXb5tiYi3ImJ1fvkWYJCknerRt5nVRl3CQ9IIKdsJkXRA3u/r9ejbzGqjkN0WSbOAicBOktqAacAggIiYCXwBOE1SO7AWODYivFti1oepWd/DHvOowGunIg+YViaYFxETUu7rI0zNLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODys71lO9huFdze4ji1cPX/0ySxNAI8Bc4DfAI/mtw8H3mhUUebwsOa0FriXTYHRVqaNf/m4oRwe1jyWAzeTBcZdwJpNs0aNGsWRRx7JvvvuyxlnnJHtcP+gIVVazuFhjdPV7khu//33Z8qUKRx55JGMHz8eSbS05D9l/EVg77pWa534d1v6qr66drrZHRk6dCiHHHIIU6ZM4YgjjmD06NHvuesLL7zAuHHjaF/fDk9RMTz8uy2VVfO7Ld7ysNrrwe7IlClTOPjggxk2bFiXizn33HNpb2+HMcCsyt1Or6ro4mwLnEw2vtufeMujr2rmtZOwO9ITI0eO5KWXXiq42Pr4F+DbjS6iDG95WONVsTvSU1deeSX3339/AcXWz9y5c7nvvvtY1ehCaiEiqprINiLvBZ4GngTOLNNGwIXAYmABML4Hyw1P3UzN8u8egqMIhr23vlGjRkVLS0vMmTMn3n777dhSTZs2LYCYBhFNOAGtkfjeL2LLox34dkTMl7QtME/SnRHxVEmbw4Fx+fRJ4Jf5X+vrTgCWbrq6ww47cN5553HyySczYMCAhpVltVf14ekRsTwi5ueXV5FtgXTeLj0auDwP44eB4ZJGVtu3NYHLyJ7dfJxzxYoVnHrqqYwZM4aWlhbmzJnDmjVruluC9VGFDphK2g24H/hYRLxVcvtvgPMj4n/z63cD34uI1m6W1cxDgo3XbGvnHTaNecxhszGPSZMmbRwkTR3zuPvuu3nggQcKKLZ+ZsyYAcA0mufbn1LVDJhWPeaxYQLeB8wDPl9m3s3AZ0qu3w3sX6ZdC9CaT40fV2jmqZn/rSf4PcE5BAdsXvv48eNj2rRp8eijj0ZHR0ePxw9GjBjR+PWeON1L48c3yk1UMeZRyJaHpEFkY+y3R8RPy8z/T2BuRMzKrz8DTIyI5d0ss/rC+rO+tHZeYtNxHnfynuM8Ro4cufE4j0mTJnV7nEdLSwsXX3xxNkR/UuVup82oruyiTMynZtTQLY+sfy4HLuimzRHArXnbA4Hf9WC5Df+0aOqpr/5bS3ALwWkEu7z3MQ0dOjSOOOKImDlzZrS1tW225fH888/HwIEDg60I/lC5r0Z/qveFiSq2PIoIj8/kT/4CskODHgMmA1OBqSUB8wtgCbAQmODwqHLqD/8Sdm++9rWvZfO/VHn5jX5j9oWpmvDwEaZ9VX9cOz3Yvdlvv/02nVVb4fwWn9tSWTW7LQ6Pvqq/r51uvr3Z6EvAFV0vwuFRmcNjS7QlrZ0AHmfToe+/y2+v8D+JOTwqc3hsibbktfMScBvZty6Tum7m8KjM4bEl8tqpyOFRWTXh4f893cySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLInDw8ySODzMLEnV4SFpjKR7JT0t6UlJZ5ZpM1HSSkmP5dMPq+3XzBprYAHLaAe+HRHzJW0LzJN0Z0Q81andAxFxZAH9mVkTqHrLIyKWR8T8/PIq4GlgdLXLNbPmVsSWx0aSdgM+DjxSZvanJD0OLAO+ExFPlrl/C9CSX10HPFFkfQXYCXit0UUAkP0Oa/PUk2mqetRk9eSaraa9U+9Y2A9dS3ofcB/w44i4vtO89wPrI2K1pMnAzyJiXIXltab+AG+tNFtNrqd7zVYPNF9N1dRTyLctkgYB1wFXdg4OgIh4KyJW55dvAQZJ2qmIvs2sMYr4tkXAfwNPR8RPu2gzIm+HpAPyfl+vtm8za5wixjw+DXwZWCjpsfy27wO7AkTETOALwGmS2oG1wLFReX/pogJqK1qz1eR6utds9UDz1ZRcT2FjHma2ZfERpmaWxOFhZkmaJjwk7SDpTknP5n+376JdR8lh7rNrUMdhkp6RtFjSWWXmD5F0dT7/kfzYlprqQU0nSnq1ZL2cUsNaLpH0iqSyx+Aoc2Fe6wJJ42tVSy9qqtvpET08XaOu66hmp5BERFNMwD8DZ+WXzwJ+0kW71TWsYQCwBNgDGAw8DnykU5uvAzPzy8cCV9d4vfSkphOBn9fpefocMB54oov5k4FbyQ5jOxB4pAlqmgj8pk7rZyQwPr+8LbCozPNV13XUw5p6vY6aZssDOBq4LL98GXBMA2o4AFgcEc9FxLvAVXldpUrrvBaYtOFr6AbWVDcRcT+wopsmRwOXR+ZhYLikkQ2uqW6iZ6dr1HUd9bCmXmum8PhgRCyH7MECH+ii3VBJrZIellR0wIwGlpZcb2PzlbyxTUS0AyuBHQuuo7c1Afxdvgl8raQxNaynkp7WW2+fkvS4pFslfbQeHXZzukbD1lFPTiHp6Toq9NyWSiTdBYwoM+vsXixm14hYJmkP4B5JCyNiSTEVUm4LovN32T1pU6Se9DcHmBUR6yRNJdsyOriGNXWn3uunJ+YDY2PT6RE3At2eHlGt/HSN64BvRcRbnWeXuUvN11GFmnq9juq65RERh0TEx8pMNwEvb9h0y/++0sUyluV/nwPmkqVoUdqA0k/tXchO5CvbRtJAYDtqu8lcsaaIeD0i1uVXLwb2r2E9lfRkHdZV1Pn0iEqna9CAdVSLU0iaabdlNnBCfvkE4KbODSRtL2lIfnknsqNbO/+/IdV4FBgnaXdJg8kGRDt/o1Na5xeAeyIfcaqRijV12l8+imyftlFmA1/Jv1E4EFi5YXe0Uep5ekTeT7ena1DnddSTmpLWUT1GoHs4IrwjcDfwbP53h/z2CcB/5Zf/ElhI9o3DQuDkGtQxmWw0eglwdn7bOcBR+eWhwDXAYuB3wB51WDeVajoPeDJfL/cCH65hLbOA5cCfyT5BTwamAlPz+QJ+kde6EJhQh/VTqabTS9bPw8Bf1rCWz5DtgiwAHsunyY1cRz2sqdfryIenm1mSZtptMbM+xOFhZkkcHmaWxOFhZkkcHmaWxOFhZkkcHmaW5P8D1L1P6R7PnTIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe9b5a505c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.imshow(arr1, cmap='gray')\n",
    "ax1.set_title(\"Starting Distribution\")\n",
    "ax2.imshow(arr2, cmap='gray')\n",
    "ax2.set_title(\"Ending Distribution\")\n",
    "plt.show()\n",
    "\n",
    "plot_flow(sig1, sig2, flow)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we see for these two distributions, 2 units of earth are distributed from the upper-left corner to the other top-row cells, another two units move along the bottom row, and a final two units in the middle row are unmoved."
   ]
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}