Image Patches.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Patches\n",
    "\n",
    "In this module, you will explore the topology of different collections of image patches.  Each image patch is a square $d \\times d$ region of pixels.  Each pixel can be thought of as a dimension, so each patch lives in $\\mathbb{R}^{d \\times d}$, and a collection of patches can be thought of as a Euclidean point cloud in $\\mathbb{R}^{d \\times d}$\n",
    "\n",
    "First, we perform all of the necessary library imports.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.offsetbox import OffsetImage, AnnotationBbox\n",
    "from ripser import ripser, plot_dgms\n",
    "from GeomUtils import getGreedyPerm\n",
    "import sys\n",
    "sys.path.append(\"DREiMac\")\n",
    "from ProjectiveCoordinates import ProjCoords, getStereoProjCodim1\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define a few functions which will help us to sample patches from an image and to plot a collection of patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getPatches(I, dim):\n",
    "    \"\"\"\n",
    "    Given an image I, return all of the dim x dim patches in I\n",
    "    :param I: An M x N image\n",
    "    :param d: The dimension of the square patches\n",
    "    :returns P: An (M-d+1)x(N-d+1)x(d^2) array of all patches\n",
    "    \"\"\"\n",
    "    #http://stackoverflow.com/questions/13682604/slicing-a-numpy-image-array-into-blocks\n",
    "    shape = np.array(I.shape*2)\n",
    "    strides = np.array(I.strides*2)\n",
    "    W = np.asarray(dim)\n",
    "    shape[I.ndim:] = W\n",
    "    shape[:I.ndim] -= W - 1\n",
    "    if np.any(shape < 1):\n",
    "        raise ValueError('Window size %i is too large for image'%dim)\n",
    "    P = np.lib.stride_tricks.as_strided(I, shape=shape, strides=strides)\n",
    "    P = np.reshape(P, [P.shape[0]*P.shape[1], dim*dim])\n",
    "    return P\n",
    "\n",
    "def imscatter(X, P, dim, zoom=1):\n",
    "    \"\"\"\n",
    "    Plot patches in specified locations in R2\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X : ndarray (N, 2)\n",
    "        The positions of each patch in R2\n",
    "    P : ndarray (N, dim*dim)\n",
    "        An array of all of the patches\n",
    "    dim : int\n",
    "        The dimension of each patch\n",
    "    \n",
    "    \"\"\"\n",
    "    #https://stackoverflow.com/questions/22566284/matplotlib-how-to-plot-images-instead-of-points\n",
    "    ax = plt.gca()\n",
    "    for i in range(P.shape[0]):\n",
    "        patch = np.reshape(P[i, :], (dim, dim))\n",
    "        x, y = X[i, :]\n",
    "        im = OffsetImage(patch, zoom=zoom, cmap = 'gray')\n",
    "        ab = AnnotationBbox(im, (x, y), xycoords='data', frameon=False)\n",
    "        ax.add_artist(ab)\n",
    "    ax.update_datalim(X)\n",
    "    ax.autoscale()\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "\n",
    "def plotPatches(P, zoom = 1):\n",
    "    \"\"\"\n",
    "    Plot patches in a best fitting rectangular grid\n",
    "    \"\"\"\n",
    "    N = P.shape[0]\n",
    "    d = int(np.sqrt(P.shape[1]))\n",
    "    dgrid = int(np.ceil(np.sqrt(N)))\n",
    "    ex = np.arange(dgrid)\n",
    "    x, y = np.meshgrid(ex, ex)\n",
    "    X = np.zeros((N, 2))\n",
    "    X[:, 0] = x.flatten()[0:N]\n",
    "    X[:, 1] = y.flatten()[0:N]\n",
    "    imscatter(X, P, d, zoom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: Patches On A Disc\n",
    "\n",
    "First, we start off by sampling patches from an image representing a disc\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First create an image of a disc\n",
    "res = 50\n",
    "R = res/2\n",
    "[I, J] = np.meshgrid(np.arange(res) ,np.arange(res))\n",
    "Im = ((I-R)**2 + (J-R)**2) < (0.5*R*R)\n",
    "Im = 1.0*Im\n",
    "plt.imshow(Im, interpolation='none', cmap='gray')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, sample all unique 5x5 patches from this image, for a collection of patches which lives in 25 dimensional Euclidean space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dim = 5\n",
    "P = getPatches(Im, dim)\n",
    "#Remove redundant patches to cut down on computation time\n",
    "toKeep = [0]\n",
    "XSqr = np.sum(P**2, 1)\n",
    "D = XSqr[:, None] + XSqr[None, :] - 2*P.dot(P.T)\n",
    "for i in range(1, D.shape[0]):\n",
    "    if np.sum(D[i, 0:i] == 0) > 0:\n",
    "        continue\n",
    "    toKeep.append(i)\n",
    "P = P[np.array(toKeep), :]\n",
    "plt.figure(figsize=(8, 8))\n",
    "plotPatches(P, zoom=3)\n",
    "ax = plt.gca()\n",
    "ax.set_facecolor((0.7, 0.7, 0.7))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's compute persistence diagrams up to H2 for this collection of patches\n",
    "\n",
    "* Based on the diagrams, what shape do the patches concentrate on?  Can you arrange the patches on that shape?\n",
    "* What happens if you get rid of the constant all black or all white patches, and you normalize the rest of the patches to have unit norm?  What topological manifold is this?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "dgms = ripser(P, maxdim=2)['dgms']\n",
    "plot_dgms(dgms)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: Oriented Line Segments\n",
    "\n",
    "We now examine the collection of patches which hold oriented, blurry line segments that are varying distances from the center of the patch.  First, let's start by setting up the patches.  Below, the \"dim\" variable sets the patch resolution, and the \"sigma\" variable sets the blurriness (a larger sigma means blurrier line segments). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def getLinePatches(dim, NAngles, NOffsets, sigma):\n",
    "    N = NAngles*NOffsets\n",
    "    P = np.zeros((N, dim*dim))\n",
    "    thetas = np.linspace(0, np.pi, NAngles+1)[0:NAngles]\n",
    "    #ps = np.linspace(-0.5*np.sqrt(2), 0.5*np.sqrt(2), NOffsets)\n",
    "    ps = np.linspace(-1, 1, NOffsets)\n",
    "    idx = 0\n",
    "    [Y, X] = np.meshgrid(np.linspace(-0.5, 0.5, dim), np.linspace(-0.5, 0.5, dim))\n",
    "    for i in range(NAngles):\n",
    "        c = np.cos(thetas[i])\n",
    "        s = np.sin(thetas[i])\n",
    "        for j in range(NOffsets):\n",
    "            patch = X*c + Y*s + ps[j]\n",
    "            patch = np.exp(-patch**2/sigma**2)\n",
    "            P[idx, :] = patch.flatten()\n",
    "            idx += 1\n",
    "    return P\n",
    "\n",
    "P = getLinePatches(dim=10, NAngles = 16, NOffsets = 16, sigma=0.25)\n",
    "plt.figure(figsize=(8, 8))\n",
    "plotPatches(P, zoom=2)\n",
    "ax = plt.gca()\n",
    "ax.set_facecolor((0.7, 0.7, 0.7))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compute persistence diagrams for this collection of patches.  This time, we will compute with both $\\mathbb{Z}/2$ coefficients and $\\mathbb{Z}/3$ coefficients up to H2.\n",
    "\n",
    "* Based on the persistence diagrams, what shape do the patches appear to concentrate on?  Can you arrange the patches on this shape to explain why?\n",
    "* What happens to the persistence diagrams when you make sigma very small and the patches become sharper, or when you make sigma close to 1 and the patches become very blurry?  Can you explain what's happening geometrically?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dgmsz2 = ripser(P, coeff=2, maxdim=2)['dgms']\n",
    "dgmsz3 = ripser(P, coeff=3, maxdim=2)['dgms']\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(121)\n",
    "plot_dgms(dgmsz2)\n",
    "plt.title(\"$\\mathbb{Z}/2$\")\n",
    "plt.subplot(122)\n",
    "plot_dgms(dgmsz3)\n",
    "plt.title(\"$\\mathbb{Z}/3$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will look at these patches using \"projective coordinates\" (finding a map to $RP^2$).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def plotProjBoundary():\n",
    "    t = np.linspace(0, 2*np.pi, 200)\n",
    "    plt.plot(np.cos(t), np.sin(t), 'c')\n",
    "    plt.axis('equal')\n",
    "    ax = plt.gca()\n",
    "    ax.arrow(-0.1, 1, 0.001, 0, head_width = 0.15, head_length = 0.2, fc = 'c', ec = 'c', width = 0)\n",
    "    ax.arrow(0.1, -1, -0.001, 0, head_width = 0.15, head_length = 0.2, fc = 'c', ec = 'c', width = 0)\n",
    "    ax.set_facecolor((0.35, 0.35, 0.35))\n",
    "\n",
    "P = getLinePatches(dim=10, NAngles = 200, NOffsets = 200, sigma=0.25)\n",
    "res = ProjCoords(P, n_landmarks=100, proj_dim=2, perc=0.9)\n",
    "X = res['X']\n",
    "idx = getGreedyPerm(X, 400)['perm']\n",
    "SFinal = getStereoProjCodim1(X[idx, :], randomSeed=-1)\n",
    "P = P[idx, :]\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "imscatter(SFinal, P, 10)\n",
    "plotProjBoundary()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Natural Image Patches\n",
    "\n",
    "We will now generate a set of patches that occur in \"natural images,\" which are essentially gradients from dark to light in different directions, which are centered at the patch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def getNaturalImagePatches(ndir, nsharp, dim):\n",
    "    N = ndir*nsharp\n",
    "    t = np.linspace(0, 2*np.pi, ndir+1)[0:ndir]\n",
    "    a, b = np.cos(t), np.sin(t)\n",
    "    t = np.linspace(0, 2*np.pi, nsharp+1)[0:nsharp]\n",
    "    c, d = np.cos(t), np.sin(t)\n",
    "    a, b, c, d = a.flatten(), b.flatten(), c.flatten(), d.flatten()\n",
    "    hdim = int((dim-1)/2)\n",
    "    xr = np.linspace(-1, 1, dim)\n",
    "    X, Y = np.meshgrid(xr, xr)\n",
    "    P = np.zeros((N, dim*dim))\n",
    "    idx = 0\n",
    "    for i in range(a.size):\n",
    "        for j in range(c.size):\n",
    "            proj = a[i]*X + b[i]*Y\n",
    "            p = c[j]*proj + d[j]*(proj**2)\n",
    "            P[idx, :] = p.flatten()\n",
    "            idx += 1\n",
    "    return P\n",
    "\n",
    "res = 15\n",
    "dim = 8\n",
    "P = getNaturalImagePatches(res, res, dim)\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plotPatches(P, zoom = 2)\n",
    "ax = plt.gca()\n",
    "ax.set_facecolor((0.15, 0.15, 0.15))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the persistent homology of this collection of patches with $\\mathbb{Z} / 2\\mathbb{Z}$ $\\mathbb{Z} / 3\\mathbb{Z}$ coefficients\n",
    "\n",
    "#### Questions\n",
    "* What topological manifold do these patches concentrate on, based on what you see in the persistence diagrams (hint: be careful that two points may be on top of each other in H1 for $\\mathbb{Z} / 2\\mathbb{Z}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 200\n",
    "dgmsz2 = ripser(P, coeff=2, maxdim=2)['dgms']\n",
    "dgmsz3 = ripser(P, coeff=3, maxdim=2)['dgms']\n",
    "\n",
    "I1 = dgmsz2[1]\n",
    "I2 = dgmsz3[1]\n",
    "I1 = I1[np.argsort(I1[:, 0]-I1[:, 1]), :]\n",
    "I2 = I2[np.argsort(I2[:, 0]-I2[:, 1]), :]\n",
    "print(\"Max 2 for Z2:\\n%s\"%I1[0:2, :])\n",
    "print(\"\\nMax 2 for Z3:\\n%s\"%I2[0:2, :])\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(121)\n",
    "plot_dgms(dgmsz2)\n",
    "plt.title(\"$\\mathbb{Z}/2$\")\n",
    "plt.subplot(122)\n",
    "plot_dgms(dgmsz3)\n",
    "plt.title(\"$\\mathbb{Z}/3$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at these points with projective coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = 500\n",
    "dim = 8\n",
    "P = getNaturalImagePatches(res, res, dim)\n",
    "\n",
    "res = ProjCoords(P, n_landmarks=100, proj_dim=2, perc=0.9)\n",
    "X = res['X']\n",
    "idx = getGreedyPerm(X, 400)['perm']\n",
    "SFinal = getStereoProjCodim1(X[idx, :], randomSeed=-1)\n",
    "P = P[idx, :]\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.subplot(121)\n",
    "plot_dgms(res['dgm1'], labels=['H1'])\n",
    "plt.subplot(122)\n",
    "imscatter(SFinal, P, dim)\n",
    "plotProjBoundary()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}