PY2-linear-transformations-2d.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Main Sources** \n",
    "    1. https://notgnoshi.github.io/linear-transformations/\n",
    "    2. https://dododas.github.io/linear-algebra-with-python/posts/16-12-29-2d-transformations.html\n",
    "    \n",
    "**Further Reading**   \n",
    "    [1] https://mmas.github.io/linear-transformations-numpy   \n",
    "    [2] [Affine Image Transformations in Python with Numpy, Pillow and OpenCV](https://stackabuse.com/affine-image-transformations-in-python-with-numpy-pillow-and-opencv/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center>2D Linear Transformations</center> #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A linear transformation of the plane $\\mathbb{R}^2$ is a geometric transformation of the form\n",
    "$$ \n",
    "  A\\vec{x} = \\vec{w}\\ \\ \\  \\Leftrightarrow \\ \\ \\ \n",
    "  \\left[ \\begin{array}{cc}\n",
    "                a & b\\\\ c & d\n",
    "           \\end{array} \n",
    "  \\right]\n",
    "   \\left[ \\begin{array}{c}\n",
    "                x\\\\ y\n",
    "           \\end{array} \n",
    "  \\right] = \n",
    "  \\left[ \\begin{array}{c}\n",
    "                u\\\\ v\n",
    "           \\end{array} \n",
    "  \\right] \n",
    "$$\n",
    "where $a, b, c$ and $d$ are real numbers.\n",
    "\n",
    "Linear transformations leave the origin fixed and preserve parallelism. As an example,\n",
    "consider the matrix\n",
    "$$\n",
    "  A = \\left[\\begin{array}{cc}\n",
    "          2 & -1 \\\\ 1 & 1\n",
    "          \\end{array}\n",
    "      \\right]\n",
    "$$\n",
    "that transforms an arbitrary vector $\\left[ \\begin{array}{c}\n",
    "                x\\\\ y\n",
    "           \\end{array} \\right]$\n",
    "to a linear combination of the column vectors of $A$:\n",
    "$$\n",
    "  \\left[ \\begin{array}{c}\n",
    "                u\\\\ v\n",
    "           \\end{array} \n",
    "  \\right] =\n",
    "  \\left[\\begin{array}{cc}\n",
    "          2 & -1 \\\\ 1 & 1\n",
    "          \\end{array}\n",
    "   \\right]       \n",
    "   \\left[ \\begin{array}{c}\n",
    "                x \\\\ y\n",
    "           \\end{array} \n",
    "  \\right] =\n",
    "  x \\left[ \\begin{array}{c}\n",
    "                2\\\\ 1\n",
    "           \\end{array} \n",
    "  \\right] +\n",
    "  y \\left[ \\begin{array}{c}\n",
    "                -1\\\\ 1\n",
    "           \\end{array} \n",
    "  \\right].\n",
    "$$\n",
    "We can get a visual feel for this transformation by looking at a regular grid of points before and after the transformation:\n",
    "\n",
    "![](img/grid-original-small.png)  ![](img/grid-transformed-small.png)\n",
    "\n",
    "Let's look at it in action!\n",
    "\n",
    "![](img/stepwise-transform.gif)\n",
    "\n",
    "This notebook describes how to use `Numpy` to perform and visualize special linear transformations, the\n",
    "**scaling**, **shearing**, **rotation** and **reflexion** of a plane. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary libraries:\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------------------------------------\n",
    "## 1. General Steps ##\n",
    "\n",
    "    [1]. Create a rectangular array of points in x-y space.\n",
    "    [2]. Map grid coordinates to colors that uniquely identify each point.\n",
    "    [3]. Experiment and plot linear transforms.\n",
    "    [4]. Create animations\n",
    "       [4.1] Generate a series of intermediate transforms that will “smoothly” transition from the original grid to the transformed grid.\n",
    "       [4.2] Plot each of the intermediate transforms and save them as individual images.\n",
    "       [4.3] Stitch images into a gif to generate an animation.\n",
    "_______________________________________    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Generate original and transformed grids ###\n",
    "\n",
    "To construct a grid of points, we generate evenly spaced vectors along the $x$ and $y$\n",
    "axes, and combine them together into a grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a grid of points in x-y space \n",
    "xvals = np.linspace(-4, 4, 9)\n",
    "yvals = np.linspace(-3, 3, 7)\n",
    "xygrid = np.column_stack([[x, y] for x in xvals for y in yvals])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the $x$-axis values span from -4 to 4 and the $y$-axis value span from -3 to 3.\n",
    "By stacking the $x$-$y$ pairs columnwise, we generate a $2\\times 63$ rectangular grid of points:\n",
    "\n",
    "$$\n",
    "  \\mathrm{xygrid} = \\left[ \\begin{array}{cccccccccccccc}\n",
    "  -4& -4& -4& -4& -4& -4& -4& -3& -3& \\cdots &  0&  1&  \\cdots& \n",
    "   4&  4 \\\\\n",
    "   -3& -2& -1&  0&  1&  2&  3& -3& -2& \\cdots & 3 & -3 & \\cdots &\n",
    "   2 & 3\n",
    "             \\end{array}\n",
    "             \\right].\n",
    "$$\n",
    "\n",
    "We define a matrix \n",
    "$A = \\left[ \\begin{array}{cc}\n",
    "              2 & -1 \\\\ 1 & 1\n",
    "            \\end{array}\n",
    "      \\right]\n",
    "$ \n",
    "and generate the transformed grid by\n",
    "multiplying the column vectors of $\\mathrm{xygrid}$ by $A$:\n",
    "$$\n",
    "  \\mathrm{uvgrid} = A\\cdot \\mathrm{xygrid}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2 -1]\n",
      " [ 1  1]]\n"
     ]
    }
   ],
   "source": [
    "# Apply linear transform\n",
    "A = np.column_stack([[2, 1], [-1, 1]])\n",
    "print(A)\n",
    "uvgrid = np.dot(A, xygrid)  # A @ xygrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Plot grids ###\n",
    "\n",
    "To plot the grid points, we will use the matplotlib function `scatter` that can apply a different color to each point. The following function transforms an $(x,y)$\n",
    "coordinate pair to an rgb color:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function assigns a unique color based on position\n",
    "def colorizer(x, y):\n",
    "    \"\"\"\n",
    "    Map x-y coordinates to a rgb color\n",
    "    \"\"\"\n",
    "    r = min(1, 1-y/3)\n",
    "    g = min(1, 1+y/3)\n",
    "    b = 1/4 + x/16\n",
    "    return (r, g, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We map this function to the $x$-$y$ coordinates to generate an array of rgb color, \n",
    "and then plot the $x$-$y$ grid points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Original grid in x-y space')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQYAAAEICAYAAAC9P1pMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX5+PFPICxCQgKEPUBABGJkERDBWhkFFCutIiooomix1qJVkVa/YpugBQSiguJaN0QBtaLsEYGMCzvKFvYlCVkgAUKWyZ6Z5/fHhCHJTSBkTgjwe96+8npxz9x7nnNzZ545dzKex0dEBKWUKqFWTQ9AKXXx0cSglLLQxKCUstDEoJSy0MSglLLQxKCUstDE4KUpU6YwduxY4/uei4+PDwcPHjTS15EjR/Dz88PpdJb7eEREBA8++GC5j/3888906dLFyDjUxcNHv8dwxqeffsprr73GoUOHaNSoEcOGDWPq1KkEBgbW9NAsfHx8OHDgAJ06dar2WBERERw8eJDPP/+82mOpi4POGIq99tprPP/888yYMYOMjAw2bNhAfHw8gwcPpqCgoNxjioqKLvAozbsczkGZp4kByMzMJDw8nLfeeoshQ4ZQp04dQkJC+Oqrr4iPj/e8U0ZERHDPPffw4IMP0qhRIz799FPLNPuzzz6jffv2NG3alFdeeYWQkBBWrVrlOf70vnFxcfj4+DBnzhzatWtHUFAQkydP9vSzadMm+vfvT2BgIK1ateLJJ5+sMEGVFRsby0033YS/vz+DBg1i3LhxlrgfffQR7dq145ZbbvG0nU4SsbGxDBgwAH9/fwYPHsyJEycqjGW32wkODvZsh4SEEBkZSffu3QkICGDEiBHk5eWVe+wTTzzBPffc49l+/vnnGThwIOVNYmfMmMHw4cNLtT311FM888wz5fY9bdo02rRpg7+/P126dGH16tXAmWs4YsQI/P396dWrF9u3b/cc9+qrr3LllVfi7+/P1Vdfzbfffluq3//+97+EhoZ6Hv/tt98ASE5OZvjw4TRr1owOHTrw5ptvVvg7uySIkhUrVkjt2rWlsLDQ8thDDz0kI0eOFBGR8PBw8fX1lW+//VacTqfk5ORIeHi4jBo1SkREdu3aJQ0bNpSff/5Z8vPz5bnnnhNfX1/54YcfPMef3jc2NlYAGTt2rOTk5Mi2bdukbt26snv3bhER2bJli6xfv14KCwslNjZWunbtKm+88YZnXIAcOHCg3PPp16+fPPfcc5Kfny8///yz+Pv7W+KOHj1aHA6H5OTkeNpOn3+/fv3k2Weflby8PPnxxx/Fz8/Pc3xZ0dHR0qZNG892+/bt5brrrpOkpCQ5efKkdO3aVd59991yj83OzparrrpKPvnkE/npp5+kadOmkpCQUO6+ycnJ0qBBAzl16pSIiBQWFkqzZs1ky5Ytln337t0rwcHBkpSU5DnngwcPeq6Br6+vfP3111JQUCAzZsyQkJAQKSgoEBGRr776SpKSksTpdMqCBQukQYMGkpyc7HmsdevWsmnTJnG5XHLgwAGJi4sTp9MpvXr1kkmTJkl+fr4cOnRIOnToIFFRUeWey6VAE4OIzJ07V1q0aFHuY88//7wMGjRIRNxPqt///velHi/5Yp80aZIniYi4n/h16tQ5a2Io+UK47rrrZP78+eWO44033pC77rrLs11RYoiPj5fatWtLdna2p23UqFGWuIcOHfI8XjIxnD7e4XB4Hr///vvPKzHMnTvXs/2Pf/xDHn/88XKPFRHZuHGjNG7cWNq1ayfz5s2rcD8RkSFDhsgHH3wgIiJLliyR0NDQcvc7cOCANGvWTH744QfPC/608PBwuf766z3bTqdTWrZsKT/99FO5ffXo0UO+++47ERG59dZbZebMmZZ9NmzYIG3bti3VNmXKFBkzZsxZz+diprcSQFBQECdOnCj3fvvo0aMEBQV5ttu2bVthP8nJyaUeb9CgAU2bNj1r7JYtW5ba3+FwALB//36GDh1Ky5YtadSoES+++OJZp/Qlx9CkSRMaNGhw1jFXdB7Jyck0btyYhg0betrat29/zrglVXRO5enbty8dO3ZERLjvvvs87bfffjt+fn74+fnxxRdfAPDwww97bus+//xzRo8eXW6fnTp1YubMmURERNC8eXNGjhxJcnKy5/GS516rVi2Cg4M9j3/22Wf07NmTwMBAAgMDiYmJ8fzeExISuPLKKy3x4uPjSU5O9hwTGBjIlClTSElJOefv6mKliQHo378/9erVY+HChaXas7OzWbFiBQMHDvS0+fj4VNhPq1atSExM9Gzn5uZy8uTJKo3piSeeoGvXrhw4cIDMzEymTJlS7r13eWNIS0sjJyfH05aQkGDZr6LzaNWqFadOnSI7O9vTduTIkSqcQeW8/fbb5Ofn07p1a6ZPn+5pX7FiBQ6HA4fDwahRowC466672LFjBzExMSxdutTTXp4HHniAX375hfj4eHx8fHj++ec9j5X8fbhcLhITE2ndujXx8fE89thjzJ49m5MnT5Kens4111zj+b23bduWQ4cOWWK1bduWDh06kJ6e7vnJyspi+fLlXv9+aoomBiAgIIDw8HCeeuopoqKiKCwsJC4ujnvvvZfg4OAK35nKuueee1iyZAnr1q2joKCA8PDwSr2Yy5OVlUWjRo3w8/Nj7969vPvuu5U6rn379vTp04eIiAgKCgpYv349S5YsqXTc08eHh4dTUFDAL7/8cl7Hn4/9+/fz0ksv8fnnnzN37lymT5/Otm3bKty/fv363HPPPTzwwAP07duXdu3albvfvn37WLNmDfn5+dSvX58rrriC2rVrex7/9ddfWbhwIUVFRcycOZN69erRr18/srOz8fHxoVmzZgB88sknxMTEeI4bO3YskZGR/Prrr4gIBw8eJD4+nr59+9KoUSOmTZtGbm4uTqeTmJgYNm/ebOg3deFpYij2z3/+kylTpjBhwgQaNWrE9ddfT9u2bVm9ejX16tWrVB9hYWG89dZbjBw5klatWuHv70/z5s0rfXxJkZGRzJs3D39/fx577DFGjBhR6WO/+OIL1q9fT9OmTXnppZcYMWLEeY1h3rx5bNy4kSZNmjBp0iQeeuih8x7/uRQVFfHggw/y/PPP06NHD6666iqmTJnC6NGjyc/Pr/C4hx9+mJ07d541Wefn5/PCCy8QFBREy5YtSU1NZcqUKZ7H77zzTr788ksaN27M3LlzWbhwIXXq1OHqq6/mueeeo3///rRo0YKdO3fyu9/9znPcvffey8SJE3nggQfw9/fnrrvuIi0tjdq1a7NkyRK2bdtGhw4dCAoKYuzYsWRkZJj5ZdUA/YJTNXI4HAQGBnLgwAE6dOhQY+MYMWIEXbt2ZdKkSTU2BlOOHDlC165dOXbsGI0aNTrv4/XLWpWjMwbDlixZQk5ODtnZ2UyYMIFu3boREhJyQcewefNmDh06hMvlIioqikWLFnHXXXdd0DFUB5fLxeuvv87IkSOrlBRU5fnW9AAuN4sWLWL06NGICH369GHBggVn/cCyOhw7doy7776bkydPEhwczLvvvsu11157QcdgWnZ2Ni1atKB9+/ZERUXV9HAue3oroZSy0FsJpZRFjdxKBAUFVeq+Ozs7u9QXbS5lei4Xn8vlPKDy5xIXF1epL8rVyFeie/fuXan9oqOjq3cgF5Cey8XncjkPkcqfS2Vfe3oroZSy0MSglLLQxKCUsjCWGJxOJ9deey1Dhw411aVSqoYYSwyzZs0iNDTUVHdKqRpkJDEkJiaybNkyYysgK6VqlpHvMTzzzDNMnz6drKysCvf54IMP+OCDDwB3IrHb7efs1+FwVGq/S4Gey8XncjkPqIZz8eJPpyLiXmLriSeeEBH331LvuOOOcx6j32O4tF0u53K5nIfIRfg9hrVr17J48WJCQkIYOXIka9asqbA4iVLq0uB1Ypg6dSqJiYnExcWxYMECbrnlFv1/3ZW6xOn3GJRSFkb/JyqbzYbNZjPZpVKqBuiMQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSll4XViyMvLo2/fvvTo0YOwsDDCw8O9HtQRjvB3/s5+9jORiaSR5nWfZeWRxwzeZDDDeJQn2cku4zEAFhLFnxjLIeL5jpXVEmM7+3iEcAbzV17jM/IpMB7jBBn8H/9lIONJ4DiJHDcew4WL94jiNiZxP6+xjr3GYwCsYTf38jYHSeFTfqmWGIc4wd/4HwN5hwiiyCTPeIxsCpjMTwziM+LJYC+VqGJdWV4sTCsiIi6XS7KyskREpKCgQPr27Svr16+v8kq1J+SEtJbWgiCR0ZGCIN2kmxRKobdDLeVPcr8ggZ6fhtJGdskeozHels8ECREkRCKj3xQkRN6XL4zG2CH75QrpJ0hPz8/dMt5ojHwpkKvlYUFsgtgkMvp9CZZ7JU0yjcYZJ+8Lcpfnx1eGS7TsNBpjkfwmPvKIIGMkMvoTQcbIS/KN0RjJki5BMlGQZzw/feV1cYnLaJxbZI4gEYJESGT0PAmQqXJY0s56zAVbJdrHxwc/Pz8ACgsLKSwsxMfHp8r9fcZnJJNcqm0nO1nGMq/GWdIe9rGYFaXassnmLT4wFgPgVd6rVJs33mQ+uWXejRayhgPEG4uxhHXsLtNfIsf5nB+MxUgjiw/K9FeEk0i+MxYDYBrLEaRU20x+II9CYzE+ZCMnyC7VtokjrOGAsRibSGINsaXaMsjnXbYY6d9HROTcu52d0+mkd+/eHDx4kHHjxjFt2jTLPmULzixYsKDcvpJI4hjHAAh2BJPolwhAO9rRjGbeDhWALBzs56ClPZAArqSDkRgAvxHjeRIGO1qQ6JeCD7XoRZixGAdJIAOHpb0L7fGjgZEYx0nnCKme7WBHExL90mhJE9oQZCRGPoXEcMTS3oB6hBJsJAZADEnkFyeBYIc/iX7uIkk9aIevoY/cjnCK42USA0AHmtDE0DXJII+DnPJsBzvqkeiXTxOuoAOBFR43YcIEtmypRPLwfkJzxqlTp8Rms8nOnWef/p1tOmMXu1D83+lbiVpSS2Il1tg4cyRHGktIqVsJJFA+kE+NxRAR+ZOMtdxKDJe/Go3xjnxZ6jYC6SlNxSa5kmcsxgFJlFpyS6lbCcQmv8gOYzFc4pLO8rdStxLIXfIvw7deT8pcQcaUupW4TiYZjbFEYkrdRiDPSF15To4ZvPU6JbnSUCaXupVAImT+OW69LtitREmBgYHYbDaioqKq3McABhBOOHWpC4AffrzHe4QQYmiUcAVXMJ8PCaIpALWoxSOM4lHMFsp5m5e5tsTsoDfdeJMIozEe424eYig+uG/fmtOE+UylPvWMxehEG97maRpSHwAffHiZR/gd3YzF8MGHL3iWdiVmhUPpwwsMNxYD4BXu5hbOFF++ihZ8yp+NxhhKGBO42TMDaUR9PuZ+WuBvLEYg9ZnLMBoXXxOAv9GHEaZmo95mrtTUVDl16pSIiOTk5MiNN94oS5Ys8TprpUqqLI1eKhmS4e0QK5QnebJeNskRSai2GCIi22W3LI9eUa0x4iVZNsgOyZeCaouRLlmyTmLkh+hV1RajSIpkk+yXA5JcbTFERPZIsiyL/t74B4IlHZNMWSex4jA4eysrRwpknRyRldGrK7V/ZWcMXteVOHr0KA8//DBOpxOXy8V9993H0KFDvU5YzWhGQxrSiEZe91WRetSjH9dVW/+ndSeUNFKqNUY7WtGOVtUaIwA/+hOGvRr+VHlabWpzHVdVW/+ndaUVx9jnmWlVhxb4G50llOcK6tCfttg5ZLRfrxND9+7d2bp1q4mxKKUuEvrNR6WUhSYGpZSFJgallIUmBqWUhSYGpZSFJgallIUmBqWUhSYGpZSFJgallIUmBqWUhSYGpZSFJgallIUmBqWUhSYGpZSFJgallIUmBqWUhdeJISEhgZtvvpnQ0FDCwsKYNWuW14Ny4eI3lpNBCnv42ev+KhJPDIt5i3UspKgairQAZHCC75lLJmlkcrJaYhRSwFqiWMxnHCln9WtTdrONxSwgl2xcuKolxjGOspTF2FlDHrnVEiOHHFbzI5lkcqKarokLF5vZyWLWsJ+4aokBcJhkFrOWbPIowmmsX69XcPL19eW1116jV69eZGVl0bt3bwYPHszVV19dpf6KKGQaQ9nOSoYQSTgjuZlHeYKPvB1qKd8xk3lM8my3JZRJLMfvLEtvn6+9bGYKj5BHNkMYz1PYmMgcOtPLWIws0gnnzyRyGIDPeYMHeZY/8ZCxGABv8yo/8j0AQ7ifabzIP5lMbWobi/ETdt7hLU/S+ZL5vMxkmtHcWIxEkpnENDLIZAgDeIrnGc/fuM7gNSmgkMm8x57iawLwB25iDMOMxQBYwBoWFr9xDqEz/+JjwnmY+sULKXvD6xlDq1at6NXL/Uv19/cnNDSUpKSkKve3kW/YXqaUWzQfc4CNXo2zpHRS+ZIppdoS2MMK3jcWA2AO/yGvRH2BXBx8xmSjMZYzz5MUTvuSt8ksUXPAW3vZ6UkKp21lI5sMzuYKKWQOH5eaiZzkBN/wtbEYAF/yLRlkeradOPmUeUZnQD+xpVRSAFjOTyRw1FiMVNL5tkx5vUMks4pfjfRv9DOGuLg4tm7dyvXXX1/lPg5VUEnnsKETBjjCLpzlVB46zDZjMdz9xVSqzbsYuy1thRSQYHBx0MPsq6B9v7EYqaSQRVY5Mcwucnq4TPUmgOOcJKucoj1Vj5FQQXuisRixHLVU1HLHMJN8jFSiAnA4HAwYMICJEydy9913Wx6vbCUqB2mkFl+8AEcwGcWVqFrTlfo0NDFUnBRyhN2WX2xjWtKYlkZiACRxkPzi++QARwsy/FKoRwPacKWxGGmkkl7mPtkHaEdnY9P8PHJJLvFkD3A0IcMvjea0ws/QKsiCEE+c5Z3bH3+jtxIppJJNDgABDn8y/LLwxZd2BqtdZZLNiXJmbMG0oC51jMQowlmqOliAox4Zfvk0pREBZ3mdVLYSlZHEUFhYyNChQ7ntttsYP378Offv06dPhYMrooBXGMQefmaIPZIo2wRuYCTPMN/bYZbyFVP5H9M92624klf4nkbFRWhM2MlaXuXPFJLPEPt41tje4f/4hDD6GYuRQRovMYaUEi/c+/gr9/C4sRgArxPBBn4EYIj9fhJse5jIDHy9/5jKYxUr+S/veRJ2AIG8whRaGlwWP54EIphGNtkMsQ9gpe1n/s5fuIGqz3LLyqeACGZzqMQ1uYV+/JURxmIAzOF7lrEBgCH2zuyypfMKj9CgRBGass722ivJ68QgIjz88MM0adKEmTNnVuqYcw2uiEI28S0p9iJCbI3pyZBqWf9/P5vYxhqCaMMN3G1sRlLSCZJYyxKwB/E72w0E0dp4jDxyWMf3nCSVntzAVQYrRJ3mwsU2NnKQvTS0N2OI7XajHzyedoR4NrERP/y4kd8bm5GUlEEma9lAgT2Xvra+tK6GehyFFLGR7SSTSlc60p0uxmMA7CaeGA7TwJ7HrbZB55yRVDYxeF2J6ueffxZAunXrJj169JAePXrIsmXLznpMZavhREdHezu8i4aey8XncjkPkcqfywWrRHXjjTciZj6mUEpdJPSbj0opC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopi4s2MRSSiZNcnORVWwzBSQa7yOdEtcUAcBCPq5oK2pyWRxoZHEAMFh0py0k+6RxCqqnYDLgXhc0kgVzSqi0GQA4ncZWzUrhJBeSSTjJOiqothgsn6aQYvyZGVvJ89NFHWbp0Kc2bNycmxvvl0fcxg11E4MPLLOVBruUt2vGAgZGecZy1bOQv5JKMD3W4kjH0ZJrRtSXzOM56/sYJtuDLeKL5iP68TX2CjMUQXGzlVQ7zLUIRDWhJXybTzGABFYB4VrKV2RSSRR0eYD8n6cy9RmNkEM96ZpBFEuBDW37Hdfyd2oZWVgYoIo8NvM9RtuHLAFbyb25gHH60MBYDYBdR7GUNLgqphx+9uZc2dDca4xgH2cwi8nBQh87s4Ae6M9hI30ZmDGPGjCEqKspEV5xgHTv4J87iJb4LSGMTD5NTwVr9VeEkn/WMIZdkAIRCDvJf4g2vRL2VcE6UqJNxgk1s42WjMWL5jkN8jRS/K+VwjPX80+i7YTbH2Mw0CovrPgjCdt4hjT3GYgBs4LXipOCOksAv7ONbozFi+JajJeqHZJDAJj40GuMYe9nN955rkI+DDXxOvsHaFYXks56vySvuU4B9rCOxnDojVWEkMdx00000adLERFcks9jSJhRxlOVG+gc4yaZybx+SDMZw97fK0pbMD0ZjJBcv6V5SPic5yU5jMY6yvtypahJrjcVwcJTMcpJ/MpuMxXD3Zy0qdJKDRl+0yeUUFXJRyLEKCvdUxXHiyq23mmwohrGCM3FxcQwdOrTCW4nKFpzJJ4Wc4oo9Po5gpLjgjB9XUsdQXUkneWSy19Jej6Y0oK2RGAAZ7PW8a/g4WiB+KdSiLgEGlxLP4Sj5pFvaA7iSWgZqGAIUkoWjRIUjH0cTxC+NBjSjHo2NxBBcpBMHZYoA1aEBfgaXd8/iKEXkA+Dj8Ef8svChFoG0BUO3kXlkkluiDN5p/gThe5aaD+ejiEKySry5+TjqIX751KchV9CowuMqW3DGXLWQc/jLX/7CX/7yF8C9tr3NZit3vzxSWck15HMcX3skRbYJ+NGZgew09kQH+In/ksIaz3Yt6nMDKwk0WJNhPwfZzjQAfO3jKbK9Tk/CuQqbsRin2MMaxpT6cLMlN/J7njAWw0kBK3kUR/E039d+P7VtUdj4hHoEGIuzhbeJLTXLqkV//kVLehqLkcBmNvAOAL72ARTZfqQzQ+jBzcZi5HCKlcygsES17gBaM5AH8DH4ef8aPuZk8SzL194ZscVyI4/jb6Bo0gVLDJVVn+bcwjr2MJVjNKQDT9OV/zOaFABu4DP28RbHWEMD2tCFp4wmBYDOjKUeTYjjG7LwozczacedRmM0JhQbH7KfueSSQktuoLPhSte1qYuNWezlC9LYSz4B3MRso0kBoBd/xZ9gkthAXfy5iqG0MPyBXVuuozZPc4g1pFOfbjzIlQaTAkADGnMLT7OXVWRxnCA60pWBRpMCwO95gL2sJZVY8riC3/GokaQAeF9w5rTY2FgJCwur1L5acObSdrmcy+VyHiLmC84YSWH3338//fv3Z9++fQQHB/PRRx+Z6FYpVUOM3ErMn2/2z3xKqZp10X7zUSlVczQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLIwkhqioKLp06UKnTp149dVXDfS4HhgE7ACGQTkrOnsvFRgHXA0MBsNLx7s5gbcAG7AfeBuqpYpTNDAauBWIgGqp4hQHTASGA4lguKaEWy7wITAWeBZKLNZr1mrc55IIfAXVUpHqIDAb+DfwKVRLtbNMYCHwGnCyOKYZXicGp9PJuHHjWLFiBbt372b+/Pns3u1N0Yt43ElhNe4L9h1wM5Dt7VDLGAHMA44CW3C/sNYbjvF68U88kA9EArMMx9gCPAfEAMeBRcAzhmPkAOOBdbiTTk5xzOOG48wCluF+ksfhTqq/GI4RDXwOJANFwArgM8MxUoD3gUNAFrAdeLc4nikCzAV2Ag6gAJgPJZb594bXiWHTpk106tSJjh07UrduXUaOHMmiRYu86HEuFFehOuMY7ie8KVvBUnjEBcwxGAPciaesLwzH+IaytRjcTxZzxU3gZ+BUmbZcMFo8Jw3KLS7zvcEY4E4MZa0HozVSN2FNAmmYnfnGYZ2FuIDfjPTu9ZqPSUlJtG17pkhLcHAwGzdutOxXtuCM3W6voMeOuN9ZweEIxm6PLG6vB1R0zPnKBsLLaQ80GAPcU2L3rYPD0QK7fTxQ23CMMKBdOe2HMfXu4X7R3O/ZcjiaYLffj9lrUgTcUU57fYMxALrgfo6Bw+GP3T4Ad6GZdZj7yM0J9C6nPQFztxT5QGfPlsNRD7u9M+7fo93r3r1ODFJOISsfH2tFn8oWnHF/rjAacGG3R2KzTQAa4J6OmyoG6wReBI6UaZ8LBovBuKepXwFgt4/HZnsdeAD4m8EYP+C+xywpGPcMy9QTPQ33uN0VnOz2+7HZvgQ+AkIMxQB4AetM58+YvSbfcXr2abcPwGb7EegGPGwwRhzwJqVnclfg/rzBTCUqdwKYyelbbLu9MzbbftyvnY5e9+71Myc4OJiEhDM1BxMTE2ndurUXPXbHfQ94ehbSFViMuaQA7nft+UCf4u0mwMvAUIMxAP4F3Ik7//rg/iB1ouEYg4G/g6f4S3fcTxiTf3BqAvyHM9fEF3gJs0kBYALQA/fvqj7u39cfDMcYCgwE6hbH6Yl7ZmdSCO7PsE6XimtVHMNUUgD3NXiguG9wX+/bMZEUAO8LzhQWFkqHDh3k8OHDkp+fL927d5eYmBgDRS+cEh29ytvhVUKGiBRWc4xciY5eU80xCkUkq5pjiIhkXYBCLTkiUlDNMQouwDVxivtcqluu8YIzXt9K+Pr6Mnv2bG677TacTiePPvooYWFhBlJWLdzv7NWt4gKg5tTHVMHUivkCftUcgwsU44oLEKMO1X9NanFhzsXkTMTNSMGZP/zhD/zhD6anfEqpmqLffFRKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWXiVGL7++mvCwsKoVasWW7ZsMTUmKMyAfZGQEweH3genyaW9i4kLjnwDm5+EXVMh19SKymWc3Aq//RtykyFte/XEyE2FmHdgcwQkrIRyFuj1mjMfDi2CzdMg/xQUmq7zUezor7DlbdgxBxzVdE0yk2DbfMg5Cane1EA5i/xs2LMKNs+H2A3gcppRDf7TAAATKElEQVSPIS6I3w6bF0GeA/LMXROvVnC65pprWLhwIY8//rip8UCRA1b3h6w94BsJv02AhAUwYA2Us/p0lf36LMSWKDRy6BMYuAoalrcUexUlLIUNTwMu8B0Pq++G/m9D8BBzMbKTYNWDUJDu3o5bBB03Q2+Di86KC376J5zY4d72bQxrnoKB74CvwWXFdn8Fu+af2T4YBTdPhsaGFjgFOL4PfooEVxH4DoAfZ0DPUXDVIHMxCnJg9RuQfdK9HbcJknbCjY+ZiwGw8VtIiHH/27czrPoABj0G9b1ffs+rGUNoaChdunTxehClxH/hTgolHbdDisHiJo44iJ1bui3/OBz8wFwMgJjXKV2SzgW73jAbY/8XZ5LCaYcXQo7Bd9tjm84khdMyYyFhtbkYhbmwd2HpNmce7P3GXAyA3YvdSaFU2yJrmzdiN55JCqcd3Q0n48zFyEg5kxROy82EQ2Zm7j4i3s87bTYbkZGR9OnTp8J9yhacWbBgQfk75iZCXgoADp9g/CTR3d6gHdRr5u1Q3YockFVOnb86AeDXwUwMgFMxnK4t4PBpgZ+kgE8tCDSxWG4xRwIUOqzt/u3Bt4GZGPnpkJN6JqRPE/wkDeo3gSsMLevvKoSMsnU+AN964B9sJgZARqInCTh8/PGTLHd7YDv3tTEh55T7VqKshk2grqFrUpgHjjPVwRw+9fCTfKh7BTQMrPCwCRMmVOq2/5y3EoMGDeLYsWOW9smTJ3PnnXeeM8BplS44k7wE1o4GwO4bia1ogru9XwwEGHpBFWbAksfBWaYUXth/oHMF46oK+wdw3F0P0+47HlvR69Di93DTOHMx9n4CO8vMdHwbwICVUKehmRjpB+GHM9Ngu+/92Irmw7XToGVfMzHECcseh9wy77Qd/gg9bWZiAGz6L8SvA8DuOwBb0Y/g3xpuHmMuxpHfYGOZGSk+8LuX3MnBhDwHLH3DfZsH2H07YyvaDz2GQsfyqmCdn3OmyFWrVhETE2P5OZ+kcF5aDYWQR0sP8eoIc0kB3DODXq9Brbpn2prb4MpHKzykSq6NgPotzmxf0Qp6llcazwudRkKz685s16oLvSaaSwoAgZ0gdDSlllvvcAe0uK7CQ86bT23oM670ZxaBHSH0HnMxALrdA/6tzmzXaQh9xpiN0bYntL32zLZPLej+R3NJAdyfI1w7pPQsp1VnCOlppv8q17goYcCAAbJ58+ZK71+pohfpOyV65UIRx2EvRnYOuakiRxaKnNhUfTGK8kSSVkn0ysUizvzqi3Nim8iRlSJ5J6svRlaiyJHVEr3q++qLUeAQSVgnkrJDxOWqnhgup8ixGIleuUykMK96YoiInEoUOfKbSHZa9cXIyRQ5slOiV62s1O6VLTjj1U3Vt99+S3BwMOvXr+eOO+7gtttuM5OtAAKugTqNoaHBe/6y6jeDtsOgqcF3vrJq14PWA6GOf+kZimlNe0DbwVDP4LtSWX5toO0tULsaz6NOQwjuD827mf0rVEk+taBFGNRp4P4Mo7oEtnHPHBo0rr4YV/hD22ugdh2j3Xr158phw4YxbNgwU2NRSl0k9JuPSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopC00MSikLrxLDP/7xD7p27Ur37t0ZNmwY6enp5z6oso7shOxTkHLYXJ9lpafCuoWwb2P1xSjMhy2rICfL/e/qsncbrPseMtKqL8axJFi7BgoKqi9GdjZsWAc7d4DLde79q8Llgm07ITsH8qqhmNFphxNg7W9w4tS5962qtExYGwP5hUa79SoxDB48mJiYGHbs2EHnzp2ZOnWq9yMSgXcegQnd3Unh753gqwjv+y3LPh/+1g1m/hn+NQReuQsKcs3GSNgPTw6A6Y/BsTh4ygaJ5Sxb7438XJj0GPz7EZj5AjwxBH5aZjYGwIKP4OnRMOsVSIyD9yLNV7za9hs8MRZenwGvhMMLEyAz02yMk2kwfiJMeR1SUuGJCbB3v9kYLhe8MQdeeB1mzYUn/wOL15iNAbByM4ybBbO+gaQT8Op8KDJT8cqrxHDrrbfi6+teHa5fv34kJiZ6P6Jfl4L90zPbIvC/SZCwy/u+T8vOgI8mQFGJd76dP8LKj83FAPg4Ak6lnNlOOwZzXjYbY8UC2FWiTkBRIXw4BXINlpCLOwgLPy+dCNYsg+2bzcVwOuG9d0q/g8fFwsL/mYsBMP8bSC5RDiE7G96fYzbGum2wftuZbZcLvlgKqQZnc+kO+CQKnCVmVb/tB/u2io85D8Y+Y/j444+5/fbbve9o94/n114Vh7ZCfo61fdcv5mIA7C7nFmWX4duWPb9a2/Jy4PAea3tV7a6g5uYuM09CAFKOQdpJa/vuGGubN/bss7YlJUOGwZnJ7nJmhSKw55C5GPsTSicFT+x4I92fsxJVZQrOTJ48mS1btrBw4UJ8KljZt9KVqLJOwHH3yTkCgvHLKJ6FtOrsXhHXhKICOFJOMdOA5tC0tZkY4L6VKP5cwRHQAr+MFKhbH4KvMhfjZEr5nyu06wS+hlYOzsl2f75QzBHQBL+MNGjWAvwDzMQQcc8Qyj4dG/pBixblH1MVR1Mg133L6Ajwxy8jC2rXhvZtzcXIyIKTGdb2Ns2hnqEVtguKIPG4Z9MRUA+/jHxo7A+NK65dWdlKVF6XqJszZw7vvfceq1evpkGDypXf6tOnT8WDy3PAi9dD4m7sQyKxRU2AsJvh36vNLif+/jOwusQUMqAZTF0NQQafIOuXwcy/g7iwDxmPbeUseHY2XG+wqG1qErw4GjJLfMB1670w9kVzMVwumPQs7HHXr7QPuR/brg0w5R2oZ7Co7f++gq9KFLWtXx8mTYYOBova7tkP/5kBhUXYhwzAFvUjPDoKhhgsauvIgf97A1JOnGnrdTW8YLio7cxvYJ17RmUf0hnb5mMw9TEIrDgxnPW1V4JXiSEqKorx48fz448/0qxZ5etKnnNwORmw+kPseUHYGufDgIehjuH1/0Vg/XewfQ00bQODx0DjlmZjABzcDj8txB7QFdu110DHbuZjpKXCqoWQlgI9b4DrB5mvyVBQAPYoOLgHe5vO2AbfCg0MVrs6beuvsGkj+PnBwMHQstW5jzlficmw5ifsDepg63o1XBNqPoYjG1ZtgORUCO0Iv+8DvrXNxnC53H+R2BmLvXVtbDffDAFnvyaVTQxeVaK68sorJTg4WHr06CE9evSQxx9/vFLHVbYaTnR0tBeju7jouVx8LpfzEKn8uVT2tedVwZmDBw3/6U0pdVHQbz4qpSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSw0MSilLDQxKKUsNDEopSwu3sTgcrmXFK9uGRlQVFS9MXJyzddgKKuw0L0IaXUScS9bXt2ycqHAbAEVi7xCcFXzNXG6IN1wrZLypJsvmuNVYvjXv/5F9+7d6dmzJ7feeivJyclmRjVvHrRvD9u2QdeusGqVmX5L2rnLvTZiYAi06AzTZ5mPkeWAB56GgO6wdRc8+Ix7LUDTpn4ELQZB4E3Q/yHYZXCZ8tNWboYuD0HjP8LOw/BlNRRQOXIcBoVDowegyUPwzznmq1EVFsG4udB4HGw9AkNnQqrhojYAH22A4Aho/CJ0mw6/VENFtS3J0Ot9aDwNtqfA7E3Guva6RN2OHTvYtm0bQ4cO5eWXDRRT2bEDRo+G08Vr9u2DO++E48fPftz5cDrhj/fDpuKaDGmn4PkI+HapuRgAz7wC8xe7ZyQi8MUimDDFbIyvf4AXZ8Op4if3hp3wp2fMvqBS0uCuf8GB4mtSUASjJsOuWHMxAEZEwmr3StRk58GM72D2crMxpiyFd9a4ZwwILNsOYz40G2N9HIz9Eo4Vz+BijsIfP4Qsg+/seUVwxzzYWlzaocgFT62AVWYSkFeJoVGjRp5/Z2dnV1hT4rx89ZX1SZ2TA4sXe9/3aes2QnyCtX2e6apH5Yx5nsHzAJgfZW07nAgbd5qL8d0vkFum7qbTBV/ZzcWIS4UN5ZSKm/ezuRgA88sp+BMVA+nlFCCqqgW/WdvSc2G5wSJAa2IhtZzZ5zwz193ruhITJ07ks88+IyAggOjo6AqXka90wZljxyDJXdzEERyM3+mZQ4cO0KSJN0M9Iye3/IpETZtASDszMQC27/F8fuEIboFfYgrU8YXuBpcrj01yFzYt6+qOcIWhJfdPZkLcmaJDjuAm+CWmQZsgaGnomhQ5YXuctd3/CuhssAjQnmTIcZcmdAT745eY5V5qv2c7qGVoyf2kjDOzhZI6BUGAoTocWQWw/0zlLkdwPfwS86F5Q2jbqMLDKltw5pzLxw8cOFDCwsIsP999912p/aZMmSL//ve/vV/COi5OpGFDEZDoyEgREGnVSsThqFTfldbnZhECz/zUaiLyy3qzMSbOECFEhBCJjnzT/e/wN8zGsG8W8blWhJ5nfq4fbTZGZrZIi2Ei2ESwSXTk+yJ+t4skpJqNM2yqCHeV/vnyF7Mx3lsjwhgRxkh05Cfufz/yodkYe1NE6k0Q4ZkzPyEvi+QXmovhdImEzhYhQoQIiY6cJ+L7sshvyWc9rLLLx3tVV6KkuLg4CQsLq9S+5xzc+vUigwdL9KxZInffLbJ3r4ERlpGSKvLIOJE2V4v0GyyyeLn5GEVFIv95S+Qqm0TPfl9k6jsiTqf5ON+tcSeDNreK/DlCJPWk+Ri740Tuekmk9T0S/cFckY27zcfIyhF5+kORdo+J9HxW5NPV5mOIiLyzWiRsokS/OVfk+a9EcgvMx7AfELHNFmkdLjJijsjhE+ZjJGaIPLhQpPVrEv3e/0S+P3jOQy5IYti/f7/n32+++aYMHz68UsdpwZlL2+VyLpfLeYhcZAVnXnjhBfbt20etWrVo37497733njfdKaUuEl4lhm+++cbUOJRSF5GL95uPSqkao4lBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSlloYlBKWWhiUEpZaGJQSllYSQxREZG4uPjw4kTJ0x0p5SqYV4nhoSEBH744QfatTO47LpSqkZ5nRieffZZpk+fbqbYjFLqouDVmo+LFy+mTZs29OjR45z7li04Y7fbz3mMw+Go1H6XAj2Xi8/lch5QDedyrmWkz1Zwpm/fvpKeni4iIu3bt5fjx48bXcL6/8flvS8Fl8u5XC7nIVIDy8evqqDS9M6dO4mNjfXMFhITE+nVqxebNm2iZcuW5jKXUuqCq/KtRLdu3UhNTfVsh4SEsGXLFoKCgowMTClVc/R7DEopC68+fCwpLi7OVFdKqRqmMwallIUmBqWUhY+IyIUOGhQUREhIyDn3O378OM2aNav+AV0Aei4Xn8vlPKDy5xIXF1ep/3WhRhJDZfXp04ctW7bU9DCM0HO5+Fwu5wHmz0VvJZRSFpoYlFIWtSMiIiJqehBn07t375oegjF6Lhefy+U8wOy5XNSfMSilaobeSiilLDQxKKUsLpnEcDksH/ePf/yDrl270r17d4YNG0Z6enpND+m8REVF0aVLFzp16sSrr75a08OpsoSEBG6++WZCQ0MJCwtj1qxZNT0krzidTq699lqGDh1qrM9LIjFcLsvHDR48mJiYGHbs2EHnzp2ZOnVqTQ+p0pxOJ+PGjWPFihXs3r2b+fPns3v37poeVpX4+vry2muvsWfPHjZs2MDbb799yZ4LwKxZswgNDTXa5yWRGC6X5eNuvfVWfH3d/99av379SExMrOERVd6mTZvo1KkTHTt2pG7duowcOZJFixbV9LCqpFWrVvTq1QsAf39/QkNDSUpKquFRVU1iYiLLli1j7NixRvu96BPD+Swfdyn5+OOPuf3222t6GJWWlJRE27ZtPdvBwcGX7IuppLi4OLZu3cr1119f00OpkmeeeYbp06dTq5bZl7Kx/+3aG4MGDeLYsWOW9smTJzNlyhRWrlxZA6OqmrOdy5133un5t6+vL6NGjbrQw6uy8v6qfanP4BwOB8OHD2fmzJk0atSopodz3pYuXUrz5s3p3bu38bUrL4rEcDktH1fRuZw2Z84cli5dyurVqy+pF1ZwcDAJCQme7cTERFq3bl2DI/JOYWEhw4cPZ9SoUdx99901PZwqWbt2LYsXL2b58uXk5eWRmZnJgw8+yOeff+5951VZeLKmnM+CsxejFStWSGhoqKSmptb0UM5bYWGhdOjQQQ4fPiz5+fnSvXt3iYmJqelhVYnL5ZLRo0fL008/XdNDMSY6OlruuOMOY/1d9J8xXE6efPJJsrKyGDx4MD179uSvf/1rTQ+p0nx9fZk9eza33XYboaGh3HfffYSFhdX0sKpk7dq1zJ07lzVr1tCzZ0969uzJ8uXLa3pYFxX9SrRSykJnDEopC00MSikLTQxKKQtNDEopC00MSikLTQxKKQtNDEopi/8H1WbXlm4MVggAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Map grid coordinates to colors\n",
    "colors = list(map(colorizer, xygrid[0], xygrid[1]))\n",
    "\n",
    "# Plot grid points \n",
    "plt.figure(figsize=(4, 4), facecolor=\"w\")\n",
    "plt.scatter(xygrid[0], xygrid[1], s=36, c=colors, edgecolor=\"none\")\n",
    "# Set axis limits\n",
    "plt.grid(True)\n",
    "plt.axis(\"equal\")\n",
    "plt.title(\"Original grid in x-y space\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, We can plot the transformed grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Transformed grid in u-v space')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ0AAAEICAYAAABF36G7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xlc1NX+x/HXAOIKgiCbKIigKWqWuF5Tci/NtLSk7Nr1lsvP6ubNuvdqmnVN7Zr3WtpGppc2LcutNK9oomW5lmYuiUsCiiIiIPt2fn8MjiwzysicYfs8e8wj5zDz/RyG4c18v9/zPceglFIIIUQFOVR1B4QQNYuEhhDCKhIaQgirSGgIIawioSGEsIqEhhDCKhIadvTdd98RHBxMkyZN+Prrr6u6O6WcPHkSg8Fgs+3FxMQQGhpq8evjxo1jzpw5Zr8WFRXFPffcY7O+CNuq9aHRpEkT083BwYGGDRua7n/yySd27cuLL77ItGnTyMjIYPjw4XatbW/h4eEcOXLklp47fvx4vvnmGxv3SNiKU1V3QLeMjAzTvwMDA1m2bBkDBw60+PiCggKcnPS8LGfPnr3hX98b0dkvW6tJfRXWq/WfNG7mxRdf5OGHHyYiIgIXFxc+/vhjfvzxR3r27Imbmxu+vr4888wz5OfnA8ZfCIPBwHvvvUdwcDDu7u4888wzpu2dOHGCvn370rRpUzw9PXnkkUcAY2DFxcVxzz330KRJEwoLC0lISGD48OE0a9aMkJAQli9ffsN+vfjii4wdO5aIiAiaNGnC7bffzqlTp5g7dy7NmzenVatWbN261bSN1NRU/vSnP+Hr64u/vz+zZ8+mqKgIgMLCQqZNm4aHhwdt2rRh8+bNN3yd9u/fT5cuXXBxcWHs2LGMGTPGtHuxdetWAgMDmTdvHj4+Pjz55JOmtmsOHDhgen5ERAS5ubkWay1btozw8PAKvd5lld3tKduPkuLj42nUqBFpaWmmtn379uHl5UVBQUG5x+/evZs777wTV1dXvL29ef7554Hru3bvv/8+fn5++Pn58Z///Mf0vBu9nwAOHz7MwIEDadasGT4+PvzrX/8CoKioiHnz5tGmTRs8PT0ZO3YsV65csfi9242qQwICAlR0dHSptpkzZ6p69eqpDRs2qMLCQpWVlaX27t2rdu/erfLz89WpU6dUSEiIWrJkiVJKqfz8fAWoESNGqNTUVHXmzBnl7u5u2u7o0aPVggULVGFhocrOzlbff/+9qVaLFi3U9u3bTfd79+6tnnrqKZWdna0OHDigPDw8VExMjMV+zZw5UzVo0EBFR0er/Px8FRERoQIDA9X8+fNVfn6+evvtt1VwcLBp+8OGDVNTpkxRmZmZKjExUd15551q2bJlSimllixZojp06KDi4+NVcnKyuuuuu5Slt0NOTo5q0aKFWrJkicrLy1Off/65cnJyUi+99JJSSqno6Gjl6Oio/vGPf6jc3FyVlZWloqOjVUBAQKnnv/HGGyovL0+tXLmy1PPLev/991W/fv0q9HqX9eijj5babsl+mHPXXXep5cuXm+4/++yzaurUqWYfGxYWpj799FOllFLp6elq9+7dSimlYmNjFaDGjRunMjMz1cGDB1WzZs1MP+sbvZ9SU1OVl5eXWrx4scrJyVFpaWlqz549SimlFi5cqHr37q0SEhJUdna2mjBhgho3bpzF78VeJDRmzlR33333DZ+3cOFCNXr0aKXU9Tfxjz/+aPr6qFGj1MKFC5VSSkVERKjJkyerhISEctspGRqnT59WTk5OKiMjw/T16dOnqz//+c8W+zVz5kw1dOhQ0/01a9YoV1dXVVhYqJRSKiUlRQHq6tWrKiEhQTVo0EDl5OSYHv/hhx+qgQMHKqWMvyzvv/++6WsbN260GBrbtm1TLVu2LNXWo0ePUqFRv359lZuba/p6yV/Wbdu2KX9/f1VUVGT6erdu3awKDUuvd1nWhsY777yjBg0apJRSqrCwUPn6+qpdu3aZfWyvXr3Uyy+/rJKTk0u1XwuN2NhYU9u0adPUxIkTzW6n5Pvpww8/VF27djX7uODgYNMfEaWUiouLU87Ozqafd1Wp87snAC1btix1//jx4wwbNgwfHx9cXV2ZPXs2ycnJpR7j4+Nj+nejRo1Mx04WLVpEfn4+YWFhdOrUiaioKLM1z58/j6enJ40bNza1BQQEcO7cOYv9AvD29jb9u2HDhjRv3hwHBwfTfTAexzl79iy5ubl4e3vj5uaGm5sbU6dO5eLFi6b6JbcfEBBg8fU5f/48/v7+pdrK9s3b2xtnZ+cbPr/k2Zkb1TPH0uttjaioKNNB8Pvuuw+AMWPG8N1333Hx4kW2b99OgwYN6N27t9nnr1ixgqNHj9KuXTu6d+/Opk2bSn297Ot5/vx54Mbvp/j4eIKDg83Wi4uL47777jP9/Dp16oTBYCApKcnq792WJDSg3KnGSZMm0bFjR06ePEl6ejqvvPIKqoIXA/v6+rJs2TISExN56623mDhxImfOnCn3OD8/P5KTk8nMzDS1xcXF0aJFC4v9skbLli1p1KgRKSkppKamkpqaSnp6Or/88oupn/Hx8aVq3+h7SkhIKNVW8rk366u559+oXmU0btyYrKws0/0LFy6Y/j1+/HgyMjLIyMjgq6++AsDDw4P+/fuzevVqPv30UyIiIixuu127dqxatYqkpCSee+45HnzwQXJyckxfL/t6+vn5ATd+P7Vs2ZJTp06Zrefv7090dLTp55eamkpOTk6pAK0KEhpmXL16laZNm9K4cWOOHTvGe++9V+Hnfv7556ZPC25ubhgMBhwdHcs9rnXr1oSFhTFjxgxyc3M5ePAgK1as4NFHH7XJ99CyZUv69evH9OnTSU9Pp6ioiJMnT7Jz504AHnroIRYvXsy5c+e4fPkyr732msVt9enTh8LCQt555x0KCgr48ssvOXDgQIX70qdPH4qKili6dCkFBQWsXr2an376qdLfozldunRh48aNXLlyhcTERN58882bPueRRx4hKiqKNWvWmA5cm/PRRx+RnJyMg4MDTZs2xWAwmD7lAfzzn/8kOzubw4cPExUVxcMPPwzc+P00YsQI4uLiWLp0KXl5eaSnp7N3714AJk+ezIwZM0wBm5SUxIYNG27pdbElCQ0zFi1aRFRUFC4uLkyaNMn0w6+IPXv20K1bNxo3bswDDzzAW2+9RatWrcw+9rPPPiM2NhYfHx9Gjx7NvHnzuPvuu231bfDxxx+TmZlJhw4dcHd3Z8yYMaa/vFOmTGHAgAF06tSJbt26MXr0aIvbqV+/PmvXruXdd9/F3d2dzz//nHvvvZf69etXqB/Xnv/+++/j7u7OmjVrGDlypE2+x7Ief/xx2rdvT0BAAEOHDmXs2LE3fc7IkSM5evQorVq1uuEp8U2bNtG+fXtcXFyYPn06n332Waldsj59+hAUFMTgwYP5xz/+Qf/+/YEbv5+aNm1KdHQ0X375JV5eXrRt25YdO3YA8Ne//pWhQ4cyYMAAXFxc6N27N/v27bvVl8ZmDKqin7uFKKFr1648++yzPPbYY1XdlSp38uRJQkJCKrwLW9PJJw1RITExMVy8eJGCggI++OADjh8/zuDBg6u6W6IKyLA9USHHjh3j4YcfJjMzkzZt2vDll1+WOpMj6g7ZPRFCWEV2T4QQVqlWuyeenp4WrxO4FZmZmaUGT9UU0m/7kn7D77//Xm4AoyXVKjQCAwPZv3+/zbYXExNjuvCpJpF+25f0G8LCwir8WNk9EUJYRUJDCGEVCQ0hhFUkNIQQVpHQEEJYRUJDCGEVCQ0hhFUkNIQQVpHQEEJYRUJDCGEVCQ0hhFUkNIQQVrFJaEyYMAEvLy86duxoaktJSWHQoEGEhIQwaNCg6rEylBCi0mwSGo8//ni5Zf0WLFjAgAEDiI2NZcCAASxYsMAWpYQQVcwmodG3b1+aNWtWqm39+vWMHz8eMK43sW7dOluUEkJUMW3zaVy8eBFfX1/AuFiOpVWhIiMjiYyMBCAhIYGYmBib9SEjI8Om27MX6bd9Sb+tZKv1Hc+cOaNCQ0NN95s2bVrq625ubjfdhqU1LW9VycWWaxLpt31Jv6373dN29sTb25vExEQAEhMT8fLy0lVKCGFH2kJjxIgRpsWPo6KiuP/++3WVEkLYkU1CIyIigl69evHbb7/h7+/PBx98wN///neio6MJCQkhOjqav//977YoJYSoYjY5ELpy5Uqz7du2bbPF5oUQ1YiMCBVCWEVCQwhhFQkNIYRVJDSEEFaR0BBCWEVCQwhhFQkNIYRVJDSEEFaR0BBCWEVCQwhhFQkNIYRVJDSEEFaR0BBCWEVCQwhhFQkNIYRVJDSEEFaR0BBCWEVCQwhhFW3rnghR3V0mi03E4kou+RRSD8eq7lKNIJ80RJ30DbG0YjF/ZB0nSaEDb3OO9KruVo0goSHqnEKKmMxGssg3tZ0khZeIqbpO1SASGqLOSSCdONLKte8ivgp6U/NIaIg6x5smuNGgXPtteFZBb2oeCQ1R5zTAiVcIL9XmgjOz6Vs1Haph5OyJqJOepgd34ssXHMWPAg4zhQDcqrpbNYKEhqiz/kAr/kArYoiRwLCC7J4IIawioSFqnKvk8CuJ5JQ4ZapTMmlcIMUutWoC2T0RNcoitvMSm8kkDw8as4QHiOBOLbWyyOEN1vAzsQC0wY/neAhPmmqpV1PIJw1RY/zAGaazgUzyALhMJn/kE+K5oqXeJ2w1BQbAKc7zLhu01KpJJDREjbGBX8u1FVDEJo5pqbef38q1HeY0OcWhVVdJaIgaw5MmFtoba6nnQqNybQ2oT706vlcvoSFqjMcIKxcQITTnPkK11BtOr3Jt99Adxzr+a1O3I1PUKN648AN/YR5bOcIFehLADAbirOltHE4XGuDMVg6QTwF/oBOD6KqlVk0ioSFqlBCas4IIu9XrSQd60sFu9WqCuv05SwhhNQkNIYRVtO+eBAYG4uLigqOjI05OTuzfv193SVHN5FPALNaynO9QKMbzB+YxWtuxiPMkEcU6jnASbzx5mKF0p7OWWgBn+YVj7CSLdHxoQxeG0qgWDwCzyzGN7du34+kpcxXUVS+yhn/xjen+Iv5HEYp/azg2UUAhr/Iul4oHfMWTyL+JYi5/IZhWNq93kdPsZa3p/jmOc5UUBjMZAwab16sOZPdEaPcB31WozRYOc8IUGNcUUcQO9mmpd4afy7Wlk8QVzmupVx0YlFJKZ4HWrVvj7u6OwWBg0qRJTJw4sdTXIyMjiYyMBCAhIYFVq1bZrHZGRgZNmpgfEFSd1bZ+HySeQgpLtTngwB0a/vJnkcMFksu1u9IYT9zNPqcyr3cmV8gjx0w9Txypd0vbrChbvk+mT59e4UMH2kPj/Pnz+Pn5kZSUxKBBg1iyZAl9+5qfISksLMymxzxiYmIIDw+32fbspbb1exorWcyWUm1T6c9SHrN5H/Ip4CnmcqXEHKAGDPyTZ2hLoNnnVOb1TiSW7/m0VJsLHgxhqvbdE1u+T6z53dO+e+Ln5weAl5cXo0aNYu/evbpLimpmAaN5moE0oQGNcGYKd/M6Y7XUqocTM5lEKMEYMOCNJ08zzmJgVJYvIXRlOA1xBQx404Y+PFJrj2eA5gOhmZmZFBUV4eLiQmZmJlu2bGH27Nk6S4pqqD71eJNHeZNH7VKvFb68xFQUyi6/vEF0JYiudqtX1bSGxsWLFxk1ahQABQUFPPLIIwwdOlRnSSFM7P0LXBcCAzSHRlBQEIcOHdJZQghhZ3LtibAohau8znp+4Djt8ed5RhKEj7Z6JznNRraQRjp3cjtDGYCTpreooojT7OYch8jEiyRi8SJESy2jq8Au4ALgC/QGXDTW00dCQ5hVSCH9mc0hfgdgB0dYw25+YTHeGmbujuUUL7GAAgoA+JVjnOEsTzPxJs+8NYfYQCw7AHDChR28w108iQ/tNVTLA5YDqcX3zwLHgSmAs4Z6esngLmHW/zhoCoxrkkjjv3yrpd5GtpgC45rv2c1lDRP65pPDKX4o06r4je02r2V0hOuBcU0qcFRTPb0kNIRZFyzMu3mh3JvfNlLNrK2qUKRpWMk9n2yKzMxknqNt1fgMK9urNwkNYdYQ7sAJx3LtwzRNQnMnt5dr86AZgRpGjTbCnab4lWv31TQDGBaPleg8hqKPhIYwqwUefMBUXIvnyXTGiVmMYaCZX25buJdB9KaH6bRlM9x5lsk4aHqL9uBRGpdY8NmH9nRgsJZa4AMM5vohxHrAUMBbUz295ECosOiP3M2D9OJX4miDD564aqvlhBPPMplxjCGNdAJphaOZTzq20hQ/7uEfpHKenznOXQzUVsuoF3AHkAx4gplV62sKCQ1xQ41pQA/a2q2eJx544mGXWgYccMcfB07apZ4xKPztVEsf2T0RQlhFPmnUQhdJ4SOiucJVRnEXYbTTWu8Ev7GffTTGhXTScdW4G1NIPvHsJZ1zeNAGP+7AoPVv3xVgP1AAdME4MEujpDNw4RQ0coWA26Fefb31boGERi1zjLPcxTNcLj59OI9PeIu/8H+M1FLvK9bzEf8FYCjDeY5n+Cfz8dHwy1VIPjEsIIXTprYW3ElvnrZ5LaPfgXfAtKLaVuAxjOGhwS/R8FuJ8SMnfoT+T0ADPYtB3SrZPall5vKRKTCumckHZJNr81pZZPE5pSdNSiONdayxeS2AePaUCgyAc/zEJTPLJ9rGJii1BGMR8BWgYQqarDT47cfSbZmpcHKP7WtVkoRGLfMrZ8q1pZJBApdsXusSSeSambUqnjib1wJI45xV7ZWXaKYtBTQEMOnJmA2jtCTb16okCY1apoeZaye8cSdQw4VmvvjS2Mz6qsGaBi01I8hsu4eF9soLMNPmjZbTpW4+4GDmFHOzFravVUkSGrXMbMaXCoh6OLGEZ7QsWuxMff7En0sNwPLFl1E8aPNaYDx+4VvmeEIb+uOuaVYuGA6lQtEZNH1vNGgMnQaUbnP3heDueupVghwIrWX8ac4xoljH91zhKvfRG3+aa6vXl3Da04GfOEABRbzOG9TTNKGuA4704S8kcZQ0zuNBkMVPH7bhA8wEDmM8e9IJLKxcbxNte4FvCFw8DY2agk8IOFS/v+sSGrVQA5wZS3+71WuOF0O4hxhitAVGSV50wMtu66s2ALrZqRbg4mm8VWPVL8aEENWafNKoRlJIYwcHcNZxdN6MOE5znnjaEkoz9P51K6KAJA4DCi864aD9E0kaEAt4gYYrZctJiYOsVGjeBuprHleRVwCHT0N2HhQV2X0XRkKjmlhNNOOZTTY5vM5E/s0kvmIxjWho81pFFLKEeewqnlDHEUceYzL3MtrmtQCuco6dvExW8WnfhnjSl9m40lJLPfgW+BRMCzSFAZNBxwVwhXmwawVcPG6871APwh6GAD1TCHDyHLy2EtIyYWhbeO4dmDkOPO23dqzsnlQDmWTzBK+QXWLMw7fsZQm2W22upF18awoMME7t9yHvkIyeMQE/84EpMACySeZn3tdSyzgjVsnAAOMw8B/NP7yyTv1wPTAAivLhwOeQX378ik2895UxMK45lwyfbtNTywIJjWrgJ46RbmYWp+3YbrW5ko5Qfob4Ioo4ZqbdFi7xa7m2JI5oqWXcJSk0035MT7kkM1fIFuYZd1dsLTMHzl4s336k/IA+nSQ0qoFA/MxONtNG02XU3hauC/FGz0CixmYmm2mibQIaS6eXvfSUa2LuMn6DhfZKaugMTc0cL/FuZvtaNyChUQ20xIeJPFCqzR1XpmlakWwAw2le5pf2DnrQVtNpzI5EUPqtZiCUCC21IBC4s0ybOxCup1xIX3Au84sc2B0aawgNBwcY3a90m6MDjDa/NrIuciC0mnibGfSjK9+wCx88+JmVBJiZx9IWXGnKfN5lCxtIJIH2dCacIVpqAfjTm/548jvfAooA7saT27TVMy4N8APGZQK8gLsBTQcKG3vAoOeMxzayU8H7Nmh1h55aAEO6QQtP+P6w8Vua/yQE6luLxhwJjWrCgIGxDGUsQ4khRltgXOOKG6P5o9YaJXnQFg+7zQDmBPQtvtlBI3foNMw+tQA6tjbeYmLsHhgguydCCCtJaNhAFtkkmL2MWo90Usgws06IDooisrhIYal5JXTKBi7bqRaQkwaqyD61Cgoh+QoUmju7o0FaDlyw/doqsntSSXNZwmu8RwaZtCOIFSykV7kDcbaRTgpLmMUhfsCAAz0YwP8xhwbFywzYWiJ7+In/kMVFnHGlIxNow/1aahknuIkCNmOc+CYIeBZ0DQDLTIY9kXA5FpzCYf8KuPMxcND0K/H9TxC1DtKugrsrjB8JvTUd+8grhClfw0e/QH4R9GkFnzwArWxzXEc+aVTCOrYwi3+TgXGwzW+c5n4mkqNpGHgkr3KoeDlBRRG7ieYT3tRSK5c0fuQlsjCOC8gjnZ9YzGVtSwluATZwfaas08BraJklC2DPe8bAAGONMzvh+CY9tRIvwdJPjIEBcCUd3vwYLiTrqTfvO1h+0BgYAN/HwTjbzaYmoVEJqyn/JrvEZXZi+ynaCshnf/GCxSXtJtrmtQAS+ZFCM+GXYKYPtlF2bVWAcxgXS7ax7Ctw2cygrAQ9g+nY84vxGpGSiopg72E99b4wE+zfxcFF2+yqSGhUgquFuRVccbF5LQccqW/mOpSG6Lk4qp6F7dbTtCuExe3a/tobHOuDwcx1KPU0LWDUyMJ2G2maadzVzHbrOUBD21wkKKFRCZN4BGecS7WF0Zme2H5f1QEHhjCmXPtQxtq8FoAPPWhcZuSoE40I0Dae4x4oXpLxuq5oWbrQuREE9CrfHjygfJst9L4DXMv8gWnaBHppOqbxlJnZvh7rbD5MboEcCK2ELnQgmg+Zx9v8TgID6M3LTNNWbyxP4YIbO9mIE04M4EEGlhlJaiuOOBPOYo6wgmSO4EorOjCexhrmGjW6HZgBrMd40VkY8JCmWkDX8dDIw7hLktcAwiZDyx56ajVpBC8/Bav/B2cSoLU/PDQUGmv4FAXwSCdj/i7dB+m5MKYD/L2PzTYvoVFJfelBXzS92cpwwIH7+CP32WlQViO86Mbf7FLLKKz4ZgcOThA60niLidEXGNe08IZn7TeYjohOxpsGsnsihLCKhIYQwiraQ2Pz5s20a9eO4OBgFixYoLtcheSQwwvMJpDOhNKLt7RNCGMUx28s4M9MpDuvMp5T/KK13hnWs4WH+Jp7+Jl/UUCWxmoZwGKMxx/OgqZTwCaXj8OOl+CrCbBrAaQn6KulFGzcAtNmwNTnYeWXUFCgr96FdHj8U/B7CXouhq/Kz0NSHWg9plFYWMjUqVOJjo7G39+fbt26MWLECDp0sNdM0uZN4Tn+y6em+0/xAo44MpkJNq+VSTovE8FVrgCQyiVO8DOL2Yq7hjke4viG/cwx3T/JSrK5SG8W2byW0VwwjUvJBeZhPE1quwNvJplJsGMOFBaPHzm/B1JOwD1vg5OG06Ubt8CHJWZPW/s1ZGbBE4/ZvhbAsPfhp+IQTEyHkcvhu6ehd2s99W6R1k8ae/fuJTg4mKCgIJydnRk7dizr16/XWfKmMsjgE1aXa3+XFVrq7WaTKTCuyeYqu9igpd4pvijXdo7t5Gi5nuMCmB3Ipud7I27H9cC4JucKnN+rp150TPm2mO/1fNrYe/Z6YFxTpCBS0zSFlWBQSmkapwtffPEFmzdvZtmyZQB89NFH7Nmzh6VLl5oeExkZSWRkJAAJCQmsWmW7eTEzMjJo0qT0+fFCCjlI+ZF4DWhAqIY5HtK5TDLny7U3wwc3C7NMmet3RV3lDAVm1ldtSggONv9gmYdxZXWjjIxmNGmSgvGThoZrRnJSjKM5y2rsBc63PqDO4usdl2AmIAzQuhUYyo4pqaSMXPjNzHq77g0hyPyEPpV5n5Q1ffp09u+v2IhYrbsn5vLIUObFnjhxIhMnTgQgLCyM8PBwm9WPiYkxu71FvMs3bC3VNoe/E65hdqcULvA0/cgvcZWoI04s4n/4WVgdzFK/K+IE5zhU5noUD7rQnym3tL2bm4hxXk6IiYkgPHwlMA0tM2Wl/g7Rz2G8uK2YY33oFwn1XW95sxZf709Ww9dlJu3tEQYT7r7lWhYVFMLUVyGuTCiumwDh5k+dVuZ9Uhlad0/8/f2Jj4833U9ISMDPT+/kMhXxX97mPobigAMNachTPMkM/qqlVjN8mM57+BavN9ocf/7CmxYDo7JCiKAd43GkIWDAm570ROcB6FfANALWARgH3KenlFsg9HgWGhav0eLSAv4wo1KBcUMPjYSB/aCek3Gqve5dYdJ4PbWcHOHrJ6B78RotHo3h9RFwv56xFpWiNMrPz1etW7dWp0+fVrm5uapz587q119/tfj4rl272rT+9u3bb/j1DJWhclSOTWtaUqSKVKZKU0Wq6KaPvVm/K6JQ5al8lVnp7VRclk36XSFFhUrlZdhsczftd16eUjn2eZ8opZRKy1Yqv+CmD7Pl623N757W3RMnJyeWLl3KkCFDKCwsZMKECYSGhuosaZXGmi72MseAgUZo+otohgP17LCKWUmahkSbY3CAevb72VHPnq8j4Krpwjkb0T6M/N577+Xee+/VXUYIYSc16tqTZJJ5kzc5xCG60Y2neAo33LTVO8E+tvJfcsjiDzxAD1376oCikDN8xgV2kkVfrnIKF9poqwcXgZUYB2TdgXFwlsa/cBcPwO+bAQUBQ8BH40rseXkQ/T848iv4+MC9w8HT0nooNhB7AZZsg/gUGNIRnuhrPEZRS9WY0Mggg970Jrb4SP0GNvAFX7CPfdTT8DH8INtYwMMUFa/WtZt1RDCLUZoOmO7jBc5inF3JiVC2MpIBrMGVEA3VLmM8YHlt5qgY4HsgUkMt4Gw07J13/X7cNgh7Hlpr+gS6cAEc+vn6/e92wGuLoJmGtUhOXIDur0BatvH+up9gx2+wcrLta1UTNebak5WsNAXGNYc4xHr0DBZbyyJTYFyznjfI1zCVXyYJnGVtqbYCMojVNOAMvuR6YFyzD/jZzGNt4OiHZto+0lMr9kTpwABISzN+8tDhzejrgXHNqj3GTx+1VI0JjXjizbbHoWHAnG0yAAASDklEQVTNTCCZc+Xaskgnm6s2r5VNIubmwswyMyjMNsysBwoYR3hqkGVmYelsPYtNc9nCvJvJmubjjEsx3x5vob0WqDGhMZjBVrVX1u2UH8DTms644mnzWu50wtnMsRlv7rJ5LaOeZtqc0DaXhY+Z7XprqtUh1PzZjtu76Kk32MzZQNeG0EPn8aiqVWNCow99mMEMHDEeYKpHPRawgI501FJvLLMIpqvpvgd+TGHpDZ5x6xxpQHcWUa/E3KJ+DKYN47TUg4HAg1yfXq8BxlmzNB0s7PI0uARcv9+kJdzxjJ5ark1h4hSoXzy1ncEAdw+A3houoAOYFA6jSwSga0OIegIaa5r/sxqoMQdCAV7lVaYwhaMc5XZux1vbyuPgigfz2MopfiaXLNrRA0eNL5cv/RnObi7zE0fI5A88ra2WMSxeBB4H4oFQ0DmGpLEPDFkBl48ACjxCjWMtdOl3N4R1h5Ox4O0NPr43f86tqucEq6fC8UTjLkmvNtCkeo+zqKwaFRoA/sX/2UsbDZMEW+JEI7zpwzFi7FTRv/hmBwYDeOr5VGhW48b6dknMuc3XeKsDaszuiRCieqhxnzQq4iR7Ocw26tGOfHKph779y0JyOccGMonHm340K3EcRI/zwDcYf3TDQMOB2VKS9sPlg+AaBL599S1bCJCTDT/GQNoVuLMXtNI8+UzseVi3B1o7QmomuNlxaHoNVutCYxUvsoZXARjK68ygOy+zk0bYZh3LkvJJYzv3kla8VOGv/JP2TKcjM21ey2gn8CTXly5cBHyEcfp/Dfa/DGfWXb/veSf0fQccnS0/51alJMOsp+FS8WnfT9+HCc/A0FG2rwXw6U744xtQWASvD4fbnoKdc6FtCz31apFatXtyibOsZX6ptrP8whbe0VLvJB+YAuOa4/yHLAtjSirvn1Bq9farwL/0lLpytHRgACT/BPGb9dRbv/J6YFzzSSRkZdq+Vn4B/HWFMTCuuZgKL39u+1q1UK0KjbMcQlFUrv2MppGOqWYmCFYUksoRDdVyATPrj6Jp8tkrxy20H9NT70xs+bbcHEjUMHHwuRRjSJT10ynb16qFalVoBHA7BjPfUmtNZ0DcMDdBigNu6Lj8vz4QbKZd01QDbu3Mt7u311Mv0Mz3Vr8B+GrYXWjRDLzM7K7eoWdipNqmVoVGcwK4v8yKYK3oxCD0XDzUhidwpfQv0W08SyMd82MCxrEVJY8nuAAv6CnVLBQC7y/d5nkHtByqp97ICGheZtxNxBPQyDZzYJZSzwn+/SdwLPH292oKcx62fa1aqNYdCH2EeYQxgl/ZhjOtGcdenDVd8u1MUwbyLQlsIIs4vOiHBxov+aYf8C2wCeOPbjjaRnECdJsDAfdC8iFwbQ1+4frOnjRrDotWwA/bi8+e9IQAjUOxH+0H3YJh7R5o5QTHl4K7hoCqhWpdaAC0pSdt6UkMMdoC4xpHGhCgc6HiclpgPINiJ17djTd7aNAQ+ttxwqa2LeBvDxjXcpXAqLBatXsihNCvzoZGEQVcYidXOGCnipeA7UCifcqlHoWk3VCUd/PHVlZeLvy6z/h/e0i4BFsPQHKafeqJUupkaKRykE0EEUM/thLGNnqRq2UFsmveADoCDwCdMS5dqEn+VYiJgOhhsCMCvv4DXP5JX72j+2HKUHhlIiSchlenGkd26vL8uxA4FgZNB/8x8OaX+moJs+pkaOzjT2SXGICVwm5+1TaK8wgwh+uDsgqAhcBuPeWOLYVLJbadmwx7nzMuZmxrRYWw5EW4WmLMw6EfYEOU7WuB8dPF659dH5SVmw/PvgUny0+YJPSpc6GRyyVSOViu/QKapoNju5XtlXThu/JtGb9D5lnb10o4DZfNzAJ26Afb1wLYsq98m1IQXbHlBIVt1LnQcMIVJzNzR+gbW2FpcJKmaxwambk826E+1G9m+1puHuBo5gSch4/tawH4Wzi93NJLTz1hVp0LDUfq047nS7UZcKQ9MzRVHAZ0KNMWhPH4hga3TQZDmenuQsZDPQ2T7Lg2g0EPlm6rVx9GaFq68I9DIKDMALCubeEeO50SFkAtHadxMx14ERfaEc9nONGENkzCg16aqjkDGzEuD/AzxgOikwBN4wI8u8GAL+Hkx5CfBi2GQKuRemoBPP4CBHWAfTHg6gbzPoSAtnpquTWB3W/DG1/Ckd+hZ3t4ahQ41t41RqqjOhkaAC0ZQ0vG2KmaG9qGe5vj3gm6vWafWg4OED7CeIuJ0RcY1/g0g/l2HNwmyqlzuydCiMqR0LipPCDLfuUKMuxXq6gIsmy/jotFObnGm6jRJDQsygeewbhr4QKMxDiqU5OrsbD9HljbEtJ+heP/0VcL4NvPYUpveLwzPH8vnNC0uhpAZjY8Phtc+xhv42ZChh2DWNiUhIZF84ElQDZQBKwH/qSnlCqC7yMguXhQVlEBHH4F4tfe+Hm36rcD8O7f4ErxGIuzx2DBBMjRMEsWwPP/gaivjDNm5RfAJ5vg2YV6agntJDQs+tRM2ybAzIxPlXXlIGSYmTUq7gvb1wLYtaF8W0YqHNypp96n35hp0zRtoNBOQsMiczOYO6HlhJODhYl6HTTNou5koZ6zpnr1zdRroGFyYmEXEhoWmZvt61G0jK9w6wjNyk7eY4AgTYOk7h5TPjia+0NnTWvHTh5dvm3Sg+XbRI1QZ8dp3NwUjBeXvQ1kAg9jnA1ckz6fwi9z4MI2oBF0WwHe/fTUatkWZqyA1Yvh/Blo3x0e/Rs4mVk42RZmTwTnerBig/FakfHDYeYTemoJ7SQ0bujp4psd1PeEbsULTMfEgH+43nodextv9uDoaAwJCYpaQXZPhBBWkdAQQlhFW2jMmTOHFi1a0KVLF7p06cKmTZt0lapal76D7+6BzR3g56chV+MMYEVF8PVb8Ld+MGMgbPtQXy2A+PPw57/BbQPggclw8OjNnyNqPa3HNKZNm8b06dN1lqhaqYdg58Dr83BePQaX98DAvXrqffYqrP339fsnD0BeDtwz0fa18vIgPAJOxxnv/3Yatu6CXzdDK1nvtC6T3ZPKOP1e+Yl7r+yDyxqm8isqgs3vl2//5j3b1wLYuP16YFxzNQOiZE7Ous6glI7JI427J//9739xdXUlLCyMRYsW4e7uXu5xkZGRREZGApCQkMCqVats1oeMjAyaNNG4nkXmGchLKd/uEgJOtz7pjfl+Kzj9i/H/JTk6QUDHW65l0eVU+N3MQtbezcHf/Mxc2l9vTaTfMH36dPbvr9i0iZUKjYEDB3LhwoVy7a+++io9e/bE09MTg8HArFmzSExMZPny5TfcXlhYWIU7XhExMTGEh4fbbHvlnFsLP5SZgat+cxgWB463vkiTxX6//hjs/bp02+AJ8MSiW65l0aXLENAHsnNKt//wJfS60+xTtL/emki/rfvdq9Qxja1bt1bocU8++STDhw+vTKnqqcUoCH0ZfltovKS9SVvjoKxKBMYNPflvyM2CQ9+CwQF63AePztFTq7kHrH4LpswyHhB1bwr//KvFwBB1h7YDoYmJifj6Gie5Xbt2LR07avgIXR10mA1tp0NeMjRsCQaDvlpNm8PMLyE1ybhb4qJhsuCShvWHM/0gIdG4W9JA07UpokbRFhovvPACBw8exGAwEBgYyHvvaTpgVx04NQKnVvar52bH2bcdHSHA3371RLWnLTQ++ugjXZsWQlQhufbEHnKz4NsP4NQ+CLgdBjwJjTQsKXDNiZMQGQVXUuGB4TBsiL5aos6R0NCtsABe6Q+xe4obPoKdH8Kru8G5oe3r7f8Z+g2HrOLp9JZ/DK/MgFnP3/h5QlSQDO7S7cBXJQKj2Nlf4IfP9dSb9+/rgXHNa2/AVTtOICxqNQkN3S6YmcYP4KKF9so6daZ8W2YmXEjSU0/UORIauoWGW9deWeF9yre18oegQD31RJ0joaFbmzAY+Y/S4zeGPgUd++upN+t5uKPz9fuuLvDBElm6UNiMHAi1h0fmQf8/w5mfIKAz+LXTV8vTAw7EwI5dxrMnA/uBi4u+eqLOkdCwF582xps9GAzmd1OEsAHZPRFCWEU+aVQnp07Bli3g5wcFBeAkPx5R/cgnjeri3XehbVv4v/+DkychLAyuXKnqXglRjoRGdZCaCn/9q3F2rmsOHYLFi6uuT0JYIKFRHRw+DNnZ5dv37CnfJkQVk9CoDtq2NX/8IjTU/n0R4iYkNKoDb2/4299Kt/n7G3dZhKhm5PB8dTF3LgwcCN98Ay1bwi+/gJmJmIWoahIa1Ul4uPEWEyOBIaot2T0RQlhFQkMIYRUJDSGEVSQ0hBBWkdAQQlhFQkMIYRUJDSGEVSQ0hBBWkdAQQlhFQkMIYRUJDSGEVSQ0hBBWkdAQQlhFQkMIYRUJDSGEVSQ0hBBWkdAQQlhFQkMIYRUJDSGEVSQ0hBBWqVRorF69mtDQUBwcHNi/f3+pr82fP5/g4GDatWvH//73v0p1UghRfVRqNvKOHTuyZs0aJk2aVKr96NGjrFq1iiNHjnD+/HkGDhzIiRMncHR0rFRnhRBVr1KfNNq3b0+7du3Kta9fv56xY8dSv359WrduTXBwMHv37q1MKSFENaFl3ZNz587Rs2dP031/f3/OnTtn9rGRkZFERkYCkJCQQExMjM36kZGRYdPt2Yv0276k39a5aWgMHDiQCxculGt/9dVXuf/++80+RylVrs1gMJh97MSJE5k4cSIAYWFhhIeH36xLFRYTE2PT7dmL9Nu+pN/WuWlobN261eqN+vv7Ex8fb7qfkJCAn5+f1dsRQlQ/Wk65jhgxglWrVpGbm8uZM2eIjY2le/fuOkoJIeysUqGxdu1a/P39+fHHHxk2bBhDhgwBIDQ0lIceeogOHTowdOhQ3nrrLTlzIkQtUakDoaNGjWLUqFFmvzZz5kxmzpxZmc0LIaohGREqhLCKhIYQwioSGkIIq0hoCCGsIqEhhLCKhIYQwioSGkIIq0hoCCGsIqEhhLCKhIYQwioSGkIIqxiUuckvqoinpyeBgYE2296lS5do3ry5zbZnL9Jv+5J+w++//05ycnKFHlutQsPWwsLCyk14XBNIv+1L+m0d2T0RQlhFQkMIYRXHOXPmzKnqTujUtWvXqu7CLZF+25f0u+Jq9TENIYTtye6JEMIqEhpCCKvUutCoDevLzpkzhxYtWtClSxe6dOnCpk2bqrpLN7R582batWtHcHAwCxYsqOruVFhgYCCdOnWiS5cuhIWFVXV3bmjChAl4eXnRsWNHU1tKSgqDBg0iJCSEQYMGceXKFft0RtUyR48eVcePH1f9+vVT+/btM7UfOXJEde7cWeXk5KjTp0+roKAgVVBQUIU9teyll15SCxcurOpuVEhBQYEKCgpSp06dUrm5uapz587qyJEjVd2tCgkICFCXLl2q6m5UyI4dO9SBAwdUaGioqe35559X8+fPV0opNX/+fPXCCy/YpS+17pOGrC9rX3v37iU4OJigoCCcnZ0ZO3Ys69evr+pu1Tp9+/alWbNmpdrWr1/P+PHjARg/fjzr1q2zS19qXWhYcu7cOVq2bGm6f6P1ZauDpUuX0rlzZyZMmGC/j523oKa9riUZDAYGDx5M165dTesJ1yQXL17E19cXAF9fX5KSkuxSV8sC0LrpXl/WHm70PUyZMoVZs2ZhMBiYNWsWzz33HMuXL6+CXt5cdXtdrbFr1y78/PxISkpi0KBB3HbbbfTt27equ1Xt1cjQqA3ry1b0e3jyyScZPny45t7cuur2ulrjWj+9vLwYNWoUe/furVGh4e3tTWJiIr6+viQmJuLl5WWXunVm96QmrS+bmJho+vfatWtLHTGvbrp160ZsbCxnzpwhLy+PVatWMWLEiKru1k1lZmZy9epV07+3bNlSrV9nc0aMGEFUVBQAUVFRFj9l25xdDrfa0Zo1a1SLFi2Us7Oz8vLyUoMHDzZ9be7cuSooKEi1bdtWbdq0qQp7eWPjxo1THTt2VJ06dVL33XefOn/+fFV36YY2btyoQkJCVFBQkJo7d25Vd6dCTp06pTp37qw6d+6sOnToUO37PXbsWOXj46OcnJxUixYt1LJly1RycrLq37+/Cg4OVv3791eXL1+2S19kGLkQwip1ZvdECGEbEhpCCKtIaAghrCKhIYSwioSGEMIqEhpCCKtIaAghrPL/2QO7lJ2bW2MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot transformed grid points\n",
    "plt.figure(figsize=(4, 4), facecolor=\"w\")\n",
    "plt.scatter(uvgrid[0], uvgrid[1], s=36, c=colors, edgecolor=\"none\")\n",
    "plt.grid(True)\n",
    "plt.axis(\"equal\")\n",
    "plt.title(\"Transformed grid in u-v space\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rotation\n",
    "theta = np.pi/3 # 60 degree clockwise rotation\n",
    "B = np.column_stack([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]])\n",
    "print(B)\n",
    "\n",
    "uvgrid2 = np.dot(B, xygrid)\n",
    "\n",
    "plt.figure(figsize=(4, 4), facecolor=\"w\")\n",
    "plt.scatter(uvgrid2[0], uvgrid2[1], s=36, c=colors, edgecolor=\"none\")\n",
    "plt.grid(True)\n",
    "plt.axis(\"equal\")\n",
    "plt.title(\"Transformed grid in u-v space\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Linear transformation on images ##\n",
    "\n",
    "Here we will do basic liear transformations on images. There are many advanced python modules of manipulating images, and you may google to find out. One example is [here](https://automatetheboringstuff.com/chapter17/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import a image manipulating library \n",
    "#import imageio\n",
    "import matplotlib.image as mpimg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Images are Data###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load an image file.\n",
    "RGB_Image = mpimg.imread(\"./img/lennon.png\")\n",
    "\n",
    "# To see basic features of the imported image.\n",
    "print(type(RGB_Image), RGB_Image.dtype, RGB_Image.shape) \n",
    "print(' ' )\n",
    "print(RGB_Image) # See the image as ndarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, channel = (RGB_Image.shape)          # X size, Y size, number of channel \n",
    "                                           # (3 for RGB image).\n",
    "\n",
    "Red_Image = RGB_Image[:,:,0]               # Store Red channel in a Numpy array.\n",
    "Green_Image = RGB_Image[:,:,1]             # Store Green channel in a Numpy array.\n",
    "Blue_Image = RGB_Image[:,:,2]              # Store Blue channel in a Numpy array.\n",
    "\n",
    "#-------------------------:\n",
    "# Show images on screen   :\n",
    "#-------------------------:\n",
    "plt.figure(figsize = (7,7))\n",
    "plt.subplot(2,2,1); plt.imshow(RGB_Image); plt.title('RGB')\n",
    "plt.subplot(2,2,2); plt.imshow(Red_Image, cmap='gray'); plt.title('RED channel')\n",
    "plt.subplot(2,2,3); plt.imshow(Green_Image, cmap='gray'); plt.title('GREEN channel')\n",
    "plt.subplot(2,2,4); plt.imshow(Blue_Image, cmap='gray'); plt.title('BLUE channel')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Linear transformation##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a linear transformation function\n",
    "def linear_transformation(src, A):\n",
    "    m, n = src.shape\n",
    "    points = np.mgrid[0:n, 0:m].reshape((2, m*n))\n",
    "    new_points = np.linalg.inv(A).dot(points).round().astype(int)\n",
    "    x, y = new_points.reshape((2, m, n), order='F')\n",
    "    indices = x + n*y\n",
    "    return np.take(src, indices, mode='wrap')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 0], [0, 0.5]]) # Scaling\n",
    "dst = linear_transformation(Red_Image, A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize = (10,10))\n",
    "plt.subplot(1,2,1); plt.imshow(Red_Image, cmap='gray'); plt.title('Original')\n",
    "plt.subplot(1,2,2); plt.imshow(dst, cmap='gray'); plt.title('Transformed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = np.array([[1, 0], [-1, 1]]) # Rotation\n",
    "dst2 = linear_transformation(Red_Image, B)\n",
    "\n",
    "plt.figure(figsize = (10,10))\n",
    "plt.subplot(1,2,1); plt.imshow(Red_Image, cmap='gray'); plt.title('Original')\n",
    "plt.subplot(1,2,2); plt.imshow(dst2, cmap='gray'); plt.title('Transformed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C = np.array([[1, 0], [0, -1]]) # Reflection about x-axis\n",
    "#C = np.array([[-1, 0], [0, 1]]) # Reflection about y-axis\n",
    "#C = np.array([[0, 1], [1, 0]]) # Reflection about y-axis\n",
    "dst3 = linear_transformation(Red_Image, C)\n",
    "plt.figure(figsize = (10,10))\n",
    "plt.subplot(1,2,1); plt.imshow(Red_Image, cmap='gray'); plt.title('Original')\n",
    "plt.subplot(1,2,2); plt.imshow(dst3, cmap='gray'); plt.title('Transformed')\n",
    "plt.show"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises ##\n",
    "\n",
    "   1. Construct a $2\\times 2$ matrix $A_1$ that performs a horizontal shear transforming $\\vec{e}_2$ into $\\vec{e}_2+2\\vec{e}_1$ (leaving $\\vec{e}_1$ unchanged). Then, write a Python code to visualize the matrix transformation, as in Section 1.\n",
    "   2. Construct a $2\\times 2$ matrix $A_2$ that performs and reflects vectors through the line $y=-x$. Then, write a Python code to visualize the transformation, as in Section 1.\n",
    "   3. Write a Python code to compute $A =A_1A_2$ and visualize the matrix transformation $A\\vec{x}$, as in Section 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further Reading ##\n",
    "\n",
    "* [Affine transform of an image](https://matplotlib.org/gallery/images_contours_and_fields/affine_image.html#sphx-glr-gallery-images-contours-and-fields-affine-image-py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Appendix. Generating Animations </center> #\n",
    "\n",
    " (Note: This **only works** on `Linux` or `OS X` and *requires* `ImageMagick` to be available at the command line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate intermediate transforms ##\n",
    "\n",
    "To create the animated version, we need a series of intermediate grids that will \n",
    "smoothly transition from the original grid to the transformed grid. One way to achieve this is\n",
    "by constructing a series of $2\\times 2$ matrices that interpolate between the identity matrix \n",
    "$I = \\left[ \\begin{array}{cc}\n",
    "              1 & 0 \\\\ 0 & 1\n",
    "            \\end{array}\n",
    "      \\right]\n",
    "$\n",
    "and the target matrix \n",
    "$A = \\left[ \\begin{array}{cc}\n",
    "              2 & -1 \\\\ 1 & 1\n",
    "            \\end{array}\n",
    "      \\right].\n",
    "$ \n",
    "\n",
    "\n",
    "Supose we want to do this in $n$ steps. Then the $j^{th}$ matrix in this sequence is:\n",
    "$$\n",
    "  A_j = I+\\frac{j}{n}(A-I) = \\left[ \\begin{array}{cc}\n",
    "                      1+j/n & -j/n\\\\ j/n & 1\n",
    "                      \\end{array}\n",
    "                      \\right]\n",
    "$$\n",
    "where $j = 0 ... n$. The matrix product\n",
    "$$ A_j\\cdot \\mathrm{xygrid} $$\n",
    "computes grid coordinates for the $j^{th}$ intermediate transform.\n",
    "The following code block generates all the intermediate grids for a given target matrix, \n",
    "and returns the results in a 3D array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To animate the transform, we generate a series of intermediates\n",
    "# Function to compute all intermediate transforms\n",
    "def stepwise_transform(A, points, nsteps=30):\n",
    "    '''\n",
    "    Generate a series of intermediate transform for the matrix multiplication\n",
    "      np.dot(A, points) # matrix multiplication\n",
    "    starting with the identity matrix, where\n",
    "      a: 2-by-2 matrix\n",
    "      points: 2-by-n array of coordinates in x-y space \n",
    "\n",
    "    Returns a (nsteps + 1)-by-2-by-n array\n",
    "    '''\n",
    "    # create empty array of the right size\n",
    "    transgrid = np.zeros((nsteps+1,) + np.shape(points))\n",
    "    # compute intermediate transforms\n",
    "    for j in range(nsteps+1):\n",
    "        intermediate = np.eye(2) + j/nsteps*(A - np.eye(2)) \n",
    "        transgrid[j] = np.dot(intermediate, points) # apply intermediate matrix transformation\n",
    "    return transgrid\n",
    "\n",
    "# Apply to x-y grid\n",
    "steps = 30\n",
    "transform = stepwise_transform(B, xygrid, nsteps=steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot intermediate transforms ##\n",
    "\n",
    "Next we plot each of the intermediate grids on a common axis. To construct the animated version, \n",
    "we need to save each of these intermediate plots as an image file. The following code block defines \n",
    "a function that generates a series of image files with the filename `frame-xx.png` and saves them \n",
    "in a subdirectory `tmp`. We apply this function to the array of intermediate grid coordinates that we generated above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a series of figures showing the intermediate transforms\n",
    "def make_plots(transarray, color, outdir=\"png-frames\", figuresize=(4,4), figuredpi=150):\n",
    "    '''\n",
    "    Generate a series of png images showing a linear transformation stepwise\n",
    "    '''\n",
    "    nsteps = transarray.shape[0]\n",
    "    ndigits = len(str(nsteps)) # to determine filename padding\n",
    "    maxval = np.abs(transarray.max()) # to set axis limits\n",
    "    # create directory if necessary\n",
    "    import os\n",
    "    if not os.path.exists(outdir):\n",
    "        os.makedirs(outdir)\n",
    "    # create figure\n",
    "    plt.ioff()\n",
    "    fig = plt.figure(figsize=figuresize, facecolor=\"w\")\n",
    "    for j in range(nsteps): # plot individual frames\n",
    "        plt.cla()\n",
    "        plt.scatter(transarray[j,0], transarray[j,1], s=36, c=color, edgecolor=\"none\")\n",
    "        plt.xlim(1.1*np.array([-maxval, maxval]))\n",
    "        plt.ylim(1.1*np.array([-maxval, maxval]))\n",
    "        plt.grid(True)\n",
    "        plt.draw()\n",
    "        # save as png\n",
    "        outfile = os.path.join(outdir, \"frame-\" + str(j+1).zfill(ndigits) + \".png\")\n",
    "        fig.savefig(outfile, dpi=figuredpi)\n",
    "    plt.ion()\n",
    "\n",
    "# Generate figures\n",
    "make_plots(transform, colors, outdir=\"png-frames\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creat animation##\n",
    "\n",
    "To stitch the image sequence into an animation, we use the `ImageMagick`, a cross-platform image manipulation library. \n",
    "The following code block performs this operation by making a system call to the `convert` script that is part of ImageMagick. It will take some time to produce the animation gif file, `animation.gif`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to gif (works on linux/os-x, requires image-magick)\n",
    "from subprocess import call\n",
    "call(\"cd png-frames && convert -delay 10 frame-*.png animation.gif\", shell=True)\n",
    "# Optional: clean up png files\n",
    "call(\"rm -f png-frames/*.png\", shell=True)"
   ]
  },
  {
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}