{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Singular Value Decomposition of an image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise aims to show you how the Singular Value Decomposition can be used in order to store images saving space.\n",
    "\n",
    "The trivial way of storing an image of resolution, say $m \\times n$, would be to store the colour of each pixel. At least in the case of black-and-white images, we can assume that the colour can be stored as some number.\n",
    "Hence, the image would be stored as an $m \\times n$ matrix and take a lot of space.\n",
    "\n",
    "In class, you talked about approximating a matrix $A$ by $A_k$, where you would only use the first $k$ singular values for the reconstruction.\n",
    "We will have a look on the images induced by the approximated matrix $A_k$ for some values of $k$, and compare the space needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading in the Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start with the code. The first lines load some libraries we use, the most important one is numpy. It allows us to work with matrices, inlcuding performing the SVD in just one line later on.\n",
    "\n",
    "Afterwards, we load the image. For this example, I took an image of \"L'Oratoire Saint-Joseph\" in Montreal, Canada, but you should be able to use any picture you like. You can upload it and adjust the path in the code below.\n",
    "\n",
    "The next lines convert the picture to grayscales (for practical purposes) and output it in a certain size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time\n",
    "\n",
    "from PIL import Image\n",
    "\n",
    "img = Image.open('loratoire_st_joseph.jpg')  ## For your own image, adjust the path here.\n",
    "imggray = img.convert('LA')\n",
    "plt.figure(figsize=(10,7))  ## This allows you to reshape the image you output.\n",
    "plt.imshow(imggray);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we convert the image into a matrix and output the image again, in order to show that we did not alter the image itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imgmat = np.array(list(imggray.getdata(band=0)), float)\n",
    "imgmat.shape = (imggray.size[1], imggray.size[0])\n",
    "imgmat = np.matrix(imgmat)\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.imshow(imgmat, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the Singular Value Decomposition and Reconstruction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As promised, one line of command is enough to get the singular value decomposition $A = P D Q$, where $A$ is the initial matrix storing our image.\n",
    "In the code below, note that 'D' is only a vector containing the diagonal entries of the matrix $D$, as the other entries are zero. The other two lines reconstruct the matrix using the first singular value only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "P, D, Q = np.linalg.svd(imgmat)\n",
    "\n",
    "reconstimg = np.matrix(P[:, :1]) * np.diag(D[:1]) * np.matrix(Q[:1, :])\n",
    "plt.imshow(reconstimg, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see what we get when we use the second, third and fourth singular value as well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "reconstimg = np.matrix(P[:, :2]) * np.diag(D[:2]) * np.matrix(Q[:2, :])\n",
    "plt.imshow(reconstimg, cmap='gray')\n",
    "title = \"n = 2\"\n",
    "plt.title(title)\n",
    "plt.show()\n",
    "\n",
    "reconstimg = np.matrix(P[:, :3]) * np.diag(D[:3]) * np.matrix(Q[:3, :])\n",
    "plt.imshow(reconstimg, cmap='gray')\n",
    "title = \"n = 3\"\n",
    "plt.title(title)\n",
    "plt.show()\n",
    "\n",
    "reconstimg = np.matrix(P[:, :4]) * np.diag(D[:4]) * np.matrix(Q[:4, :])\n",
    "plt.imshow(reconstimg, cmap='gray')\n",
    "title = \"n = 4\"\n",
    "plt.title(title)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we cannot see much yet, but already for only four singular values the rather dark church and statue in front of the brighter sky are adumbrated.\n",
    "But now your part starts. The following code shows a for-loop that runs from 10 to 100 in steps of 10.\n",
    "Use the code above to complete the loop below.\n",
    "If you don't write anything, you will get an error when running the cell, as python does not expect an empty loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i in range(10, 51, 10): ## The syntax is range(start index, end index, step-width)\n",
    "    ## Use the code from above to fill the loop.\n",
    "    ## Note that in python, you indicate the part that belongs to the loop by alignment.\n",
    "    ## Thus, all your code should start where the hashes in this comment lines are.\n",
    "    ## (See Programming Exercise 1 for a few more things on python.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see how much space we save. The following command gives us the number of entries in D. As already mentioned, it is only a vector, as we only store the diagonal entries of $D$. We also output the size of the matrices $P$ and $Q$ (which you know already from the image)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"We have %d singular values.\" % D.shape)\n",
    "print(\"P is of size\", P.shape, \".\") \n",
    "print(\"Q is of size\", Q.shape, \".\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "As we store the whole singular value decomposition, we do not really save space in this exercise. But as you saw in the first theoretical exercise of the 10th series, we do not have to compute the whole matrices $P$ and $Q$ if we know that we only want to reconstruct the rank $k$ approximation. How many numbers do you have to store for the initial matrix of the picture? How many numbers do you have to store if you want to reconstruct the rank $k$ approximation only?\n",
    "\n",
    "Use the following Cell to find a large enough $i$ so that you are satisfied with the quality of the image.\n",
    "After that, finish the code to check how much percent of the initial size you have to store.\n",
    "\n",
    "If you are interested in keeping a fairly high quality, i.e. use a lot of the singular values, does it still make sense to use the singular value decomposition instead of the original matrix?\n",
    "At which point does the SVD start to use more space?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "i = 150 ## Number of used singular values\n",
    "reconstimg = np.matrix(P[:, :i]) * np.diag(D[:i]) * np.matrix(Q[:i, :])\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.imshow(reconstimg, cmap='gray')\n",
    "title = \"n = %s\" % i\n",
    "plt.title(title)\n",
    "plt.show()\n",
    "\n",
    "## This is the original image once again, to have direct comparison.\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.imshow(imgmat, cmap='gray');\n",
    "plt.title(\"Original Image\")\n",
    "plt.show()\n",
    "\n",
    "## The following gives you the resolution of your image, stored in a and b.\n",
    "## Use these numbers to calculate how much numbers your approximation needs.\n",
    "a,b = imgmat.shape\n",
    "\n",
    "numbers_apx = 1    ## How many numbers do you have to store for the approximated image?\n",
    "numbers_orig = 1   ## How many numbers do you have to store for the original image?\n",
    "\n",
    "print(\"For the original image, we have to store %d numbers.\" % numbers_orig)\n",
    "print(\"For this quality, we have to store %d numbers.\" % numbers_apx)\n",
    "print(\"Hence, %f percent.\" % (numbers_apx/numbers_orig))"
   ]
  },
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}