{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Audio File Matching\n",
    "\n",
    "This notebook continues the audio file matching problem. Be sure to have song.wav and clip.wav in the same directory as the notebook.\n",
    "\n",
    "In this notebook, we will look at the problem of searching for a small audio clip inside a song.\n",
    "\n",
    "The song \"Mandelbrot Set\" by Jonathan Coulton is licensed under <a href=\"http://creativecommons.org/licenses/by-nc/3.0/\">CC BY-NC 3.0</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the next block of code before preceding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import wave\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io.wavfile\n",
    "import operator\n",
    "from IPython.display import Audio\n",
    "%matplotlib inline\n",
    "\n",
    "given_file = 'song.wav'\n",
    "target_file = 'clip.wav'\n",
    "rate_given,  given_signal  = scipy.io.wavfile.read(given_file)\n",
    "rate_target, target_signal = scipy.io.wavfile.read(target_file)\n",
    "given_signal  = given_signal[:2000000].astype(float)\n",
    "target_signal = target_signal.astype(float)\n",
    "def play_clip(start, end, signal=given_signal):\n",
    "    scipy.io.wavfile.write('temp.wav', rate_given, signal[start:end].astype(np.int16))\n",
    "    return Audio(url='temp.wav', autoplay=True)\n",
    "\n",
    "def run_comparison(target_signal, given_signal, idxs=None):\n",
    "    # Run everything if not called with idxs set to something\n",
    "    if idxs is None:\n",
    "        idxs = [i for i in range(len(given_signal)-len(target_signal))]\n",
    "    return idxs, [vector_compare(target_signal, given_signal[i:i+len(target_signal)])\n",
    "                for i in idxs]\n",
    "\n",
    "play_clip(0, len(given_signal))\n",
    "\n",
    "#scipy.io.wavfile.write(target_file, rate_given, (-0.125*given_signal[1380000:1380000+70000]).astype(np.int16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will load the song into the variable `given_signal` and load the short clip into the variable `target_signal`. Your job is to finish code that will identify the short clip's location in the song. The clip we are trying to find will play after executing the following block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Audio(url=target_file, autoplay=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have defined the function 'vector_compare' to compare vectors by the \"cosine similarity measure.\"\n",
    "Run the following code that compares example vectors. Because the song has a lot of data, you should use the provided examples from the previous parts of the problem before running the later code. Do you results here make sense given your answers to previous parts of the problem?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def vector_compare(desired_vec, test_vec):\n",
    "    \"\"\"This function compares two vectors, returning a number.\n",
    "    The test vector with the largest magnitude return value is regarded as being closest to the desired vector.\"\"\"\n",
    "    return np.dot(desired_vec.T, test_vec)/(np.linalg.norm(desired_vec)*np.linalg.norm(test_vec)+1.e-10)\n",
    "\n",
    "print(\"PART A:\")\n",
    "print(vector_compare(np.array([1,1,1]), np.array([1,1,1])))\n",
    "print(vector_compare(np.array([1,1,1]), np.array([-1,-1,-1])))\n",
    "print(\"PART C:\")\n",
    "print(vector_compare(np.array([1,2,3]), np.array([1,2,3])))\n",
    "print(vector_compare(np.array([1,2,3]), np.array([2,3,4])))\n",
    "print(vector_compare(np.array([1,2,3]), np.array([3,4,5])))\n",
    "print(vector_compare(np.array([1,2,3]), np.array([4,5,6])))\n",
    "print(vector_compare(np.array([1,2,3]), np.array([5,6,7])))\n",
    "print(vector_compare(np.array([1,2,3]), np.array([6,7,8])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1\n",
    "Run the following code that runs `vector_compare` on every subsequence in the song- it will probably take at least 5 minutes. How do you interpret this plot to find where the clip is in the song?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "t0 = time.time()\n",
    "idxs, song_compare = run_comparison(target_signal, given_signal)\n",
    "t1 = time.time()\n",
    "plt.plot(idxs, song_compare)\n",
    "print (\"That took %(time).2f minutes to run\" % {'time':(t1-t0)/60.0} )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2\n",
    "Fill in the following code to use `song_compare` to print the index of `given_signal` where `target_signal` begins. Then, verify that your answer is correct by playing the song at that index using the `play_clip` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "index = np.????(np.abs(np.array(song_compare)))\n",
    "print (index)\n",
    "play_clip(index,index+len(target_signal))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hint: Have you heard of an argmax?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 7 Image Stitching\n",
    "\n",
    "This section of the notebook continues the image stiching problem. Be sure to have a `figures` folder in the same directory as the notebook. The `figures` folder should contain the files:\n",
    "\n",
    "    Berkeley_banner_1.jpg\n",
    "    Berkeley_banner_2.jpg\n",
    "    stacked_pieces.jpg\n",
    "    lefthalfpic.jpg\n",
    "    righthalfpic.jpg\n",
    "    \n",
    "Note: This structure is present in the provided HW2 zip file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the next block of code before proceeding\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.matlib\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from numpy import pi, cos, exp, sin\n",
    "import matplotlib.image as mpimg\n",
    "import matplotlib.transforms as mtransforms\n",
    "\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "#loading images\n",
    "image1=mpimg.imread('figures/Berkeley_banner_1.jpg')\n",
    "image1=image1/255.0\n",
    "image2=mpimg.imread('figures/Berkeley_banner_2.jpg')\n",
    "image2=image2/255.0\n",
    "image_stack=mpimg.imread('figures/stacked_pieces.jpg')\n",
    "image_stack=image_stack/255.0\n",
    "\n",
    "\n",
    "image1_marked=mpimg.imread('figures/lefthalfpic.jpg')\n",
    "image1_marked=image1_marked/255.0\n",
    "image2_marked=mpimg.imread('figures/righthalfpic.jpg')\n",
    "image2_marked=image2_marked/255.0\n",
    "\n",
    "def euclidean_transform_2to1(transform_mat,translation,image,position,LL,UL):\n",
    "    new_position=np.round(transform_mat.dot(position)+translation)\n",
    "    new_position=new_position.astype(int)\n",
    "\n",
    "    \n",
    "    if (new_position>=LL).all() and (new_position<UL).all():\n",
    "        values=image[new_position[0][0],new_position[1][0],:]\n",
    "    else:\n",
    "        values=np.array([2.0,2.0,2.0])\n",
    "    \n",
    "    return values\n",
    "\n",
    "def euclidean_transform_1to2(transform_mat,translation,image,position,LL,UL):\n",
    "    transform_mat_inv=np.linalg.inv(transform_mat)\n",
    "    new_position=np.round(transform_mat_inv.dot(position-translation))\n",
    "    new_position=new_position.astype(int)\n",
    "    \n",
    "    if (new_position>=LL).all() and (new_position<UL).all():\n",
    "        values=image[new_position[0][0],new_position[1][0],:]\n",
    "    else:\n",
    "        values=np.array([2.0,2.0,2.0])\n",
    "    \n",
    "    return values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will stick to a simple example and just consider stitching two images (if you can stitch two pictures, then you could conceivably stitch more by applying the same technique over and over again).\n",
    "\n",
    "Professor Alon decided to take an amazing picture of the Campanile overlooking the bay. Unfortunately, the field of view of his camera was not large enough to capture the entire scene, so he decided to take two pictues and stich them together. \n",
    "\n",
    "The next block will display the two images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "\n",
    "plt.figure(figsize=(20,40))\n",
    "\n",
    "plt.subplot(311)\n",
    "plt.imshow(image1)\n",
    "\n",
    "plt.subplot(312)\n",
    "plt.imshow(image2)\n",
    "\n",
    "plt.subplot(313)\n",
    "plt.imshow(image_stack)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you apply Marcela's algorithm on the two images you get the following result (run the next block):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20,30))\n",
    "\n",
    "plt.subplot(211)\n",
    "plt.imshow(image1_marked)\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.imshow(image2_marked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see Marcela's algorithm was able to find four common points between the two images. These points expressed in the coordinates of the first image and second image are\n",
    "\n",
    "$\\begin{align} \n",
    "&&\\vec{p_1}=\\begin{bmatrix} 200 \\\\ 700 \\end{bmatrix}\n",
    "&&\\vec{p_2}=\\begin{bmatrix} 310 \\\\ 620 \\end{bmatrix} \n",
    "&&\\vec{p_3}=\\begin{bmatrix} 390 \\\\ 660 \\end{bmatrix}\n",
    "&&\\vec{p_4}=\\begin{bmatrix} 460 \\\\ 630 \\end{bmatrix}\\\\\n",
    "&&\\vec{q_1}=\\begin{bmatrix} 162.2976 \\\\ 565.8862 \\end{bmatrix}\n",
    "&&\\vec{q_2}=\\begin{bmatrix} 285.4283 \\\\ 458.7469 \\end{bmatrix} \n",
    "&&\\vec{q_3}=\\begin{bmatrix} 385.2465 \\\\ 498.1973 \\end{bmatrix}\n",
    "&&\\vec{q_4}=\\begin{bmatrix} 465.7892 \\\\ 455.0132 \\end{bmatrix}\n",
    "\\end{align}$\n",
    "\n",
    "It should be noted that in relation to the image the positive x-axis is down and the positive y-axis is right. This will have no bearing as to how you solve the problem, however it helps in interpreting what the numbers mean relative to the image you are seeing.\n",
    "\n",
    "Using the points determine the parameters $R_{11},R_{12},R_{21},R_{22},T_x,T_y$ that map the points from the first image to the points in the second image by solving an appropriate system of equations. Hint: you do not need all the points to recover the parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Note that the following is a general template for solving for 6 unknowns from 6 equations represented as Az = b.\n",
    "# You do not have to use the following code exactly. \n",
    "# All you need to do is to find parameters R_11, R_12, R_21, R_22, T_x, T_y. \n",
    "# If you prefer finding them another way it is fine.\n",
    "\n",
    "# fill in the entries\n",
    "A = np.array([[?,?,?,?,?,?],\n",
    "              [?,?,?,?,?,?],\n",
    "              [?,?,?,?,?,?],\n",
    "              [?,?,?,?,?,?],\n",
    "              [?,?,?,?,?,?],\n",
    "              [?,?,?,?,?,?]])\n",
    "\n",
    "# fill in the entries\n",
    "b = np.array([[?],[?],[?],[?],[?],[?]])\n",
    "\n",
    "A = A.astype(float)\n",
    "b = b.astype(float)\n",
    "\n",
    "# solve the linear system for the coefficiens\n",
    "z = np.linalg.solve(A,b)\n",
    "\n",
    "#Parameters for our transformation\n",
    "R_11 = z[0,0]\n",
    "R_12 = z[1,0]\n",
    "R_21 = z[2,0]\n",
    "R_22 = z[3,0]\n",
    "T_x  = z[4,0]\n",
    "T_y  = z[5,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stitch the images using the transformation you found by running the code below.\n",
    "\n",
    "### Note that it takes about 40 seconds for the block to finish running on a modern laptop.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "matrix_transform=np.array([[R_11,R_12],[R_21,R_22]])\n",
    "translation=np.array([T_x,T_y])\n",
    "\n",
    "#Creating image canvas (the image will be constructed on this)\n",
    "num_row,num_col,blah=image1.shape\n",
    "image_rec=1.0*np.ones((int(num_row),int(num_col),3))\n",
    "\n",
    "#Reconstructing the original image\n",
    "\n",
    "LL=np.array([[0],[0]]) #lower limit on image domain\n",
    "UL=np.array([[num_row],[num_col]]) #upper limit on image domain\n",
    "\n",
    "for row in range(0,int(num_row)):\n",
    "    for col in range(0,int(num_col)):\n",
    "        #notice that the position is in terms of x and y, so the c  \n",
    "        position=np.array([[row],[col]])       \n",
    "        if image1[row,col,0] > 0.995 and image1[row,col,1] > 0.995 and image1[row,col,2] > 0.995:\n",
    "            temp = euclidean_transform_2to1(matrix_transform,translation,image2,position,LL,UL)\n",
    "            image_rec[row,col,:] = temp\n",
    "        else:\n",
    "            image_rec[row,col,:] = image1[row,col,:]\n",
    "        \n",
    "\n",
    "plt.figure(figsize=(20,20))\n",
    "plt.imshow(image_rec)\n",
    "plt.axis('on')\n",
    "plt.show()"
   ]
  }
 ],
 "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"
  },
  "name": "graphs_for_SOE.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 0
}