III_07_Image_compression_Change_of_basis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This notebook is part of lecture 31 *Change of basis and image compression* in the OCW MIT course 18.06 by Prof Gilbert Strang [1]\n",
    "+ Created by me, Dr Juan H Klopper\n",
    "    + Head of Acute Care Surgery\n",
    "    + Groote Schuur Hospital\n",
    "    + University Cape Town\n",
    "    + <a href=\"mailto:juan.klopper@uct.ac.za\">Email me with your thoughts, comments, suggestions and corrections</a> \n",
    "<a rel=\"license\" href=\"https://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br /><span xmlns:dct=\"https://purl.org/dc/terms/\" href=\"https://purl.org/dc/dcmitype/InteractiveResource\" property=\"dct:title\" rel=\"dct:type\">Linear Algebra OCW MIT18.06</span> <span xmlns:cc=\"https://creativecommons.org/ns#\" property=\"cc:attributionName\">IPython notebook [2] study notes by Dr Juan H Klopper</span> is licensed under a <a rel=\"license\" href=\"https://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>.\n",
    "\n",
    "+ [1] <a href=\"https://ocw.mit.edu/courses/mathematics/18-06sc-linear-algebra-fall-2011/index.htm\">OCW MIT 18.06</a>\n",
    "+ [2] Fernando Pérez, Brian E. Granger, IPython: A System for Interactive Scientific Computing, Computing in Science and Engineering, vol. 9, no. 3, pp. 21-29, May/June 2007, doi:10.1109/MCSE.2007.53. URL: https://ipython.org"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='https://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='https://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('https://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #ddd;\n",
       "    color: #000000;\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 2.2em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(0, 80, 120);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-weight: 400;\n",
       "    font-size: 1.9em;\n",
       "    line-height: 100%;\n",
       "    color: rgb(200,100,0);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: italic;\n",
       "    color: rgb(94,127,192);\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "        font-size: 100%;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML, Image\n",
    "css_file = 'style.css'\n",
    "HTML(open(css_file, 'r').read())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import init_printing, Matrix, symbols, sqrt, Rational\n",
    "from warnings import filterwarnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "init_printing(use_latex = 'mathjax')\n",
    "filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image compression and change of basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lossy image compression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Consider a 2<sup>9</sup> &times; 2<sup>9</sup> monochrome image\n",
    "+ Every pixel in this 512&times;512 image can take a value of 0 &le; *x*<sub>i</sub> < 255 (this is 8-bit)\n",
    "+ This make **x** a vector in &#8477;<sup>n</sup>, with *n* = 512<sup>2</sup> (for color images this would be 3&times;n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$262144$$"
      ],
      "text/plain": [
       "262144"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Just look at what 512 square is\n",
    "512 ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ This is a very large, unwieldy basis\n",
    "+ Consider the standard basis\n",
    "$$ \\begin{bmatrix} 1 \\\\ 0 \\\\ 0 \\\\ \\vdots  \\\\ 0 \\end{bmatrix},\\begin{bmatrix} 0 \\\\ 1 \\\\ 0 \\\\ \\vdots  \\\\ 0 \\end{bmatrix},\\cdots ,\\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ \\vdots  \\\\ 1 \\end{bmatrix} $$\n",
    "+ Consider now the better basis\n",
    "$$  \\begin{bmatrix} 1 \\\\ 1 \\\\ 1 \\\\ \\vdots  \\\\ 1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ \\vdots  \\\\ 1 \\\\ -1 \\\\ \\vdots  \\\\ -1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ -1 \\\\ 1 \\\\ -1 \\\\ \\vdots  \\end{bmatrix},\\cdots  $$\n",
    "+ Indeed, there are many options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ JPEG uses an 8 &times; 8 Fourier basis\n",
    "    + This means that an image is broken up into 64 &times; 64 pixel blocks\n",
    "    + See the lectures on the Fourier basis\n",
    "    $$ \\begin{bmatrix} 1 \\\\ 1 \\\\ 1 \\\\ \\vdots  \\\\ 1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ W \\\\ { W }^{ 2 } \\\\ \\vdots  \\\\ { W }^{ n-1 } \\end{bmatrix},\\cdots  $$\n",
    "    + This gives us a vector **x** in &#8477;<sup>64</sup> (i.e. with 64 coefficients)\n",
    "    + Up until this point the compression is lossless\n",
    "    + Now comes the compression (of which there are many such as thresholding)\n",
    "    + Thresholding\n",
    "        + Get rid of values more or less than set values (now there a less coefficients)\n",
    "        $$ \\hat{x}=\\sum{\\hat{c}_{i}{v}_{i}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Video is a sequence of images that are highly correlated (not big changes from one image to the next) and you can predict future changes from previous changes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ There are newer basis such as *wavelets*\n",
    "+ Here is an example\n",
    "$$  \\begin{bmatrix} 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ -1 \\\\ -1 \\\\ -1 \\\\ -1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ 1 \\\\ -1 \\\\ -1 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\end{bmatrix},\\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 1 \\\\ 1 \\\\ 1 \\\\ 1 \\end{bmatrix},\\begin{bmatrix} 1 \\\\ -1 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\end{bmatrix},\\begin{bmatrix} 0 \\\\ 0 \\\\ 1 \\\\ -1 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\end{bmatrix},\\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 1 \\\\ -1 \\\\ 0 \\\\ 0 \\end{bmatrix},\\begin{bmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 0 \\\\ 1 \\\\ -1 \\end{bmatrix} $$\n",
    "+ Every vector in &#8477;<sup>8</sup> is a linear combination of these 8 basis vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's do some linear algebra\n",
    "+ Consider only a top row of 8 pixels\n",
    "    + The standard vector of the values will be as follows (with 0 &le; *p*<sub>i</sub> < 255)\n",
    "    $$ \\begin{bmatrix} { p }_{ 1 } \\\\ { p }_{ 2 } \\\\ { p }_{ 3 } \\\\ { p }_{ 4 } \\\\ { p }_{ 5 } \\\\ { p }_{ 6 } \\\\ { p }_{ 7 } \\\\ { p }_{ 8 } \\end{bmatrix} $$\n",
    "+ We have to write this as a linear combination of the wavelet basis vectors *w*<sub>i</sub> (the lossless step)\n",
    "$$ {P}={c}_{1}{w}_{1}+{c}_{2}{w}_{2}+\\dots+{c}_{8}{w}_{8} $$\n",
    "+ In vector form we have the following\n",
    "$$ P=\\begin{bmatrix} \\vdots  & \\cdots  & \\vdots  \\\\ { w }_{ 1 } & \\cdots  & { w }_{ 8 } \\\\ \\vdots  & \\cdots  & \\vdots  \\end{bmatrix}\\begin{bmatrix} { c }_{ 1 } \\\\ \\vdots  \\\\ { c }_{ 8 } \\end{bmatrix} \\\\ P=Wc \\\\ c={W}^{-1}{P}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's bring some reality to this\n",
    "    + For fast computation, W must be as easy to invert as possible\n",
    "    + There is great competition to come up with *better* compression matrices\n",
    "    + A *good* matrix must have the following\n",
    "        + Be fast, i.e. the fast Fourier transform (FFT)\n",
    "            + The wavelet basis above is fast\n",
    "                + The basis vectors are orthogonal (and can be made orthonormal)\n",
    "                + **If they are orthonormal then the inverse is equal to the transpose**\n",
    "        + Good compression\n",
    "            + If we threw away some of the *p*<sub>i</sub> values, we would just have a dark image\n",
    "            + We we threw away, say the last two *c*<sub>i</sub> values (last two basis vectors) that won't lose us so much quality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Change of basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's look at this change in basis\n",
    "+ Above, we had the following\n",
    "$$ x=Wc $$\n",
    "+ Here W is the matrix that takes us from the vector **x** in the old basis to the vector **c** in the new basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Consider any transformation T (such as a rotation transformation)\n",
    "    + With respect to *v*<sub>1</sub>,...,*v*<sub>8</sub> it has a matrix A\n",
    "    + With respect to *w*<sub>1</sub>,...,*w*<sub>8</sub> it has a matrix B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Turns out that matrices A and B are similar\n",
    "$$ B={M}^{-1}AM $$\n",
    "+ Here M is the matrix that transforms the basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ What is A then, using the basis *v*<sub>1</sub>,...,*v*<sub>8</sub>?\n",
    "+ We know T completely from T(*v*<sub>i</sub>)...\n",
    "    + ... because if every **x**=&Sigma;*c*<sub>i</sub>*v*<sub>i</sub>\n",
    "    + ... then T(**x**)=&Sigma;*c*<sub>i</sub>T(*v*<sub>i</sub>)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Constructing A\n",
    "    + Write down all the transformations\n",
    "    $$ T\\left( { v }_{ 1 } \\right) ={ a }_{ 11 }{ v }_{ 1 }+{ a }_{ 21 }{ v }_{ 2 }+\\dots +{ a }_{ 81 }{ v }_{ 8 }\\\\ T\\left( { v }_{ 2 } \\right) ={ a }_{ 12 }{ v }_{ 1 }+{ a }_{ 22 }{ v }_{ 2 }+\\dots +{ a }_{ 82 }{ v }_{ 8 }\\\\ \\vdots \\\\ T\\left( { v }_{ 8 } \\right) ={ a }_{ 18 }{ v }_{ 1 }+{ a }_{ 28 }{ v }_{ 2 }+\\dots +{ a }_{ 88 }{ v }_{ 8 } $$\n",
    "    + Now we know A\n",
    "    $$ A=\\begin{bmatrix} { a }_{ 11 } & \\cdots  & { a }_{ 18 } \\\\ \\vdots  & \\cdots  & \\vdots  \\\\ { a }_{ 81 } & \\cdots  & { a }_{ 88 } \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Let's consider the linear transformation T(*v*<sub>i</sub>)=&lambda;<sub>i</sub>\n",
    "+ This makes A the following\n",
    "$$ A=\\begin{bmatrix} { \\lambda  }_{ 1 } & 0 & \\cdots  & \\cdots  & 0 \\\\ 0 & { \\lambda  }_{ 2 } & 0 & \\cdots  & \\vdots  \\\\ \\vdots  & 0 & \\ddots  & \\cdots  & \\vdots  \\\\ \\vdots  & \\vdots  & \\vdots  & \\ddots  & 0 \\\\ 0 & \\cdots  & \\cdots  & 0 & { \\lambda  }_{ 8 } \\end{bmatrix} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example problems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example problem 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ The vector space of all polynomials in *x* (of degree &le; 2) has the basis 1, *x*, *x*<sup>2</sup>\n",
    "+ Consider a different basis *w*<sub>1</sub>, *w*<sub>2</sub>, *w*<sub>3</sub> whose values at *x* = -1, 0, and 1 are given by the following\n",
    "$$ x=-1\\rightarrow 1{ w }_{ 1 }+{ 0w }_{ 2 }+{ 0w }_{ 3 }\\\\ x=0\\rightarrow 0{ w }_{ 1 }+1{ w }_{ 2 }+{ 0w }_{ 3 }\\\\ x=1\\rightarrow 0{ w }_{ 1 }+{ 0w }_{ 2 }+{ 1w }_{ 3 } $$\n",
    "+ Express *y*(*x*)=-*x*+5 in the new basis\n",
    "+ Find the change of basis matrices\n",
    "+ Find the matrix of taking derivatives in both of the basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ y\\left( x \\right) =5-x\\\\ y\\left( x \\right) =\\alpha { w }_{ 1 }+\\beta { w }_{ 2 }+\\gamma { w }_{ 3 } \\\\ y\\left( -1 \\right) =6 \\\\ y\\left( 0 \\right) =5\\\\ y\\left( 1 \\right) =4\\\\ \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 0 & 1 \\end{bmatrix}\\begin{bmatrix} \\alpha  \\\\ \\beta  \\\\ \\gamma  \\end{bmatrix}=\\begin{bmatrix} 6 \\\\ 5 \\\\ 4 \\end{bmatrix} \\\\ \\alpha =6,\\beta =5,\\gamma =4 \\\\ y=6{w}_{1}+5{w}_{2}+4{w}_{3} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ For the second part let's look at what happens at *x* for the various values at 1 (which is *x*<sup>0</sup>), *x*, and *x*<sup>2</sup>\n",
    "    + For -1 we have 1, -1, and 1\n",
    "    + For 0 we have 1, 0, and 0\n",
    "    + For 1 we have 1, 1, and 1\n",
    "+ From this we can conclude the following\n",
    "$$ 1={w}_{1}+{w}_{2}+{w}_{3} \\\\ x=-{w}_{1}+{w}_{3} \\\\ {x}^{2}={w}_{1}+{w}_{3} $$\n",
    "+ Now we have the following matrix\n",
    "$$ A=\\begin{bmatrix}1&-1&1\\\\1&0&0\\\\1&1&1\\end{bmatrix} $$\n",
    "+ This converts the first basis to the second\n",
    "+ To convert the second basis to the original we just need A<sup>-1</sup>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & 1 & 0\\\\- \\frac{1}{2} & 0 & \\frac{1}{2}\\\\\\frac{1}{2} & -1 & \\frac{1}{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡ 0    1    0 ⎤\n",
       "⎢             ⎥\n",
       "⎢-1/2  0   1/2⎥\n",
       "⎢             ⎥\n",
       "⎣1/2   -1  1/2⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix([[1, -1, 1], [1, 0, 0], [1, 1, 1]])\n",
    "A.inv()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Now for derivative matrices\n",
    "    + For the original basis, this is easy\n",
    "    $$ {D}_{x}=\\begin{bmatrix}0&1&0\\\\0&0&2\\\\0&0&0\\end{bmatrix} $$\n",
    "    + For the second basis we need the following\n",
    "    $$ {D}_{w}=AD{A}^{-1} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\left[\\begin{matrix}- \\frac{3}{2} & 2 & - \\frac{1}{2}\\\\- \\frac{1}{2} & 0 & \\frac{1}{2}\\\\\\frac{1}{2} & -2 & \\frac{3}{2}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡-3/2  2   -1/2⎤\n",
       "⎢              ⎥\n",
       "⎢-1/2  0   1/2 ⎥\n",
       "⎢              ⎥\n",
       "⎣1/2   -2  3/2 ⎦"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dx = Matrix([[0, 1, 0], [0, 0, 2], [0, 0, 0]])\n",
    "Dw = A * Dx * A.inv()\n",
    "Dw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "+ Just to conclude we can write the values for *w*<sub>i</sub> from the inverse of A (the columns)\n",
    "$$ {w}_{1}=\\frac{-1}{2}{x}+\\frac{1}{2}{x}^{2} \\\\ {w}_{2}=1-{x}^{2} \\\\ {w}_{3}=\\frac{1}{2}x+\\frac{1}{2}{x}^{2} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}