{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perspective Transform\n",
    "Students: Sam van Kampen, Robin Wacanno\n",
    "\n",
    "Studentnrs: 11874716, 11741163\n",
    "\n",
    "$\\newcommand{\\T}[0]{^\\intercal}$\n",
    "\n",
    "Voor een generieke transformatie worden punten $(x_i, y_i)$ naar punten $(x_i', y_i')$ omgezet met de volgende matrixoperatie, waarbij $s$ een scalar is :\n",
    "\\begin{split}s\\begin{pmatrix}x'\\\\y'\\\\1\\end{pmatrix} =&\n",
    "\\begin{pmatrix}a & b & c\\\\d & e & f\\\\g & h & i\\end{pmatrix}\n",
    "\\begin{pmatrix}x\\\\y\\\\1\\end{pmatrix} \\\\\n",
    "\\widetilde x' \\sim& P\\,\\widetilde x\\end{split}\n",
    "\n",
    "Deze operatie kan herschreven worden als:\n",
    "\\begin{split}x' =& \\frac{ax+by+c}{gx+hy+i}\\\\\n",
    "y' =& \\frac{dx+ey+f}{gx+hy+i}\\end{split}\n",
    "\n",
    "Door de bovenstaande vergelijkingen te herschrijven zodat aan een kant enkel een 0 komt te staan kan deze herschreven worden tot de volgende matrixoperatie:\n",
    "\\begin{split}\\begin{pmatrix}\n",
    "x & y & 1 & 0 & 0 & 0 & -x'x & -x'y & -x'\\\\\n",
    "0 & 0 & 0 & x & y & 1 & -y'x & -y'y & -y'\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}a\\\\b\\\\c\\\\d\\\\e\\\\f\\\\g\\\\h\\\\i \\end{pmatrix} =\n",
    "\\begin{pmatrix}0\\\\0 \\end{pmatrix}\\end{split}\n",
    "\n",
    "Het bovenstaande is echter voor 1 punt, terwijl we voor de transformatie een matrix moeten vaststellen die 4 gegeven punten $(x_i, y_i)$ naar 4 andere gegeven punten $(x_i', y_i')$ omzet. Om dit voor 4 punten gedaan te krijgen moet de bovenstaande operatie tot het volgende herschreven worden:\n",
    "\\begin{split}\\begin{pmatrix}\n",
    "x_1 & y_1 & 1 & 0 & 0 & 0 & -x_1'x_1 & -x_1'y_1 & -x_1'\\\\\n",
    "0 & 0 & 0 & x_1 & y_1 & 1 & -y_1'x_1 & -y_1'y_1 & -y_1'\\\\\n",
    "x_2 & y_2 & 1 & 0 & 0 & 0 & -x_2'x_2 & -x_2'y_2 & -x_2'\\\\\n",
    "0 & 0 & 0 & x_2 & y_2 & 1 & -y_2'x_2 & -y_2'y_2 & -y_2'\\\\\n",
    "x_3 & y_3 & 1 & 0 & 0 & 0 & -x_3'x_3 & -x_3'y_3 & -x_3'\\\\\n",
    "0 & 0 & 0 & x_3 & y_3 & 1 & -y_3'x_3 & -y_3'y_3 & -y_3'\\\\\n",
    "x_4 & y_4 & 1 & 0 & 0 & 0 & -x_4'x_4 & -x_4'y_4 & -x_4'\\\\\n",
    "0 & 0 & 0 & x_4 & y_4 & 1 & -y_4'x_4 & -y_4'y_4 & -y_4'\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}a\\\\b\\\\c\\\\d\\\\e\\\\f\\\\g\\\\h\\\\i \\end{pmatrix} =&\n",
    "\\begin{pmatrix}0\\\\0\\\\0\\\\0\\\\0\\\\0\\\\0\\\\0 \\end{pmatrix} \\\\\n",
    "M p =& 0\\end{split}\n",
    "\n",
    "Gezien de kans op een daadwerkelijke nulvector nihil is gezien het voorkomen van ruis, moet gezorgd worden dat de uitkomst van de operatie minimaal is. We zoeken echter niet naar een triviale uitkomst, dus:\n",
    "\\min \\|M p\\| \\quad\\text{s.t.}\\quad \\|p\\| = 1\n",
    "\n",
    "Laten we nu de SVD van $M$ nemen:\n",
    "\\begin{split}\\min \\|UDV\\T p\\| \\quad\\text{s.t.}\\quad \\|p\\| = 1\\\\\n",
    "\\min \\|DV\\T p\\| \\quad\\text{s.t.}\\quad \\|p\\| = 1\\end{split}\n",
    "\n",
    "gezien $U$ normbehoudend is omdat deze orthogonaal is. Als nu $q = V^Tp$:\n",
    "\\begin{split}\\min \\|D q\\| \\quad\\text{s.t.}\\quad \\|V q\\| = 1\\\\\n",
    "\\min \\|D q\\| \\quad\\text{s.t.}\\quad \\|q\\| = 1\\\\\\end{split}\n",
    "daar ook $V$ orthogonaal is.\n",
    "\n",
    "Hieruit volgt dat $p = Vq$, dus de laatste rij uit $V\\T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage.transform import warp\n",
    "from skimage.io import imread"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perspectiveTransform(f, x1, y1, x2, y2, x3, y3, x4, y4, width, height):\n",
    "    w = width - 1\n",
    "    h = height - 1\n",
    "    M = np.array([\n",
    "        [x1, y1, 1, 0, 0, 0, 0, 0, 0],\n",
    "        [0, 0, 0, x1, y1, 1, 0, 0, 0],\n",
    "        [x2, y2, 1, 0, 0, 0, -w*x2, -w*y2, -w],\n",
    "        [0, 0, 0, x2, y2, 1, 0, 0, 0],\n",
    "        [x3, y3, 1, 0, 0, 0, 0, 0, 0],\n",
    "        [0, 0, 0, x3, y3, 1, -h*x3, -h*y3, -h],\n",
    "        [x4, y4, 1, 0, 0, 0, -w*x4, -w*y4, -w],\n",
    "        [0, 0, 0, x4, y4, 1, -h*x4, -h*y4, -h]\n",
    "    ])\n",
    "    U, D, Vt = np.linalg.svd(M)\n",
    "    trans = Vt[-1].reshape(3, 3)\n",
    "    return warp(f, np.linalg.inv(trans))[:h, :w]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd14285ed30>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 132272 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = imread(\"data/flyeronground.png\")\n",
    "plt.imshow(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd1424d5ba8>"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 68244 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f_trans = perspectiveTransform(f, 570, 188, 823, 172, 348, 555, 595, 588, 330, 434)\n",
    "plt.imshow(f_trans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fd1423e15f8>"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 70080 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f_trans = perspectiveTransform(f, 149, 590, 104, 377, 392, 394, 303, 207, 400, 250)\n",
    "plt.imshow(f_trans)"
   ]
  },
  {
   "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
