{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Affine transformation\n",
    "**Students:** Sam van Kampen, Robin Wacanno\n",
    "\n",
    "**Student numbers:** 11874716, 11741163"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The affine transformation matrix $A$ is of the form\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix}a & b & c\\\\d & e & f\\\\0 & 0 & 1\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "and maps input points $x_i, y_i$ to output points $x_{i}', y_{i}'$ with\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}x_i'\\\\y_i'\\\\1\\end{pmatrix} =\n",
    "\\begin{pmatrix}a & b & c\\\\d & e & f\\\\0 & 0 & 1\\end{pmatrix}\n",
    "\\begin{pmatrix}x_i\\\\y_i\\\\1\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "working this out and rewriting so that the unknown vector is on the right hand side, which makes the system easily solvable as a system of linear equations (given enough points):\n",
    "\n",
    "$$\n",
    "\\begin{split}\\begin{pmatrix}x_i'\\\\ y_i'\\end{pmatrix} =\n",
    "\\begin{pmatrix}x_i& y_i& 1& 0& 0& 0\\\\0& 0& 0& x_i& y_i& 1 \\end{pmatrix}\n",
    "\\begin{pmatrix}a\\\\b\\\\c\\\\d\\\\e\\\\f \\end{pmatrix}\\end{split}\n",
    "$$\n",
    "\n",
    "This system does not have enough point correspondences to be invertible, so we need to go from one point to three (non-collinear) points, yielding:\n",
    "\n",
    "\\begin{split}\\begin{pmatrix}x_1'\\\\ y_1'\\\\x_2'\\\\ y_2'\\\\x_3'\\\\ y_3'\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "x_1& y_1& 1& 0& 0& 0\\\\0& 0& 0& x_1& y_1& 1\\\\\n",
    "x_2& y_2& 1& 0& 0& 0\\\\0& 0& 0& x_2& y_2& 1\\\\\n",
    "x_3& y_3& 1& 0& 0& 0\\\\0& 0& 0& x_3& y_3& 1 \\end{pmatrix}\n",
    "\\begin{pmatrix}a\\\\b\\\\c\\\\d\\\\e\\\\f \\end{pmatrix}\\end{split}\n",
    "\n",
    "This system can be solved easily using numpy etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def affineTransform(f, x1, y1, x2, y2, x3, y3, width, height, crop=True):\n",
    "    M = np.array([[x, y, 1, *[0]*6, x, y, 1] for x, y in ((x1, y1), (x2, y2), (x3, y3))]).reshape(6, 6)\n",
    "    new_points = np.array([0, 0, width - 1, 0, 0, height - 1])\n",
    "    \n",
    "    M_inverse = np.linalg.inv(M)\n",
    "    \n",
    "    affine_transform_parameters = M_inverse @ new_points\n",
    "    \n",
    "    transform_matrix = np.array([*affine_transform_parameters, 0, 0, 1]).reshape(3, 3)\n",
    "    \n",
    "    # skimage.transform.warp does a backward transform, so we need the inverse transformation matrix.\n",
    "    inv_transform_matrix = np.linalg.inv(transform_matrix)\n",
    "    \n",
    "    # Warp, then crop to (height, width)\n",
    "    warped_img = skimage.transform.warp(f, inv_transform_matrix)\n",
    "    \n",
    "    return warped_img[:height,:width] if crop else warped_img   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test the transform code, we load some test images and warp them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0830bb8780>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 7664 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import skimage\n",
    "\n",
    "f = skimage.io.imread('data/sheared.png')\n",
    "# Corner points corresponding to new corners:\n",
    "f_points = [6, 87, 490, 87, 150, 401]\n",
    "\n",
    "imshow(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0830b58630>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 6188 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(affineTransform(f, *f_points, 486, 314))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to show I'm not cheating and cropping to a rectangular area containing all blue instead of properly transforming:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0830ab89e8>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 7500 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(affineTransform(f, *f_points, 486, 314, crop=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also do a rotation around the center, for good measure. Instead of using my own image I'll just rotate Lena."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0830a94908>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 155924 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# They removed Lena from the default test images!\n",
    "# (I understand the reasons, but I just think it's a nice tradition to use her image)\n",
    "f2 = skimage.io.imread('data/lena.png')\n",
    "\n",
    "imshow(f2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's rotate Lena 45 degrees, clockwise. To do this, the measurements we give have to be the inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "f2_rotated = affineTransform(f2, -106, 256, 256, -106, 256, 618, 512, 512)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f08301fd278>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 128624 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "imshow(f2_rotated)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It seems to work fine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perspective Transform\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": 10,
   "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": 11,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f08301e27f0>"
      ]
     },
     "execution_count": 12,
     "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": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f0830143358>"
      ]
     },
     "execution_count": 13,
     "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f083011e710>"
      ]
     },
     "execution_count": 14,
     "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
