homographies.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image Stitching:  The final step\n",
    "Reading: Szeliski 6.1, 9.1\n",
    "\n",
    "## Homographies\n",
    "\n",
    "We have successfully made our way through all the steps of image stitching save one.  We have detected keypoints from two images (called image 1 and image 2), generated descriptors, matched them together, and now all we need to do is to determine the way in which we need to move one image (let's call it image 2) such that its pixel intensity information is in *the same coordinate system as image 1*.  \n",
    "\n",
    "As it turns out, under a very specific assumption about the motion of the camera between images, this transformation is a linear operation called a *homography*.  To see exactly what this operation does, let's reference the following image:\n",
    "\n",
    "<img src=\"homography.jpg\">\n",
    "\n",
    "In this image, the blue points are *things* out in the real world, perhaps ones that we've identified as corners and come up with descriptors for.  We'll assume that we've already solved the problem of correspondence: if we see a corner in both images, we have successfully matched it.  These corners reflect light and send out a ray (blue line) which impacts the effective focal plane of the camera.  In the figure, these impacts are denoted by an X in a color corresponding to the image whose focal plane was impacted.  What we need to do is to determine where corners that impacted the second image would have impacted the focal plane of the first image *if that first image had its sensor synthetically extended*, i.e. what are the coordinates of that black question mark?  We can determine this mapping by looking at the points for which we have correspondences.\n",
    "\n",
    "As it turns out, the transformation from one 2D plane embedded in 3D to another obeys a linear transformation called a *homography*.  A homography is a linear transformation that preserves straight lines (and nothing else).  It occupies the most expressive position on a continuum that ranges from simple rotation (governed by one parameter), to rotation-translations (which preserve scale and relative orientation and position of points), to affine transformations (which preserve scale and parallelity), to homographies (which only preserve straight lines).  We can define it as\n",
    "$$\n",
    "\\begin{bmatrix} x' \\\\ y' \\\\ w' \\end{bmatrix} = \\begin{bmatrix} h_{11} & h_{12} & h_{13}\\\\\n",
    "                                                               h_{21} & h_{22} & h_{23}\\\\\n",
    "                                                               h_{31} & h_{32} & h_{33} \\end{bmatrix} \n",
    "                                                               \\begin{bmatrix} u \\\\\n",
    "                                                                               v \\\\\n",
    "                                                                               1 \\end{bmatrix}\n",
    "$$                                                               \n",
    "Note that this makes use of generalized coordinates: after the transformation, we need to divide: $u' = x'/w'$ and $v' = y'/w'$.  What does a homography do?  Here is an example of a random homography on the unit square:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGx5JREFUeJzt3Xt01dWd9/H3NwmIlIsiGWG4GKpQS72ABsRLBQeVA+3A1FtBEakIBmufdpzpLB1djrXTNVO7nsdZzrIKIkUdFa1jK1oUfVS0o6LGEVSoWh5Ga7xARKVYVAS+zx/7WGNIOL8k55z9yy+f11pZOZffOef7W0k+bPbev73N3RERkWypiF2AiIgUn8JdRCSDFO4iIhmkcBcRySCFu4hIBincRUQySOEuIpJBCncRkQxSuIuIZFBVrA/u37+/19TUxPp4EZFO6bnnnnvX3asLHRct3Gtqaqivr4/18SIinZKZvZ7kOHXLiIhkkMJdRCSDFO4iIhmkcBcRySCFu4hIBhUMdzNbbGabzOylVp43M7vGzNab2QtmdkTxyxQRkbZI0nJfAuT28PxkYHj+ax5wXcfLEhGRjig4z93dHzezmj0cMg242cN+favMbB8zG+jubxepxi9Y+NxCbnvxtlK8tYhEcuahZzLvyHmxy8iUYvS5DwLeaHK/If/YbsxsnpnVm1l9Y2Njuz7sthdvY/U7q9v1WhFJn9XvrFaDrQSKcYWqtfBYi7tuu/tCYCFAbW1tu3fmHjVgFCtnr2zvy0UkRSYsmRC7hEwqRsu9ARjS5P5g4K0ivK+IiLRTMcJ9GTArP2tmHLClVP3tIiKSTMFuGTO7HZgA9DezBuCfgG4A7n49sByYAqwHtgHfKVWxIiKSTJLZMjMKPO/Ad4tWkYiIdJiuUBURySCFu4hIBincRUQySOEuIpJBCncRkQxSuIuIZJDCXUQkgxTuIiIZpHAXEckghbuISAYp3EVEMkjhLiKSQQp3EZEMUriLiGSQwl1EJIMU7iIiGaRwFxHJIIW7iEgGKdxFRDJI4S4ikkEKdxGRDFK4i4hkUFXsAkSkc9i4ES67DLp3h/79P/+qrv7i/R49YlcqoHAXkYRefRUWLYK99oLt28G95eO+9KUvhv2e/iHo3x927aikompneU+mC1C4i0giY8eG4J41C/793+H99+Hdd8NXY+Pnt5t//f734fsf/9jaOz9MVc+tDP9Jsn8MBg0KdcieKdxFJJG99oKJE+H++6Gi4vOwTWr7dti8efd/DK56aDE7PuzLmOpTefddePNNWLMmHPPxx7u/T79+sGED9O1bvHPLIoW7iCSWy8GyZaGL5itfadtru3eHgQPDV1N39rwZgNtmn7rba7Zt++L/Dl5+GX7wA7j1VrjggvaeRdeg2TIiktikSeH7Aw+U5/N69oShQ+GII8Jnf//7cOSRcP31rff5S6BwF5HEvvxlGDGifOHekro6ePFFeOqpeDV0Bgp3EWmTXA5WroSPPorz+dOnQ58+cN11cT6/s0gU7maWM7NXzGy9mV3cwvNDzexRM3vezF4wsynFL1VE0iCXCwOdjz8e5/N79YKzz4Zf/jL0xUvLCoa7mVUC1wKTgZHADDMb2eywy4A73X00MB34ebELFZF0GD8+zJyJ3TXzySdw003xaki7JC33scB6d9/g7tuBpcC0Zsc40Cd/uy/wVvFKFJE06dkzBHzMcD/kEDjuOFiwAHbtildHmiUJ90HAG03uN+Qfa+oKYKaZNQDLge8VpToRSaVcLkxLfO21eDXU1YULpB55JF4NaZYk3K2Fx5pPQpoBLHH3wcAU4BYz2+29zWyemdWbWX1jY2PbqxWRVMjlwvcVK+LVcOqpsN9+YVqk7C5JuDcAQ5rcH8zu3S5zgDsB3P0poAew27Vr7r7Q3Wvdvba6urp9FYtIdAcfHOafx+ya6dEDvvMd+PWv4S11BO8mSbg/Cww3s2Fm1p0wYLqs2TF/ACYCmNlXCeGuprlIRpmF1vvDD4dlBWI5/3zYuRMWL45XQ1oVDHd33wFcCKwAfkeYFbPWzK40s6n5w/4OmGtma4Dbgdnuun5MJMtyOdi6Ne7FRAcdBCedBAsXhpCXzyVaW8bdlxMGSps+dnmT2+uAY4tbmoik2cSJUFUVumbGj49XR11d6H9fvhz++q/j1ZE2ukJVRNqlTx849ti4/e4QAn3gQA2sNqdwF5F2y+Vg9Wp4++14NXTrBnPnhqWIY07NTBuFu4i0WxqmRAKcd14Y5L3hhrh1pInCXUTa7fDDYcCA+F0zQ4bAN78ZtgGMOXsnTRTuItJuZmGd9QcfjD9bpa4ONm0K895F4S4iHZTLhf1Un302bh2TJkFNjQZWP6NwF5EOOemk0IKP3TVTUREuanr00bDuTVencBeRDtlvPxg7Nn64Q1iOoFu3sFpkV6dwF5EOy+XgmWdg8+a4dey/P5xyCixZEm+nqLRQuItIh+VyYcPqhx6KXUkYWP3gA7jzztiVxKVwF5EOGzMG9t03HV0z48eHVSu7+sCqwl1EOqyyEk4+OYR77J2RzELrfdWqcPVsV6VwF5GiyOVg40Z44YXYlcCsWWG9967cele4i0hRTJoUvqeha2bffWH6dLj11rAscVekcBeRohg4EEaNSke4A8yfDx9+GAK+K1K4i0jR5HLwxBPwxz/GriQM8o4eHbpmuuLWQQp3ESmaXA527IBHHoldyecDq2vWhMHVrkbhLiJFc/TR0Lt3WFs9Dc48M9TTFQdWFe4iUjTdu4ft9x54IB1dIb16wdlnwx13wHvvxa6mvBTuIlJUuRz84Q/pWbzr/PPhk0/gpptiV1JeCncRKao0TYkEOOwwOOaYrjewqnAXkaKqqQmX/6cl3CEMrL76algOuKtQuItI0eVy8NhjsG1b7EqC00+Hfv261sCqwl1Eii6XC/3cjz0Wu5KgR4+w1vuvfgXvvBO7mvJQuItI0R1/fAjUNHXNzJsX5uAvXhy7kvJQuItI0e29N0yYkK5wHzEiTNNcsCD+Zt7loHAXkZLI5cIg5oYNsSv5XF1dmKaZpn90SkXhLiIlkcuF7ytWxK2jqWnTYMCArjGwqnAXkZIYMQKGDUtXK7lbNzjvPPjNb+D112NXU1oKdxEpCbPQen/4Ydi+PXY1n5s7N9S2aFHsSkpL4S4iJZPLwZ/+FJYBTouhQ2HKlBDun34au5rSSRTuZpYzs1fMbL2ZXdzKMWeY2TozW2tmtxW3TBHpjE44IXSFpKlrBsJGHu+8A/fcE7uS0ikY7mZWCVwLTAZGAjPMbGSzY4YDlwDHuvvXgB+UoFYR6WR694bjjktfuE+aBAcckO2B1SQt97HAenff4O7bgaXAtGbHzAWudff3Adx9U3HLFJHOKpcLm2a/+WbsSj5XWRkuanr44TBdM4uShPsg4I0m9xvyjzU1AhhhZk+Y2Sozy7X0RmY2z8zqzay+sbGxfRWLSKeSximRAOeeC1VV4aKmLEoS7tbCY80XzqwChgMTgBnAIjPbZ7cXuS9091p3r62urm5rrSLSCR16aNg8O21dMwMGwLe+BUuWwEcfxa6m+JKEewMwpMn9wcBbLRxzj7t/6u7/A7xCCHsR6eI+mxL50ENhbZc0mT8/7NB0112xKym+JOH+LDDczIaZWXdgOrCs2TG/Bk4AMLP+hG6aFF10LCIx5XLwwQfwzDOxK/miCRPCxVZZHFgtGO7uvgO4EFgB/A64093XmtmVZjY1f9gKYLOZrQMeBX7o7ptLVbSIdC4nnggVFenrmjEL6808+WQY9M2SRPPc3X25u49w9wPd/Sf5xy5392X52+7uF7n7SHc/1N2XlrJoEelc+vWDo45KX7gDnHNOWJ44a613XaEqImWRy0F9PaRtoly/fvDtb8Mtt8DWrbGrKR6Fu4iURS4XNqh+6KHYleyurg4+/BBuvz12JcWjcBeRsqithf7909k1c9RRcPjhoWvGm0/07qQU7iJSFhUVcPLJ4WKmXbtiV/NFnw2sPv98+mb0tJfCXUTKJpeDTZtg9erYlezurLOgV6/sDKwq3EWkbE4+OXxPY9dM794wcyYsXQrvvx+7mo5TuItI2ey/PxxxRDrDHULXzMcfw803x66k4xTuIlJWuVy4aGjLltiV7O7ww2HcuGwMrCrcRaSscjnYuTMst5tGdXXw8svw2GOxK+kYhbuIlNW4cdCnD9x/f+xKWnbGGbDvvp1/YFXhLiJl1a1bWGvmgQfS2fWx994wezbcfTds3Bi7mvZTuItI2eVy0NAA69bFrqRl558fNs/+xS9iV9J+CncRKbtJk8L3tM6a+cpXwubeCxaE8YHOSOEuImU3dCiMHJnecIewkcdrr8GDD8aupH0U7iISRS4Hjz8OOz/pEbuUFk2bFubld9aBVYW7iESRy8H27fDBy6Nil9Ki7t1hzhy47z54443Y1bSdwl1Eovj616FnT3jvxbGxS2nV3LlhRs+iRbEraTuFu4hE0aNHGLR8P8XhXlMDU6bADTeE2TOdicJdRKLJ5eCjTYP5aOOg2KW0qq4O3n4b7r03diVto3AXkWhyufD9vZfGxC1kDyZPhiFDOt/AqsJdRKI56CDo8Rdv8t6acbFLaVVlJcybF7YHXL8+djXJKdxFJKrq2sd478VxnHcebNsWu5qWzZkDVVXhoqbOQuEuIlENO+VGhn7zFhYvhrFjYe3a2BXtbuBA+Ju/CcsRfPxx7GqSUbiLSFRWuZNhp97IAw9AYyOMGRNCNG2LitXVwebN8J//GbuSZBTuIpIKJ58c9lYdNw7OPRdmzYKtW2NX9bkTToDhwzvPwKrCXURSY+DAMHD5ox/BbbdBbW16NtOuqAirRf7Xf8FLL8WupjCFu4ikSmUlXH552Klp69bQkr/uunR008yeDXvt1Tla7wp3EUmlCRNgzZrQHXLBBfDtb8ffd3W//cJOTTffDB9+GLeWQhTuIpJa1dXwm9/AT38adkYaPRqefTZuTXV14X8US5fGraMQhbuIpFpFBfzDP+SXB94Jxx4L//Zv8bppjj4aDj00/V0zicLdzHJm9oqZrTezi/dw3Glm5mZWW7wSRUTgmGPg+efDQl5/+7dh3vl775W/DrOwkcdzz8X/X8SeFAx3M6sErgUmAyOBGWY2soXjegP/C3i62EWKiAD06we/+lVoud9/P4waBU88Uf46zjoLvvSldLfek7TcxwLr3X2Du28HlgLTWjjux8BVQCe5fktEOiMz+P734cknoVs3GD8e/vVfYdeu8tXQp08I+Ntvhw8+KN/ntkWScB8ENN2HpCH/2J+Z2WhgiLvfV8TaRERaVVsL//3fcOqpcMklobtm06byfX5dHXz0EdxyS/k+sy2ShLu18NifhzLMrAK4Gvi7gm9kNs/M6s2svrGxMXmVIiIt6Ns3zFq5/npYuTJ006xcWZ7PHj06rIVz/fXpmIPfXJJwbwCGNLk/GHiryf3ewCHASjN7DRgHLGtpUNXdF7p7rbvXVldXt79qEZE8s3Dl6NNPh+6SiRPDFa47d5b+s+fPh3Xr4Le/Lf1ntVWScH8WGG5mw8ysOzAdWPbZk+6+xd37u3uNu9cAq4Cp7l5fkopFRFpw+OFQXx/6wq+4Ak46Cd56q+DLOuSMM2CffdI5sFow3N19B3AhsAL4HXCnu681syvNbGqpCxQRSapXr3D16C9+EVryo0bBihWl+7yePeGcc+Cuu8rb359Eonnu7r7c3Ue4+4Hu/pP8Y5e7+7IWjp2gVruIxDR7dpiDvv/+YSu/Sy6BHTtK81nnnx82z16ypDTv3166QlVEMmnkyNB6nzs3TJWcMAHeeKPgy9rsq18N771gQXmnYxaicBeRzOrZExYuDMsHr1kTumnuvbf4n1NXBxs2hOWK00LhLiKZN2NGmBN/wAEwdSpcdBFs31689//Wt8IiZ2kaWFW4i0iXMHw4PPUUfO97cPXVcNxxobVdDN27h020770XGhqK854dpXAXkS5jr73gmmvCPqi//324EOmuu4rz3vPmhT73RYuK834dpXAXkS7nlFPCCpMHHwynnx42A/m4g6tiDRsWZubccEPpZua0hcJdRLqkmppwZenf/33Yxm/cOHj11Y69Z11duHDqvhSssqVwF5Euq3t3+NnPQhg3NMARR8Ctt7b//aZMgcGD0zGwqnAXkS7vG9+A1atDH/zMmWFwdNu2tr9PVVWYV79iRfEGa9tL4S4iQmhxP/ooXHppWL5gzBhYu7bt73PeeVBZGS5qiknhLiKSV1UF//zPoeX97rsh4BcvbtuSvn/5lzBtWnjdJ5+UrtZCFO4iIs2cdFK4ovWYY0IXzdlnw9atyV9fVxf+cbj77tLVWIjCXUSkBQMGhBb8j38cttM78sjQL5/ExIlw4IFxB1YV7iIiraishMsuC33xf/pTmC75858X7qapqAit98cfb1+/fTEo3EVECjj++NBq/6u/gu9+N1z4VGhj7Nmzw1TLWAOrCncRkQSqq8N8+KuugnvuCdMmn3mm9eP79w//CNx8c2j1l5vCXUQkoYoK+OEPQ3fLrl1h8bGrr269m6auDrZsgTvuKG+doHAXEWmzo48O3TTf+EZYPnjqVNi8effjjj0Wvva1OAOrCncRkXbYd98w1fGaa+DBB8NGIE888cVjzGD+/LDl33PPlbc+hbuISDuZhfXhn3wyLCc8fjz8y798cbu9mTPDjlDlbr0r3EVEOujII8NOT6edBv/4jzB5MmzcGJ7r2xfOPDNs9bdlS/lqUriLiBRBnz7hYqcFC8KA66hR8Mgj4bm6urAQ2X/8R/nqUbiLiBSJWdiR6emnYZ994MQT4YorQtCPGRPWjW/LOjUdoXAXESmyww4Lg6izZsGPfhRCfurUcLVq80HXUlG4i4iUQK9esGRJ+HrmmTDQCuUbWFW4i4iU0DnnhGmQBx4Y7t96K7z9duk/V+EuIlJiBx8c+uG//vVw/7DDSv+ZCncRkTLYe+8wi+bcc2HdutJ/XlXpP0JERD5z443l+Ry13EVEMkjhLiKSQYnC3cxyZvaKma03s4tbeP4iM1tnZi+Y2cNmdkDxSxURkaQKhruZVQLXApOBkcAMMxvZ7LDngVp3Pwy4C7iq2IWKiEhySVruY4H17r7B3bcDS4FpTQ9w90fdfVv+7ipgcHHLFBGRtkgS7oOAN5rcb8g/1po5wP0tPWFm88ys3szqGxsbk1cpIiJtkiTcrYXHWlz6xsxmArXAz1p63t0Xunutu9dWV1cnr1JERNokyTz3BmBIk/uDgbeaH2RmJwKXAuPd/ZPilCciIu2RpOX+LDDczIaZWXdgOrCs6QFmNhpYAEx1903FL1NERNqiYLi7+w7gQmAF8DvgTndfa2ZXmtnU/GE/A3oBvzSz1Wa2rJW3ExGRMki0/IC7LweWN3vs8ia3TyxyXSIi0gG6QlVEJIMU7iIiGaRwFxHJIIW7iEgGKdxFRDJI4S4ikkEKdxGRDFK4i4hkkMJdRCSDFO4iIhmkcBcRySCFu4hIBincRUQySOEuIpJBCncRkQxSuIuIZJDCXUQkgxTuIiIZpHAXEckghbuISAYp3EVEMkjhLiKSQQp3EZEMUriLiGSQwl1EJIMU7iIiGaRwFxHJIIW7iEgGKdxFRDJI4S4ikkGJwt3Mcmb2ipmtN7OLW3h+LzO7I//802ZWU+xCRUQkuYLhbmaVwLXAZGAkMMPMRjY7bA7wvrsfBFwN/LTYhYqISHJJWu5jgfXuvsHdtwNLgWnNjpkG3JS/fRcw0cyseGWKiEhbVCU4ZhDwRpP7DcBRrR3j7jvMbAuwH/BuMYps6rHXHwNgwpIJxX5rEYlg9TurGTVgVOwyMidJy72lFri34xjMbJ6Z1ZtZfWNjY5L6RCTjRg0YxZmHnhm7jMxJ0nJvAIY0uT8YeKuVYxrMrAroC7zX/I3cfSGwEKC2tna38E/C/6ldLxMR6VKStNyfBYab2TAz6w5MB5Y1O2YZcE7+9mnAI+6uFBYRiaRgyz3fh34hsAKoBBa7+1ozuxKod/dlwI3ALWa2ntBin17KokVEZM+SdMvg7suB5c0eu7zJ7Y+B04tbmoiItJeuUBURySCFu4hIBincRUQySOEuIpJBCncRkQyyWNPRzawReL2dL+9PCZY2KDOdQzroHNJB55DcAe5eXeigaOHeEWZW7+61sevoCJ1DOugc0kHnUHzqlhERySCFu4hIBnXWcF8Yu4Ai0Dmkg84hHXQORdYp+9xFRGTPOmvLXURE9iDV4Z6FjbkTnMNFZrbOzF4ws4fN7IAYde5JoXNoctxpZuZmlpoZA59Jcg5mdkb+Z7HWzG4rd42FJPhdGmpmj5rZ8/nfpykx6myNmS02s01m9lIrz5uZXZM/vxfM7Ihy11hIgnM4K1/7C2b2pJkdXu4a/8zdU/lFWF74/wFfBroDa4CRzY65ALg+f3s6cEfsuttxDicAPfO353fGc8gf1xt4HFgF1Mauux0/h+HA88C++ft/EbvudpzDQmB+/vZI4LXYdTer73jgCOClVp6fAtxP2NltHPB07JrbcQ7HNPkdmhzzHNLccs/CxtwFz8HdH3X3bfm7qwg7XaVJkp8DwI+Bq4CPy1lcQknOYS5wrbu/D+Dum8pcYyFJzsGBPvnbfdl9x7So3P1xWtihrYlpwM0erAL2MbOB5akumULn4O5PfvY7ROS/5zSHe0sbcw9q7Rh33wF8tjF3WiQ5h6bmEFouaVLwHMxsNDDE3e8rZ2FtkOTnMAIYYWZPmNkqM8uVrbpkkpzDFcBMM2sg7L/wvfKUVjRt/XtJu6h/z4k264ikaBtzR5S4PjObCdQC40taUdvt8RzMrAK4GphdroLaIcnPoYrQNTOB0Nr6rZkd4u4flLi2pJKcwwxgibv/bzM7mrA72iHuvqv05RVF2v+eEzOzEwjhflysGtLccm/LxtzsaWPuiJKcA2Z2InApMNXdPylTbUkVOofewCHASjN7jdBXuixlg6pJf5fucfdP3f1/gFcIYZ8WSc5hDnAngLs/BfQgrHfSWST6e0k7MzsMWARMc/fNsepIc7hnYWPugueQ79JYQAj2tPXzQoFzcPct7t7f3WvcvYbQzzjV3evjlNuiJL9LvyYMbmNm/QndNBvKWuWeJTmHPwATAczsq4RwbyxrlR2zDJiVnzUzDtji7m/HLqotzGwocDdwtru/GrWY2KPPBUampwCvEmYJXJp/7EpCeED45f0lsB54Bvhy7JrbcQ7/F9gIrM5/LYtdc1vPodmxK0nZbJmEPwcD/g+wDngRmB675nacw0jgCcJMmtXAybFrblb/7cDbwKeEVvocoA6oa/IzuDZ/fi+m9Peo0DksAt5v8vdcH6tWXaEqIpJBae6WERGRdlK4i4hkkMJdRCSDFO4iIhmkcBcRySCFu4hIBincRUQySOEuIpJB/x9i7fXqMPdt/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "X = np.array([[0,0,1],\n",
    "              [1,0,1],\n",
    "              [1,1,1],\n",
    "              [0,1,1],\n",
    "              [0,0,1]])\n",
    "\n",
    "H = np.random.rand(3,3)\n",
    "#H/= H[2,2]\n",
    "\n",
    "Xprime = (H @ X.T).T\n",
    "Xprime/=Xprime[:,2][:,np.newaxis]\n",
    "\n",
    "plt.plot(X[:,0],X[:,1],'g-')\n",
    "plt.plot(Xprime[:,0],Xprime[:,1],'b-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do we solve for this thing given our point correspondences?  We need an equation for the eight unknowns in the matrix $H$.  By doing some algebra, we find that we can generate a linear system of equations for the coefficients of $H$:\n",
    "$$\n",
    "A \\vec{h} = 0,\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\vec{h} = \\begin{bmatrix} h_{11}, h_{12}, \\ldots, h_{32}, h_{33} \\end{bmatrix}\n",
    "$$\n",
    "and \n",
    "$$\n",
    "A = \\begin{bmatrix}  0 & 0 & 0 & -u_1 & -v_1 & -1 & v_1' u_1 & v'_1 v_1& v'_1 \\\\\n",
    "                     u_1 & v_1 & 1 &  0   & 0    &  0 & -u'_1 u_1 & -u'_1 v_1 & -u'_1 \\\\\n",
    "                     \\vdots & \\ddots & & & & & & & \\vdots \\\\\n",
    "                     0 & 0 & 0 & -u_n & -v_n & -1 & v_n' u_n & v'_n v_n& v'_n \\\\\n",
    "                     u_n & v_n & 1 &  0   & 0    &  0 & -u'_n u_n & -u'_n v_n & -u'_n \\end{bmatrix}\n",
    "$$                     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that for this matrix problem to be well-posed, $n$ (the number of correspondences) must equal 4.  Even still, this system is singular: this is because $h$ is only really defined up to a constant.  Thus, the entries of $h$ are in the nullspace of $A$.  Fortunately, finding this nullspace is easy: take the singular value decomposition of $A$, i.e.\n",
    "$$\n",
    "SVD(A) = U \\Sigma V^T.\n",
    "$$\n",
    "If the singular values (the entries in $\\Sigma$) are in descending order (as they usually are), then the nullspace of $A$ is given by the corresponding vector in $V$, or the bottom row in $V^T$.  \n",
    "\n",
    "As it turns out this SVD-based method for finding the nullspace is also applicable when $n>4$, e.g. we have more correspondences than the requisite 4.  In the case of this overdetermined system, there won'd be a nullspace, but the optimal solution for the homography matrix is still given by the column in $V$ that corresponds to the smallest singular valuem, and there is no change to the routine (except that $A$ has more rows). \n",
    "\n",
    "**Your assignment is to generate code that takes a set of point correspondences $\\{(u_i,v_i, u'_i,v'_i)\\}$ for $i$ in $n$, and solves for the homography matrix between them.  You can directly apply this to the matches that you've found via your stitching code so far, but you will find that this procedure is extremely sensitive to outliers.  To determine if your code is working properly, use the synthetic correspondences generated in the homography demo above to recover the matrix $H$ from $X$ and $X'$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_random = np.random.rand(8,9)\n",
    "U,Sigma,Vt = np.linalg.svd(A_random)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}