02_SVD_Eigenfaces_empty.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import scipy.io\n",
    "from matplotlib.image import imread\n",
    "import skimage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compression with SVD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [16, 8]\n",
    "\n",
    "\n",
    "A = imread(os.path.join('.','DATA','dog.jpg'))\n",
    "X = np.mean(A, -1); # Convert RGB to grayscale\n",
    "\n",
    "img = plt.imshow(X)\n",
    "img.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using [`np.linalg.svd`](https://docs.scipy.org/doc/numpy-1.10.4/reference/generated/numpy.linalg.svd.html) compress the above picture with the $r$ first singular values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code\n",
    "\n",
    "j = 0\n",
    "for r in (5, 20, 100):\n",
    "    # Construct approximate image\n",
    "    Xapprox = # your code\n",
    "    plt.figure(j+1)\n",
    "    j += 1\n",
    "    img = plt.imshow(Xapprox)\n",
    "    img.set_cmap('gray')\n",
    "    plt.axis('off')\n",
    "    plt.title('r = ' + str(r))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [Eigenface](https://en.wikipedia.org/wiki/Eigenface)\n",
    "\n",
    "dataset: 38 individuals under 9 poses and 64 lighting conditions. Each image is 192x168 = 32256 elements.\n",
    "\n",
    "We use the first 36 people as our training data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [10, 10]\n",
    "#plt.rcParams.update({'font.size': 18})\n",
    "\n",
    "mat_contents = scipy.io.loadmat(os.path.join('.','DATA','allFaces.mat'))\n",
    "faces = mat_contents['faces']\n",
    "m = int(mat_contents['m'])\n",
    "n = int(mat_contents['n'])\n",
    "nfaces = np.ndarray.flatten(mat_contents['nfaces'])\n",
    "\n",
    "allPersons = np.zeros((n*6,m*6))\n",
    "count = 0\n",
    "\n",
    "for j in range(6):\n",
    "    for k in range(6):\n",
    "        allPersons[j*n : (j+1)*n, k*m : (k+1)*m] = np.reshape(faces[:,np.sum(nfaces[:count])],(m,n)).T\n",
    "        count += 1\n",
    "        \n",
    "img = plt.imshow(allPersons)\n",
    "img.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for person in range(len(nfaces)):\n",
    "    subset = faces[:,sum(nfaces[:person]) : sum(nfaces[:(person+1)])]\n",
    "    allFaces = np.zeros((n*8,m*8))\n",
    "    \n",
    "    count = 0\n",
    "    \n",
    "    for j in range(8):\n",
    "        for k in range(8):\n",
    "            if count < nfaces[person]:\n",
    "                allFaces[j*n:(j+1)*n,k*m:(k+1)*m] = np.reshape(subset[:,count],(m,n)).T\n",
    "                count += 1\n",
    "                \n",
    "    img = plt.imshow(allFaces)\n",
    "    img.set_cmap('gray')\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We use the first 36 people for training data\n",
    "trainingFaces = faces[:,:np.sum(nfaces[:36])]\n",
    "avgFace = np.mean(trainingFaces,axis=1) # size n*m by 1\n",
    "\n",
    "# Compute eigenfaces on mean-subtracted training data\n",
    "X = trainingFaces - np.tile(avgFace,(trainingFaces.shape[1],1)).T\n",
    "U, S, VT = np.linalg.svd(X,full_matrices=0)\n",
    "\n",
    "fig1 = plt.figure()\n",
    "ax1 = fig1.add_subplot(121)\n",
    "img_avg = ax1.imshow(np.reshape(avgFace,(m,n)).T)\n",
    "img_avg.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "\n",
    "ax2 = fig1.add_subplot(122)\n",
    "img_u1 = ax2.imshow(np.reshape(U[:,0],(m,n)).T)\n",
    "img_u1.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the eigenface recosntruction for various order r."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Now show eigenface reconstruction of image that was omitted from test set\n",
    "\n",
    "testFace = faces[:,np.sum(nfaces[:36])] # First face of person 37\n",
    "plt.imshow(np.reshape(testFace,(m,n)).T)\n",
    "plt.set_cmap('gray')\n",
    "plt.title('Original Image')\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "r_list = [25, 50, 100, 200, 400, 800, 1600]\n",
    "\n",
    "for r in r_list:\n",
    "    reconFace = #your code here\n",
    "    img = plt.imshow(np.reshape(reconFace,(m,n)).T)\n",
    "    img.set_cmap('gray')\n",
    "    plt.title('r = ' + str(r))\n",
    "    plt.axis('off')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importance of data alignment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [8,8]\n",
    "plt.rcParams.update({'font.size': 18})\n",
    "\n",
    "n = 1000\n",
    "q = int(n/4)\n",
    "X = np.zeros((n,n))\n",
    "X[(q-1):(3*q),(q-1):(3*q)] = 1\n",
    "\n",
    "plt.imshow(X)\n",
    "plt.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "plt.title('Unrotated Matrix')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_rot = skimage.transform.rotate(X,10)\n",
    "X_rot[np.nonzero(X_rot)] = 1\n",
    "\n",
    "\n",
    "plt.imshow(X_rot)\n",
    "plt.set_cmap('gray')\n",
    "plt.axis('off')\n",
    "plt.title('Rotated Matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, VT = np.linalg.svd(X,full_matrices=0)\n",
    "\n",
    "\n",
    "plt.semilogy(S,'-o',color='k')\n",
    "plt.ylim((10**(-16),10**(4)+1))\n",
    "plt.yticks(np.power(10,np.arange(-16,4,4, dtype=float)))\n",
    "plt.xticks(np.arange(0,1000,250))\n",
    "plt.grid()\n",
    "plt.title('Unrotated Matrix: Spectrum')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U_rot, S_rot, VT_rot = np.linalg.svd(X_rot,full_matrices=0)\n",
    "\n",
    "plt.semilogy(S_rot,'-o',color='k')\n",
    "plt.ylim((10**(-16),10**(4)+1))\n",
    "plt.yticks(np.power(10,np.arange(-16,4,4, dtype=float)))\n",
    "plt.xticks(np.arange(0,1000,250))\n",
    "plt.grid()\n",
    "plt.title('Rotated Matrix: Spectrum')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [12,6]\n",
    "plt.rcParams.update({'font.size': 18})\n",
    "\n",
    "n = 1000\n",
    "q = int(n/4)\n",
    "X = np.zeros((n,n))\n",
    "X[(q-1):(3*q),(q-1):(3*q)] = 1\n",
    "\n",
    "nAngles = 12 # Sweep through 12 different angles, from 0:4:44 degrees\n",
    "cm_np = np.array([[0,0,2/3],\n",
    "                 [0,0,1],\n",
    "                 [0,1/3,1],\n",
    "                 [0,2/3,1],\n",
    "                 [0,1,1],\n",
    "                 [1/3,1,2/3],\n",
    "                 [2/3,1,1/3],\n",
    "                 [1,1,0],\n",
    "                 [1,2/3,0],\n",
    "                 [1,1/3,0],\n",
    "                 [1,0,0],\n",
    "                 [2/3,0,0]])\n",
    "\n",
    "\n",
    "cmap = plt.cm.jet\n",
    "cmap.set_bad(alpha = 0.0)\n",
    "\n",
    "U, S, VT = np.linalg.svd(X,full_matrices=0)\n",
    "\n",
    "fig1 = plt.figure()\n",
    "ax1 = fig1.add_subplot(121)\n",
    "plt.imshow(np.ma.masked_where(X == 0, X),vmin=0,vmax=nAngles)\n",
    "plt.set_cmap(cmap)\n",
    "plt.axis('off')\n",
    "\n",
    "\n",
    "ax2 = fig1.add_subplot(122)\n",
    "plt.semilogy(S,'-o',color=tuple(cm_np[0]))\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xrot = X\n",
    "\n",
    "fig, axs = plt.subplots(1,2)\n",
    "\n",
    "for j in range(nAngles):\n",
    "    Xrot = skimage.transform.rotate(X,j*4) #rotate by theta = j*4 degrees\n",
    "    Xrot[np.nonzero(Xrot)] = j\n",
    "    \n",
    "    U, S, VT = np.linalg.svd(Xrot)\n",
    "\n",
    "    axs[0].imshow(np.ma.masked_where(Xrot == 0, Xrot),vmin=0,vmax=nAngles)\n",
    "    plt.set_cmap(cmap)\n",
    "    axs[0].axis('off')\n",
    "\n",
    "    axs[1].semilogy(S,'-o',color=tuple(cm_np[j]))\n",
    "    axs[1].axis('on')\n",
    "    axs[1].grid(1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ml-icfp",
   "language": "python",
   "name": "ml-icfp"
  },
  "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}