{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Solving 2D homography from 4 points correspondences"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Small experiment to solve a 2D homography turning any quadrilateral into another quadrilateral"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import numpy as np\n",
"import numpy.linalg as la\n",
"import pylab as pl\n",
"%matplotlib inline\n",
"np.set_printoptions(suppress=True, precision=5)"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Source and target quadrilaterals\n",
"x = np.array([\n",
" [2, 2],\n",
" [2, 4],\n",
" [4, 4],\n",
" [4, 2]\n",
"])\n",
"\n",
"xp = np.array([\n",
" [2, 2],\n",
" [2.5, 4],\n",
" [3.5, 3.5],\n",
" [4, 2]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.PathCollection at 0x10f96ef60>"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAl0AAAEyCAYAAADAyGU5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGShJREFUeJzt3W+MnGd57/HvL7ZTcBI1UrMHrDiOOVJeNCAIYeRGJAIS\nKcih0KgSL4zcIKFWVjggQQ9qDyUSR7zIKyTEgUMb7QFUEKZRdJK0UZRAgxqppCgm69QkOAnIJ80/\nK1I2UPKnRiCH67yYJ3RYdr2z3vE9z3i/H2k0M/dzzex97zPPpZ9n5xmnqpAkSdKpdca0JyBJkrQR\nGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDWye9gSWc955\n59XOnTunPQ1JjRw8ePD5qpqb9jwmwf4lbTzj9rBehq6dO3eysLAw7WlIaiTJk9Oew6TYv6SNZ9we\n5p8XJUmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYMXZIkSQ2MHbqSbEry\nr0nuXGZbknwhyZEkDyW5dGTb7iQ/6rZ9clITl6Rx2b8k9cFa3un6GPDoCtuuAS7qLvuAv4FhowO+\n1G2/GPhAkotPerZL7N8PO3fCGWcMr/fvn9QzS1qrnh+PvetfG0bPXxhSS2OFriTbgT8EvrxCybXA\n12vofuDcJNuAXcCRqnq8qn4J3NzVrtv+/bBvHzz5JFQNr/ft83iWpqHPx2Mf+9eG0ecXhjQF477T\n9XngL4FfrbD9fODpkfvPdGMrja/bDTfAsWO/OXbs2HBcUls9Px571782jJ6/MKTWVg1dSd4LPFdV\nB0/lRJLsS7KQZGFxcXHV+qeeWtu4pFOnr8djX/vXhtHXF4Y0JeO803U58EdJnmD49vpVSb6xpOYo\ncMHI/e3d2Erjv6Wq5qtqUFWDubm5VSe1Y8faxiWdOj0+HnvZvzaMHr8wpGlYNXRV1V9V1faq2gns\nAf6pqv5kSdkdwAe7s4AuA16oqmeBB4CLkrwhyZnd4++YxMRvvBG2bv3Nsa1bh+OS2urr8djX/rVh\n9PWFIU3JSX9PV5Lrk1zf3b0LeBw4Avwf4L8BVNVx4KPAtxmeOXRLVR1e14w7e/fC/DxceCEkw+v5\n+eG4pLZm7Xicdv/aMGbthSGdYqmqac/htwwGg1pYWJj2NCQ1kuRgVQ2mPY9JsH9JG8+4PcxvpJck\nSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLU\ngKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFD\nlyRJUgOGLkmSpAYMXZIkSQ1sXq0gyWuAfwZ+p6v/v1X1P5fU/AWwd+Q5fx+Yq6qfJnkCeAl4BThe\nVYPJTV+SVmb/ktQnq4Yu4BfAVVX1cpItwH1J7q6q+18tqKrPAp8FSPI+4M+r6qcjz3FlVT0/yYlL\n0hjsX5J6Y9XQVVUFvNzd3dJd6gQP+QDwd+ufmiStj/1LUp+M9ZmuJJuSHAKeA+6pqgMr1G0FdgO3\njgwX8J0kB5PsO8HP2JdkIcnC4uLi+CuQpBOwf0nqi7FCV1W9UlWXANuBXUnetELp+4B/WfLW/BXd\nY68BPpLkHSv8jPmqGlTVYG5ubg1LkKSV2b8k9cWazl6sqp8B9zL81+By9rDkrfmqOtpdPwfcDuxa\n+zQlaX3sX5KmbdXQlWQuybnd7dcCVwOPLVP3u8A7gX8YGTsryTmv3gbeDfxwMlOXpBOzf0nqk3HO\nXtwGfC3JJoYh7ZaqujPJ9QBVdVNX98fAP1bVf4w89nXA7Ule/VnfrKpvTWz2knRi9i9JvZHhyT39\nMhgMamFhYdrTkNRIkoOny3dg2b+kjWfcHuY30kuSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVID\nhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxd\nkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAZWDV1JXpPk+0l+\nkORwks8sU/OuJC8kOdRdPj2ybXeSHyU5kuSTk16AJK3E/iVNyP79sHMnnHHG8Hr//mnPaCZtHqPm\nF8BVVfVyki3AfUnurqr7l9R9t6reOzqQZBPwJeBq4BnggSR3VNUjk5i8JK3C/iWt1/79sG8fHDs2\nvP/kk8P7AHv3Tm9eM2jVd7pq6OXu7pbuUmM+/y7gSFU9XlW/BG4Grj2pmUrSGtm/pAm44Yb/DFyv\nOnZsOK41GeszXUk2JTkEPAfcU1UHlil7e5KHktyd5I3d2PnA0yM1z3Rjy/2MfUkWkiwsLi6uYQmS\ntDL7l7ROTz21tnGtaKzQVVWvVNUlwHZgV5I3LSl5ENhRVW8Gvgj8/VonUlXzVTWoqsHc3NxaHy5J\ny7J/Seu0Y8faxrWiNZ29WFU/A+4Fdi8Zf/HVt/Cr6i5gS5LzgKPABSOl27sxSWrK/iWdpBtvhK1b\nf3Ns69bhuNZknLMX55Kc291+LcMPlT62pOb1SdLd3tU970+AB4CLkrwhyZnAHuCOyS5BkpZn/5Im\nYO9emJ+HCy+EZHg9P++H6E/COGcvbgO+1p3JcwZwS1XdmeR6gKq6CXg/8OEkx4GfA3uqqoDjST4K\nfBvYBHy1qg6fioVI0jLsX9Ik7N1ryJqADHtLvwwGg1pYWJj2NCQ1kuRgVQ2mPY9JsH9JG8+4Pcxv\npJckSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAl\nSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5Ik\nqQFDlyRJUgOGLkmSpAZWDV1JXpPk+0l+kORwks8sU7M3yUNJHk7yvSRvGdn2RDd+KMnCpBcgSSux\nf0nqk81j1PwCuKqqXk6yBbgvyd1Vdf9Izb8B76yqf09yDTAP/MHI9iur6vnJTVuSxmL/ktQbq4au\nqirg5e7ulu5SS2q+N3L3fmD7pCYoSSfL/iWpT8b6TFeSTUkOAc8B91TVgROU/ylw98j9Ar6T5GCS\nfSf4GfuSLCRZWFxcHGdakrQq+5ekvhgrdFXVK1V1CcN/Ae5K8qbl6pJcybBp/Y+R4Su6x14DfCTJ\nO1b4GfNVNaiqwdzc3JoWIUkrsX9J6os1nb1YVT8D7gV2L92W5M3Al4Frq+onI4852l0/B9wO7FrP\nhCXpZNi/JE3bOGcvziU5t7v9WuBq4LElNTuA24DrqurHI+NnJTnn1dvAu4EfTm76krQy+5ekPhnn\n7MVtwNeSbGIY0m6pqjuTXA9QVTcBnwZ+D/jrJADHq2oAvA64vRvbDHyzqr41+WVI0rLsX5J6I8OT\ne/plMBjUwoJfiSNtFEkOdkFn5tm/pI1n3B7mN9JLkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElS\nA4YuSZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYM\nXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGVg1dSV6T5PtJ\nfpDkcJLPLFOTJF9IciTJQ0kuHdm2O8mPum2fnPQCJGkl9i9JfTLOO12/AK6qqrcAlwC7k1y2pOYa\n4KLusg/4G4Akm4AvddsvBj6Q5OIJzV2SVmP/ktQbq4auGnq5u7ulu9SSsmuBr3e19wPnJtkG7AKO\nVNXjVfVL4OauVpJOOfuXpD4Z6zNdSTYlOQQ8B9xTVQeWlJwPPD1y/5lubKXx5X7GviQLSRYWFxfH\nnb8knZD9S1JfjBW6quqVqroE2A7sSvKmSU+kquaralBVg7m5uUk/vaQNyv4lqS/WdPZiVf0MuBfY\nvWTTUeCCkfvbu7GVxiWpKfuXpGkb5+zFuSTndrdfC1wNPLak7A7gg91ZQJcBL1TVs8ADwEVJ3pDk\nTGBPVytJp5z9S1KfbB6jZhvwte5MnjOAW6rqziTXA1TVTcBdwHuAI8Ax4EPdtuNJPgp8G9gEfLWq\nDk9+GZK0LPuXpN5I1dITeaZvMBjUwsLCtKchqZEkB6tqMO15TIL9S9p4xu1hfiO9JElSA4YuSZKk\nBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAYMXZIkSQ0Y\nuiRJkhowdEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkNGLokSZIaMHRJ\nkiQ1sHm1giQXAF8HXgcUMF9V/2tJzV8Ae0ee8/eBuar6aZIngJeAV4DjVTWY3PQlaWX2L0l9smro\nAo4Dn6iqB5OcAxxMck9VPfJqQVV9FvgsQJL3AX9eVT8deY4rq+r5SU5cksZg/5LUG6v+ebGqnq2q\nB7vbLwGPAuef4CEfAP5uMtOTpJNn/5LUJ2v6TFeSncBbgQMrbN8K7AZuHRku4DtJDibZd4Ln3pdk\nIcnC4uLiWqYlSauyf0matrFDV5KzGTajj1fViyuUvQ/4lyVvzV9RVZcA1wAfSfKO5R5YVfNVNaiq\nwdzc3LjTkqRV2b8k9cFYoSvJFoYNa39V3XaC0j0seWu+qo52188BtwO7Tm6qkrR29i9JfbFq6EoS\n4CvAo1X1uRPU/S7wTuAfRsbO6j68SpKzgHcDP1zvpCVpHPYvSX0yztmLlwPXAQ8nOdSNfQrYAVBV\nN3Vjfwz8Y1X9x8hjXwfcPux7bAa+WVXfmsTEJWkM9i9JvbFq6Kqq+4CMUfe3wN8uGXsceMtJzk2S\n1sX+JalP/EZ6SZKkBgxdkiRJDRi6JEmSGjB0SZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOG\nLkmSpAYMXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2S\nJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUwKqhK8kFSe5N8kiSw0k+tkzNu5K8kORQd/n0yLbd\nSX6U5EiST056AZK0EvuXpD7ZPEbNceATVfVgknOAg0nuqapHltR9t6reOzqQZBPwJeBq4BnggSR3\nLPNYSToV7F+SemPVd7qq6tmqerC7/RLwKHD+mM+/CzhSVY9X1S+Bm4FrT3aykrQW9i9JfbKmz3Ql\n2Qm8FTiwzOa3J3koyd1J3tiNnQ88PVLzDCs0vCT7kiwkWVhcXFzLtCRpVfYvSdM2duhKcjZwK/Dx\nqnpxyeYHgR1V9Wbgi8Dfr3UiVTVfVYOqGszNza314ZK0IvuXpD4YK3Ql2cKwYe2vqtuWbq+qF6vq\n5e72XcCWJOcBR4ELRkq3d2OS1IT9S1JfjHP2YoCvAI9W1edWqHl9V0eSXd3z/gR4ALgoyRuSnAns\nAe6Y1OQl6UTsX5L6ZJyzFy8HrgMeTnKoG/sUsAOgqm4C3g98OMlx4OfAnqoq4HiSjwLfBjYBX62q\nwxNegyStxP4lqTcy7C39MhgMamFhYdrTkNRIkoNVNZj2PCbB/iVtPOP2ML+RXpIkqQFDlyRJUgOG\nLkmSpAYMXZIkSQ0YuiRJkhowdEmSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2S\nJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmS\nGlg1dCW5IMm9SR5JcjjJx5ap2ZvkoSQPJ/lekreMbHuiGz+UZGHSC5Ckldi/JPXJ5jFqjgOfqKoH\nk5wDHExyT1U9MlLzb8A7q+rfk1wDzAN/MLL9yqp6fnLTlqSx2L8k9caqoauqngWe7W6/lORR4Hzg\nkZGa74085H5g+4TnKUlrZv+S1Cdr+kxXkp3AW4EDJyj7U+DukfsFfCfJwST71jpBSZoE+5ekaRvn\nz4sAJDkbuBX4eFW9uELNlQyb1hUjw1dU1dEk/wW4J8ljVfXPyzx2H7APYMeOHWtYgiSdmP1LUh+M\n9U5Xki0MG9b+qrpthZo3A18Grq2qn7w6XlVHu+vngNuBXcs9vqrmq2pQVYO5ubm1rUKSVmD/ktQX\n45y9GOArwKNV9bkVanYAtwHXVdWPR8bP6j68SpKzgHcDP5zExCVpNfYvSX0yzp8XLweuAx5Ocqgb\n+xSwA6CqbgI+Dfwe8NfDHsfxqhoArwNu78Y2A9+sqm9NdAWStDL7l6TeGOfsxfuArFLzZ8CfLTP+\nOPCW336EJJ169i9JfeI30kuSJDVg6JIkSWrA0CVJktSAoUuSJKkBQ5ckSVIDhi5JkqQGDF2SJEkN\nGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA4YuSZKkBgxdkiRJDRi6JEmSGjB0\nSZIkNWDokiRJasDQJUmS1IChS5IkqQFDlyRJUgOGLkmSpAZWDV1JLkhyb5JHkhxO8rFlapLkC0mO\nJHkoyaUj23Yn+VG37ZOTnPz+/bBzJ5xxxvB6//5JPruktejj8djn/gX085cmbVQtjseqOuEF2AZc\n2t0+B/gxcPGSmvcAdwMBLgMOdOObgP8H/FfgTOAHSx+73OVtb3tbreYb36jaurUK/vOydetwXFJb\n6z0egYVapS+czKWv/WsivzRJk7PO43HcHrbqO11V9WxVPdjdfgl4FDh/Sdm1wNe7n30/cG6SbcAu\n4EhVPV5VvwRu7mrX7YYb4Nix3xw7dmw4Lqmtvh6Pfe1fQH9/adJG1Oh4XNNnupLsBN4KHFiy6Xzg\n6ZH7z3RjK40v99z7kiwkWVhcXFx1Lk89tbZxSafOLByPfepfwGz80qSNotHxOHboSnI2cCvw8ap6\ncaKzAKpqvqoGVTWYm5tbtX7HjrWNSzp1+n489q1/Af3/pUkbSaPjcazQlWQLw4a1v6puW6bkKHDB\nyP3t3dhK4+t2442wdetvjm3dOhyX1Fafj8c+9i+g3780aaNpdDyOc/ZigK8Aj1bV51YouwP4YHcW\n0GXAC1X1LPAAcFGSNyQ5E9jT1a7b3r0wPw8XXgjJ8Hp+fjguqa2+Ho997V9Af39p0kbU6HjM8EP3\nJyhIrgC+CzwM/Kob/hSwA6Cqbuoa2/8GdgPHgA9V1UL3+PcAn2d4JtBXq2rV2DgYDGphYeGkFiRp\n9iQ5WFWDU/C89i9Jp9y4PWzzagVVdR/DU6lPVFPAR1bYdhdw12o/R5Imzf4lqU/8RnpJkqQGDF2S\nJEkNGLokSZIaMHRJkiQ1YOiSJElqwNAlSZLUgKFLkiSpAUOXJElSA6t+I/00JFkEnlzDQ84Dnj9F\n0+kL13h6cI3Lu7CqxvyfovvtJPoX+Lo4XbjG08Mp62G9DF1rlWThVPwXIn3iGk8PrlHL2Qi/M9d4\nenCN6+OfFyVJkhowdEmSJDVwuoSu+WlPoAHXeHpwjVrORviducbTg2tch9PiM12SJEl9d7q80yVJ\nktRrhi5JkqQGZiZ0Jbkgyb1JHklyOMnHlqlJki8kOZLkoSSXTmOuJ2vMNb4ryQtJDnWXT09jricr\nyWuSfD/JD7o1fmaZmlnfj+Oscab346uSbEryr0nuXGbbTO/HSbOH/bpmZl/79q9f18zsPhw1lf5V\nVTNxAbYBl3a3zwF+DFy8pOY9wN1AgMuAA9Oe9ylY47uAO6c913WsMcDZ3e0twAHgstNsP46zxpne\njyPr+O/AN5dby6zvx1Pwu7KH1Wy/9u1fs78Pl6yjef+amXe6qurZqnqwu/0S8Chw/pKya4Gv19D9\nwLlJtjWe6kkbc40zrds3L3d3t3SXpWdzzPp+HGeNMy/JduAPgS+vUDLT+3HS7GGzz/51+phW/5qZ\n0DUqyU7grQwT+KjzgadH7j/DjB7wJ1gjwNu7tzvvTvLGphObgO4t3UPAc8A9VXXa7ccx1ggzvh+B\nzwN/Cfxqhe0zvx9PFXvY7L727V+/NrP7sDOV/jVzoSvJ2cCtwMer6sVpz+dUWGWNDwI7qurNwBeB\nv289v/Wqqleq6hJgO7AryZumPadJG2ONM70fk7wXeK6qDk57LrPGHjbbr337FzDj+3Ca/WumQleS\nLQwP5P1VddsyJUeBC0bub+/GZsZqa6yqF19967eq7gK2JDmv8TQnoqp+BtwL7F6yaeb346tWWuNp\nsB8vB/4oyRPAzcBVSb6xpOa02Y+TYg87LV77gP1rxvfh1PrXzISuJAG+AjxaVZ9boewO4IPdWQeX\nAS9U1bPNJrlO46wxyeu7OpLsYrgPf9JuluuTZC7Jud3t1wJXA48tKZv1/bjqGmd9P1bVX1XV9qra\nCewB/qmq/mRJ2Uzvx0mzh/26ZmZf+/avX9fM7D6E6favzet9goYuB64DHu7+1gzwKWAHQFXdBNzF\n8IyDI8Ax4ENTmOd6jLPG9wMfTnIc+Dmwp6pm6UOO24CvJdnE8EC9paruTHI9nDb7cZw1zvp+XNZp\nth8nzR42+699+9fs78MVtdiP/jdAkiRJDczMnxclSZJmmaFLkiSpAUOXJElSA4YuSZKkBgxdkiRJ\nDRi6JEmSGjB0SZIkNfD/AXwJvqU6nc4WAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10f632358>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pl.figure(figsize=(10, 5))\n",
"pl.subplot(121)\n",
"pl.scatter(x[:,0], x[:,1], c='b')\n",
"pl.subplot(122)\n",
"pl.scatter(xp[:,0], xp[:,1], c='r')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solving for the homography\n",
"\n",
"See _Hartley & Zissermann \"Multiple view geometry in computer vision\"_, section 4.1 and equation 4.3 for reference.\n",
"\n",
"We want to solve for the homography $H$ such that\n",
"$$x_i' = Hx_i$$\n",
"Which can be rewritten as a cross product\n",
"$$x_i' \\times Hx_i = 0$$\n",
"\n",
"With the $j-th$ row of H denoted as $h^{jT}$, then we have\n",
"\n",
"$$Hx_i = \\left( \\begin{array}{c} h^{1T}x_i \\\\ h^{2T}x_i \\\\ h^{3T}x_i\\end{array} \\right)$$\n",
"\n",
"So now the cross-product above can be rewritten (assuming $x_i'(k)$ denotes the k-th coordinate of $x_i'$)\n",
"\n",
"$$\n",
"x_i' \\times Hx_i = \\left(\n",
" \\begin{array}{c}\n",
" x_i'(1) h^{3T}x_i - x_i'(2) h^{2T}x_i \\\\\n",
" x_i'(2) h^{1T}x_i - x_i'(0) h^{3T}x_i \\\\\n",
" x_i'(0) h^{2T}x_i - x_i'(1) h^{1T}x_i\n",
" \\end{array} \\right)\n",
"$$\n",
"\n",
"Since $h^{jT}x_i = x_i^T h^j$, we can rewrite as\n",
"\n",
"$$\n",
"\\left(\n",
"\\begin{array}{c c c}\n",
"0^T & -x_i'(2)x_i^T & x_i'(1)x_i^T \\\\\n",
"x_i'(2)x_i^T & 0^T & -x_i'(0)x_i^T \\\\\n",
"-x_i'(1)x_i^T & x_i'(0)x_i^T & 0^T\n",
"\\end{array}\n",
"\\right)\n",
"\\\n",
"\\left(\n",
"\\begin{array}{c}\n",
"h^1 \\\\\n",
"h^2 \\\\\n",
"h^3\n",
"\\end{array}\n",
"\\right)\n",
"\\\n",
"=0\n",
"$$\n",
"\n",
"Note that this is a 3x9 matrix multiplying a 9x1 vector containing the unknown (the row of H).\n",
"\n",
"The third row can be obtaine as a combination of first and second and therefore we can drop it, so we have two remaining equations per point\n",
"$$\n",
"\\left(\n",
"\\begin{array}{c c c}\n",
"0^T & -x_i'(2)x_i^T & x_i'(1)x_i^T \\\\\n",
"x_i'(2)x_i^T & 0^T & -x_i'(0)x_i^T\n",
"\\end{array}\n",
"\\right)\n",
"\\\n",
"\\left(\n",
"\\begin{array}{c}\n",
"h^1 \\\\\n",
"h^2 \\\\\n",
"h^3\n",
"\\end{array}\n",
"\\right)\n",
"\\\n",
"=0\n",
"$$\n",
"\n",
"Which we'll solve as an $Ah = 0$ problem"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-0.22664 -0.17434 0.62761]\n",
" [-0.06973 -0.33124 0.62761]\n",
" [-0.03487 -0.06102 0.1046 ]]\n"
]
}
],
"source": [
"def compute_A_for_point(x, xp):\n",
" \"\"\"\n",
" Computes two rows of A for a single point correspondence x and xp\n",
" \"\"\"\n",
" assert len(x) == len(xp) == 2\n",
" # convert to homogeneous\n",
" x = np.r_[x, 1]\n",
" xp = np.r_[xp, 1]\n",
" \n",
" return np.array([\n",
" np.r_[np.zeros(3), -xp[2] * x , xp[1] * x],\n",
" np.r_[xp[2] * x , np.zeros(3), -xp[0] * x]\n",
" ])\n",
"\n",
"def compute_homography(x, xp):\n",
" \"\"\"\n",
" Given two arrays of corresponding points as 4x2 arrays, compute the\n",
" homography H such that xp = Hx\n",
" \"\"\"\n",
" A = np.concatenate([compute_A_for_point(*pair) for pair in zip(x, xp)], axis=0)\n",
" U, s, V = la.svd(A)\n",
" # Solution is the last row of V (corresponding to the smallest eigenvalue)\n",
" H = np.reshape(V[-1], (3, 3))\n",
" return H\n",
" \n",
"compute_A_for_point(x[0], xp[0])\n",
"H = compute_homography(x, xp)\n",
"\n",
"print(H)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-- 0\n",
"x = [2 2]\n",
"xp = [ 2. 2.]\n",
"Hx = [ 2. 2.]\n",
"-- 1\n",
"x = [2 4]\n",
"xp = [ 2.5 4. ]\n",
"Hx = [ 2.5 4. ]\n",
"-- 2\n",
"x = [4 4]\n",
"xp = [ 3.5 3.5]\n",
"Hx = [ 3.5 3.5]\n",
"-- 3\n",
"x = [4 2]\n",
"xp = [ 4. 2.]\n",
"Hx = [ 4. 2.]\n"
]
}
],
"source": [
"def hdot(M, p):\n",
" \"\"\"Convenience function to apply a 3x3 transform to a 2D point\"\"\"\n",
" pp = M.dot(np.r_[p, 1]) # to homogeneous\n",
" return pp[:-1] / pp[-1] # back to inhomogeneous\n",
"\n",
"for i in range(4):\n",
" print(\"-- %d\" % i)\n",
" print(\"x = %s\" % x[i])\n",
" print(\"xp = %s\" % xp[i])\n",
" print(\"Hx = %s\" % (hdot(H, x[i])))\n",
" assert np.allclose(hdot(H, x[i]), xp[i])"
]
}
],
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}