{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###Problem 8.1 (b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Load A matrix and b vector\n",
    "phantom1_data = io.loadmat('phantom1.mat')\n",
    "A = double(phantom1_data['A'])\n",
    "b = phantom1_data['b']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Find the min-norm solution\n",
    "#BEGIN SOLUTION HERE\n",
    "#image_guess = \n",
    "#END SOLUTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Plot the recovered image\n",
    "image_reconstructed = np.reshape(image_guess, (50,50)).T\n",
    "plt.imshow(image_reconstructed, cmap='gray', interpolation='nearest', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Load original image to compare\n",
    "image = double(phantom1_data['image'])\n",
    "plt.imshow(image, cmap='gray', interpolation='nearest', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Repeat with a different matrix and image\n",
    "#Load A matrix and b vector\n",
    "phantom2_data = io.loadmat('phantom2.mat')\n",
    "A = double(phantom2_data['A'])\n",
    "b = phantom2_data['b']\n",
    "\n",
    "#Find the min-norm solution\n",
    "#BEGIN SOLUTION HERE\n",
    "image_guess = \n",
    "#END SOLUTION\n",
    "\n",
    "#Plot the recovered image\n",
    "image_reconstructed = np.reshape(image_guess, (50,50)).T\n",
    "plt.imshow(image_reconstructed, cmap='gray', interpolation='nearest', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Load original image to compare\n",
    "image = double(phantom2_data['image'])\n",
    "figure; plt.imshow(image, cmap='gray', interpolation='nearest', vmin=0, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Better Recovery?\n",
    "Now, notice how the resulting images are fairly noisy. However, there are more advanced methods that help get more accurate data even after using significantly less data. \n",
    "\n",
    "Run cell below to see the recovered image changes as a smaller percentage of the data is used for recovery. Why does this work significantly better than the min-norm solution?\n",
    "\n",
    "In this case, the image is converted to a form of the frequency domain, and the most important frequency components are selected. Therefore, the sampling is done in the frequency domain, not the spatial domain. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = io.loadmat('wavelet.mat')\n",
    "image = data['im']\n",
    "figure(figsize=(10,8)); subplot(2,3,1); title('Original Image')\n",
    "plt.imshow(image, cmap='gray', interpolation='nearest', vmin=0, vmax=1)\n",
    "\n",
    "#0.2% of Data\n",
    "recon = abs(data['im_rec_002'])\n",
    "subplot(2,3,6); title('0.2% of Data')\n",
    "plt.imshow(recon, cmap='gray', interpolation='nearest', vmin=0, vmax=1)\n",
    "\n",
    "#1% of Data\n",
    "recon = abs(data['im_rec_01'])\n",
    "subplot(2,3,5); title('1% of Data')\n",
    "plt.imshow(recon, cmap='gray', interpolation='nearest', vmin=0, vmax=1)\n",
    "\n",
    "#10% of Data\n",
    "recon = abs(data['im_rec_1'])\n",
    "subplot(2,3,4); title('10% of Data')\n",
    "plt.imshow(recon, cmap='gray', interpolation='nearest', vmin=0, vmax=1)\n",
    "\n",
    "#25% of Data\n",
    "recon = abs(data['im_rec_25'])\n",
    "subplot(2,3,3); title('25% of Data')\n",
    "plt.imshow(recon, cmap='gray', interpolation='nearest', vmin=0, vmax=1)\n",
    "\n",
    "#50% of Data\n",
    "recon = abs(data['im_rec_50'])\n",
    "subplot(2,3,2); title('50% of Data')\n",
    "plt.imshow(recon, cmap='gray', interpolation='nearest', vmin=0, vmax=1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}