cd_check_poly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CD Test (Polygons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Template for running CD algorithms through parameterizations. For these runs, focus on parameterizing CD algorithms, see at what points affine offset and sun azimuth directions begin to make TPR drop. \n",
    "\n",
    "\n",
    "Steps: \n",
    "1. generate base DEM\n",
    "2. create after DEM by filling depressions (use defaults for now)\n",
    "3. add craters to after image (use defaults for now)\n",
    "4. hillshade both images (use azi param to control light direction)\n",
    "5. apply noise to both images\n",
    "6. apply affine to both images and use after affine on truth data (stick with offset differences for now, as it tends to be the most common registration error)\n",
    "7. run before and after through a CD algorithm \n",
    "8. compute TPR and FPR \n",
    "    * TPR = TP/TP+FP\n",
    "    * FPR = FP/FP+TP\n",
    "    * IF function returns points, TP = point equals to or is contained by set of \"truth\" points/polys 1 pixel buffer. \n",
    "    * ELSE function returns polygons, TP = polygon has >50% intersection with \"truth\" poly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/setuptools/distutils_patch.py:26: UserWarning: Distutils was imported before Setuptools. This usage is discouraged and may exhibit undesirable behaviors or errors. Please use Setuptools' objects directly or at least import Setuptools first.\n",
      "  \"Distutils was imported before Setuptools. This usage is discouraged \"\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "os.environ['ISISROOT'] = '/usgs/cpkgs/anaconda3_linux/envs/isis3.9.0'\n",
    "\n",
    "from numpy import random\n",
    "from scipy.ndimage.interpolation import zoom\n",
    "import scipy\n",
    "\n",
    "import os\n",
    "import math\n",
    "\n",
    "from PIL import Image\n",
    "import numpy as np\n",
    "from scipy.misc import imresize\n",
    "\n",
    "import shapely\n",
    "from shapely.geometry import Point, MultiPolygon\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "\n",
    "from autocnet.utils.utils import bytescale, cartesian\n",
    "from autocnet.examples import get_path\n",
    "from autocnet.cg import change_detection as cd\n",
    "from autocnet.cg import cg\n",
    "\n",
    "def cartesian_product(*arrays):\n",
    "    la = len(arrays)\n",
    "    dtype = \"object\"\n",
    "    arr = np.empty([len(a) for a in arrays] + [la], dtype=dtype)\n",
    "    for i, a in enumerate(np.ix_(*arrays)):\n",
    "        arr[...,i] = a\n",
    "    return arr.reshape(-1, la)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate before and after images "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)\n",
      "Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)\n"
     ]
    }
   ],
   "source": [
    "# Generate base DEM \n",
    "\n",
    "before_dem = cd.generate_dem(alpha=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "before_dem, before_polys, after_dem, after_polys = cd.generate_boulder_field(before_dem, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(MultiPolygon(after_polys).geoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)\n",
      "Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)\n",
      "Warning! No geotransform defined. Choosing a standard one! (Top left cell's top let corner at <0,0>; cells are 1x1.)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/krodriguez/repos/autocnet/autocnet/cg/change_detection.py:614: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  dmask[[*cluster.T]] = True\n"
     ]
    }
   ],
   "source": [
    "# Compute Depressions\n",
    "\n",
    "after_dem, change_poly = cd.compute_depression(after_dem, alpha=0.5, return_polygon=True)\n",
    "change_poly = shapely.ops.unary_union([change_poly, MultiPolygon(after_polys)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hillshade \n",
    "azi_before = 190 \n",
    "azi_after = 210\n",
    "\n",
    "before_image = cd.hillshade(before_dem, azi=azi_before).mean(axis=2)\n",
    "after_image = cd.hillshade(after_dem, azi=azi_after).mean(axis=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Add noise\n",
    "noise_stddev_before = np.std(before_image)/10\n",
    "noise_stddev_after = np.std(after_image)/10\n",
    "\n",
    "before_image += np.random.normal(0, noise_stddev_before, before_image.shape)\n",
    "after_image += np.random.normal(0, noise_stddev_after, after_image.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# TODO: Add affine distortion "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = (after_image-before_image).astype(float)\n",
    "\n",
    "fig, axs = plt.subplots(2,2, figsize=(10,10))\n",
    "\n",
    "axs[0][0].imshow(before_image, cmap=\"bone\")\n",
    "axs[0][0].set_title('\"Before\" Hillshade')\n",
    "axs[0][1].imshow(after_image, cmap=\"bone\")\n",
    "axs[0][1].set_title('\"After\" Hillshade')\n",
    "\n",
    "axs[1][0].imshow(after_image, cmap='bone')\n",
    "\n",
    "if isinstance(change_poly, MultiPolygon):\n",
    "    for geom in change_poly.geoms:    \n",
    "        xs, ys = geom.exterior.xy    \n",
    "        axs[1][0].fill(ys, xs, alpha=1, fc='black', ec='none')\n",
    "else: \n",
    "    xs, ys = change_poly.exterior.xy    \n",
    "    axs[1][0].fill(ys, xs, alpha=1, fc='black', ec='none')\n",
    "axs[1][0].set_title(\"Change Polygons\")\n",
    "    \n",
    "axs[1][1].imshow(diff, cmap=\"magma\")\n",
    "axs[1][1].set_title('Diff')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate CD parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9X3 Matrix\n"
     ]
    }
   ],
   "source": [
    "# Change Based on Algorithm\n",
    "\n",
    "nfeatures = [500, 1000, 2000]\n",
    "scalefactor = [1.5]\n",
    "nlevels = [3]\n",
    "extractor_params = [dict(zip([\"nfeatures\", \"scaleFactor\", \"nlevels\"], ps)) for ps in cartesian_product(nfeatures, scalefactor, nlevels)] \n",
    "\n",
    "min_samples = [5, 15, 20]\n",
    "max_eps = [1, 20, 30]\n",
    "xi = [.5]\n",
    "eps = [.5]\n",
    "cluster_params =  [dict(zip([\"min_samples\", \"max_eps\", \"eps\", \"xi\"], ps)) for ps in cartesian_product(min_samples, max_eps, eps, xi)] \n",
    "\n",
    "print(f\"{len(cluster_params)}X{len(extractor_params)} Matrix\")\n",
    "parameters = cartesian_product(extractor_params, cluster_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run OKB with Params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/work/users/krodriguez/anaconda3/envs/autocnetdev/lib/python3.7/site-packages/sklearn/cluster/_optics.py:503: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n"
     ]
    }
   ],
   "source": [
    "results = [cd.okbm_detector(before_image, \n",
    "                   after_image, \n",
    "                   extractor_kwargs=params[0], cluster_kwargs=params[1]) for params in parameters]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "126.5625"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results[5][0].geometry[2].area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Polygonal Results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# OKB, needs to be within 1 pixel of truth data \n",
    "\n",
    "nx = len(extractor_params)\n",
    "ny = len(cluster_params)\n",
    "\n",
    "fig, axs = plt.subplots(ny, nx, sharex='col', sharey='row')\n",
    "fig.set_figheight(50)\n",
    "fig.set_figwidth(20)\n",
    "\n",
    "for ax in axs.flat:\n",
    "    ax.set(xlabel='x-label', ylabel='y-label')\n",
    "    \n",
    "xs, ys = np.mgrid[:ny, :nx]\n",
    "cmap = matplotlib.cm.get_cmap('plasma')\n",
    "\n",
    "for i, coord in enumerate(zip(xs.flatten(),ys.flatten())):\n",
    "    x,y = coord\n",
    "    \n",
    "    # unpack results \n",
    "    df, diff = results[i]\n",
    "    \n",
    "    # geopandas likes to swap x,y\n",
    "    sdf = df.geometry.map(lambda polygon: shapely.ops.transform(lambda x, y: (y, x), polygon))\n",
    "    \n",
    "    # get true positives and false negatives\n",
    "    true_positives = 0\n",
    "    for geom in sdf.geometry:\n",
    "        for truth_geom in change_poly.geoms:\n",
    "            if geom.intersection(truth_geom).area >= geom.area/5:\n",
    "                true_positives += 1\n",
    "                break\n",
    "        \n",
    "    false_positives = df.shape[0] - true_positives\n",
    "            \n",
    "    ratio = true_positives/df.shape[0] if df.shape[0] is not 0 else np.nan\n",
    "    \n",
    "    axs[x,y].set_title(f\"{i} | tp: {true_positives}, fp: {false_positives}, score: {ratio}\")\n",
    "    axs[x,y].set(xlabel=f\"ex {parameters[i][1]}\", ylabel=f\"cl {parameters[i][0]}\")\n",
    "    \n",
    "    axs[x,y].imshow(after_image, alpha=0.2, cmap=\"bone\")\n",
    "    for geom in change_poly.geoms:    \n",
    "        xs, ys = geom.exterior.xy    \n",
    "        axs[x,y].fill(ys, xs, alpha=0.5, fc='black', ec='none')\n",
    "    \n",
    "    df.plot(ax=axs[x,y])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "AutocnetDev",
   "language": "python",
   "name": "autocnetdev"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}