{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "matplotlib.rcParams['image.cmap'] = 'gray'\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Going interplanetary\n",
    "\n",
    "The sun is one of the most spherical objects in our solar system.\n",
    "According to an [article in Scientific American](http://www.scientificamerican.com/gallery/well-rounded-sun-stays-nearly-spherical-even-when-it-freaks-out/):\n",
    "\n",
    "> Earth's closest star is one of the roundest objects humans have\n",
    "> measured. If you shrank the sun down to beach ball size, the\n",
    "> difference between its north-south and the east-west diameters would\n",
    "> be thinner than the width of a human hair, says Jeffery Kuhn, a\n",
    "> physicist and solar researcher at the University of Hawaii at\n",
    "> Manoa. \"Not only is it very round, but it's too round,\" he adds. The\n",
    "> sun is more spherical and more invariable than theories predict.\n",
    "\n",
    "If the sun is spherical, we should be able to fit a circle to a 2D\n",
    "slice of it!  Your task is to do just that, using RANSAC and scikit-image's CircleModel.\n",
    "\n",
    "Let's start by loading an example image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import io\n",
    "\n",
    "image = io.imread('../../images/superprom_prev.jpg')\n",
    "\n",
    "f, ax = plt.subplots(figsize=(8, 8))\n",
    "ax.imshow(image);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this specific image, we got a bit more than we bargained for in the\n",
    "form of magnificently large solar flares.  Let's see if some *canny\n",
    "edge detection* will help isolate the sun's boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import feature, color\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "edges = feature.canny(color.rgb2gray(image), sigma=2)\n",
    "ax.imshow(edges, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The edges look good, but there's a lot going on inside the sun.  We\n",
    "use RANSAC to fit a robust circle model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from skimage.measure import ransac, CircleModel\n",
    "\n",
    "points = np.array(np.nonzero(edges)).T\n",
    "\n",
    "model_robust, inliers = ransac(points, CircleModel, min_samples=3,\n",
    "                               residual_threshold=2, max_trials=5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The parameters of the circle are center x, y and radius:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_robust.params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the results, drawing a circle on the sun, and also\n",
    "highlighting inlier vs outlier edge pixels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage import draw\n",
    "\n",
    "cy, cx, r = model_robust.params\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(15, 8))\n",
    "\n",
    "ax0.imshow(image)\n",
    "ax1.imshow(image)\n",
    "\n",
    "ax1.plot(points[inliers, 1], points[inliers, 0], 'b.', markersize=1)\n",
    "ax1.plot(points[~inliers, 1], points[~inliers, 0], 'g.', markersize=1)\n",
    "ax1.axis('image')\n",
    "\n",
    "circle = plt.Circle((cx, cy), radius=r, facecolor='none', linewidth=2)\n",
    "ax0.add_patch(circle);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circular fit is, indeed, excellent, and rejects all the inner\n",
    "squiggly edges generated by solar turbulence!\n",
    "\n",
    "Note a general principle here: algorithms that aggregate across an\n",
    "entire path are often robust against noise.  Here, we have *high\n",
    "uncertainty* in the solar edge, but also know that only the solar edge\n",
    "pixels contribute coherently to the full circular path around the\n",
    "solar edge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: CardShark\n",
    "\n",
    "Your small start-up, CardShark, run from your garage over nights and\n",
    "evenings, takes photos of credit cards and turns them into machine\n",
    "readable information.\n",
    "\n",
    "The first step is to identify where in a photo the credit card is\n",
    "located.\n",
    "\n",
    "1. Load the photo `../../images/credit_card.jpg`\n",
    "2. Using RANSAC and LineModelND shown above, find the first most\n",
    "   prominent edge of the card\n",
    "3. Remove the datapoints belonging to the most prominent edge, and\n",
    "   repeat the process to find the second, third, and fourth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "\n",
    "image = io.imread('../../images/credit_card.jpg')\n",
    "ax.imshow(image);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.measure import LineModelND\n",
    "\n",
    "f, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "edges = feature.canny(color.rgb2gray(image), sigma=3)\n",
    "edge_pts = np.array(np.nonzero(edges), dtype=float).T\n",
    "edge_pts_xy = edge_pts[:, ::-1]\n",
    "\n",
    "for i in range(4):\n",
    "    model_robust, inliers = ransac(edge_pts_xy, LineModelND, min_samples=2,\n",
    "                                   residual_threshold=1, max_trials=1000)\n",
    "\n",
    "    x = np.arange(800)\n",
    "    plt.plot(x, model_robust.predict_y(x))\n",
    "\n",
    "    edge_pts_xy = edge_pts_xy[~inliers]\n",
    "\n",
    "plt.imshow(edges);"
   ]
  }
 ],
 "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}