Light_Field_Blind_Motion_Deblurring.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import tensorflow as tf\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from scipy import io\n",
    "from scipy import interpolate\n",
    "from scipy import ndimage\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameters: see parameters_to_use.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#camera parameters\n",
    "cf = 3.19 #sensor crop factor\n",
    "ml_pitch = 0.02 #distance between microlenses: 20 um\n",
    "num_sa = 14.0 #number of subaperture images\n",
    "num_sa_d = 8.0 #number of desired subaperture images (cropped to only consider those that lie within camera aperture)\n",
    "lfsize = [8, 8, 200, 200] #size of 4D light field\n",
    "f35 = 30.0 #focal length (35 mm equivalent)\n",
    "f_num = 2.0 #f-number\n",
    "f = f35/cf\n",
    "aperture_size = f/f_num\n",
    "aperture_res = (aperture_size/num_sa)/ml_pitch\n",
    "\n",
    "#forward model parameters\n",
    "path_order = 3 #implemented for either 2 (quadratic) or 3 (cubic) Bezier curves\n",
    "num_exp_pts = 7 #number of exposure points to average along motion curve (try {5, 7, 10})\n",
    "\n",
    "#blind deblurring parameters\n",
    "lam_init = 8e-3 #initial regularization weight\n",
    "num_iters = 1000 #optimization iterations\n",
    "eta_pts = 1.0 #step size for control points\n",
    "eta_lf = 0.4 #step size for sharp light field\n",
    "eps_init = 0.12 #initial parameter for evolving L0 approximation (try 0.1 or 0.15 for real data, 0.12 for synthetic data)\n",
    "lam_decay = 0.997 #regularization weight decay\n",
    "eps_decay = 0.997 #evolving L0 approximation decay\n",
    "lam_min = 0.1*lam_init #minimum regularization weight (try 0.25*lam_init for real data, 0.1*lam_init for synthetic data)\n",
    "eps_min = 0.04 #minimum L0 approximation parameter \n",
    "\n",
    "#final light field restoration parameters\n",
    "lam_lf = 0.004 #TV regularization weight (try {0.002, 0.04})\n",
    "iters_lf = 100 #TV regularization iterations\n",
    "eta_lf = 0.04 #TV regularization step size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize_lf(lf):\n",
    "    #normalize to between 0 and 1\n",
    "    return ((lf-np.amin(lf))/(np.amax(lf)-np.amin(lf)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def meshgrid_4D(v_vals, u_vals, y_vals, x_vals):\n",
    "    #4D meshgrid\n",
    "    with tf.name_scope('meshgrid_4D'):\n",
    "        \n",
    "        size = [tf.size(v_vals), tf.size(u_vals), tf.size(y_vals), tf.size(x_vals)]\n",
    "        v_vals_e = tf.expand_dims(tf.expand_dims(tf.expand_dims(v_vals, 1), 2), 3)\n",
    "        u_vals_e = tf.expand_dims(tf.expand_dims(tf.expand_dims(u_vals, 1), 2), 3)\n",
    "        y_vals_e = tf.expand_dims(tf.expand_dims(tf.expand_dims(y_vals, 1), 2), 3)\n",
    "        x_vals_e = tf.expand_dims(tf.expand_dims(tf.expand_dims(x_vals, 1), 2), 3)\n",
    "        v = tf.tile(tf.reshape(v_vals_e, [-1, 1, 1, 1]), [1, size[1], size[2], size[3]])\n",
    "        u = tf.tile(tf.reshape(u_vals_e, [1, -1, 1, 1]), [size[0], 1, size[2], size[3]])\n",
    "        y = tf.tile(tf.reshape(y_vals_e, [1, 1, -1, 1]), [size[0], size[1], 1, size[3]])\n",
    "        x = tf.tile(tf.reshape(x_vals_e, [1, 1, 1, -1]), [size[0], size[1], size[2], 1])\n",
    "        \n",
    "        return v, u, y, x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gather_4D(params, indices):\n",
    "    #gather values from tensor params at indices\n",
    "    with tf.name_scope('gather_4D'):\n",
    "        \n",
    "        flat_params = tf.reshape(params, [-1])\n",
    "        multipliers = [lfsize[3]*lfsize[2]*lfsize[1], lfsize[3]*lfsize[2], lfsize[3], 1]\n",
    "        ind_0 = tf.multiply(tf.slice(indices, [0,0], [-1,1]), multipliers[0])\n",
    "        ind_1 = tf.multiply(tf.slice(indices, [0,1], [-1,1]), multipliers[1])\n",
    "        ind_2 = tf.multiply(tf.slice(indices, [0,2], [-1,1]), multipliers[2])\n",
    "        ind_3 = tf.multiply(tf.slice(indices, [0,3], [-1,1]), multipliers[3])\n",
    "        flat_indices = tf.add_n([ind_0, ind_1, ind_2, ind_3])\n",
    "                \n",
    "        return tf.gather(flat_params, flat_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lf_reparam_tf(lf, coord):\n",
    "    #tensorflow forward model to reparameterize light field by moving the two parameterization planes by coord\n",
    "    with tf.name_scope('lf_reparam_tf'):\n",
    "        \n",
    "        z_multiple = 20.0 #scale z coordinates to similar magnitude as x/y coordinates for optimization\n",
    "                \n",
    "        #create and reparameterize light field grid\n",
    "        v_vals = tf.multiply(aperture_res, tf.subtract(tf.cast(tf.range(lfsize[0]), tf.float32), \n",
    "                                             tf.divide(tf.subtract(tf.cast(lfsize[0], tf.float32), 1.0), 2.0)))\n",
    "        u_vals = tf.multiply(aperture_res, tf.subtract(tf.cast(tf.range(lfsize[1]), tf.float32), \n",
    "                                             tf.divide(tf.subtract(tf.cast(lfsize[1], tf.float32), 1.0), 2.0)))\n",
    "        y_vals = tf.subtract(tf.cast(tf.range(lfsize[2]), tf.float32), \n",
    "                                tf.divide(tf.subtract(tf.cast(lfsize[2], tf.float32), 1.0), 2.0))\n",
    "        x_vals = tf.subtract(tf.cast(tf.range(lfsize[3]), tf.float32), \n",
    "                                tf.divide(tf.subtract(tf.cast(lfsize[3], tf.float32), 1.0), 2.0))\n",
    "    \n",
    "        v, u, y, x = meshgrid_4D(v_vals, u_vals, y_vals, x_vals)\n",
    "        \n",
    "        v_r = tf.add(tf.subtract(v, tf.multiply(coord[2]/z_multiple, tf.subtract(v, y))), coord[1])\n",
    "        u_r = tf.add(tf.subtract(u, tf.multiply(coord[2]/z_multiple, tf.subtract(u, x))), coord[0])\n",
    "        y_r = tf.add(tf.subtract(y, tf.multiply(coord[2]/z_multiple, tf.subtract(v, y))), coord[1])\n",
    "        x_r = tf.add(tf.subtract(x, tf.multiply(coord[2]/z_multiple, tf.subtract(u, x))), coord[0])\n",
    "        \n",
    "        v_r = tf.add(tf.div(v_r, aperture_res), tf.div(tf.subtract(tf.to_float(lfsize[0]), 1.0), 2.0))\n",
    "        u_r = tf.add(tf.div(u_r, aperture_res), tf.div(tf.subtract(tf.to_float(lfsize[1]), 1.0), 2.0))\n",
    "        y_r = tf.add(y_r, tf.divide(tf.subtract(tf.to_float(lfsize[2]), 1.0), 2.0))\n",
    "        x_r = tf.add(x_r, tf.divide(tf.subtract(tf.to_float(lfsize[3]), 1.0), 2.0))\n",
    "        \n",
    "        v_r = tf.reshape(v_r, [-1,1])\n",
    "        u_r = tf.reshape(u_r, [-1,1])\n",
    "        y_r = tf.reshape(y_r, [-1,1])\n",
    "        x_r = tf.reshape(x_r, [-1,1])\n",
    "    \n",
    "        v_r_1 = tf.cast(tf.floor(v_r), tf.int32)\n",
    "        v_r_2 = v_r_1 + 1\n",
    "        u_r_1 = tf.cast(tf.floor(u_r), tf.int32)\n",
    "        u_r_2 = u_r_1 + 1\n",
    "        y_r_1 = tf.cast(tf.floor(y_r), tf.int32)\n",
    "        y_r_2 = y_r_1 + 1\n",
    "        x_r_1 = tf.cast(tf.floor(x_r), tf.int32)\n",
    "        x_r_2 = x_r_1 + 1\n",
    "        \n",
    "        v_r_1 = tf.clip_by_value(v_r_1, 0, lfsize[0]-1)\n",
    "        v_r_2 = tf.clip_by_value(v_r_2, 0, lfsize[0]-1)\n",
    "        u_r_1 = tf.clip_by_value(u_r_1, 0, lfsize[1]-1)\n",
    "        u_r_2 = tf.clip_by_value(u_r_2, 0, lfsize[1]-1)\n",
    "        y_r_1 = tf.clip_by_value(y_r_1, 0, lfsize[2]-1)\n",
    "        y_r_2 = tf.clip_by_value(y_r_2, 0, lfsize[2]-1)\n",
    "        x_r_1 = tf.clip_by_value(x_r_1, 0, lfsize[3]-1)\n",
    "        x_r_2 = tf.clip_by_value(x_r_2, 0, lfsize[3]-1)\n",
    "        \n",
    "        #interpolate reparameterized points (quadrilinear)\n",
    "        interp_pts_1 = tf.concat([v_r_1, u_r_1, y_r_1, x_r_1], 1)\n",
    "        interp_pts_2 = tf.concat([v_r_2, u_r_1, y_r_1, x_r_1], 1)\n",
    "        interp_pts_3 = tf.concat([v_r_1, u_r_2, y_r_1, x_r_1], 1)\n",
    "        interp_pts_4 = tf.concat([v_r_1, u_r_1, y_r_2, x_r_1], 1)\n",
    "        interp_pts_5 = tf.concat([v_r_1, u_r_1, y_r_1, x_r_2], 1)\n",
    "        interp_pts_6 = tf.concat([v_r_2, u_r_2, y_r_1, x_r_1], 1)\n",
    "        interp_pts_7 = tf.concat([v_r_2, u_r_1, y_r_2, x_r_1], 1)\n",
    "        interp_pts_8 = tf.concat([v_r_2, u_r_1, y_r_1, x_r_2], 1)\n",
    "        interp_pts_9 = tf.concat([v_r_1, u_r_2, y_r_2, x_r_1], 1)\n",
    "        interp_pts_10 = tf.concat([v_r_1, u_r_2, y_r_1, x_r_2], 1)\n",
    "        interp_pts_11 = tf.concat([v_r_1, u_r_1, y_r_2, x_r_2], 1)\n",
    "        interp_pts_12 = tf.concat([v_r_2, u_r_2, y_r_2, x_r_1], 1)\n",
    "        interp_pts_13 = tf.concat([v_r_2, u_r_2, y_r_1, x_r_2], 1)\n",
    "        interp_pts_14 = tf.concat([v_r_2, u_r_1, y_r_2, x_r_2], 1)\n",
    "        interp_pts_15 = tf.concat([v_r_1, u_r_2, y_r_2, x_r_2], 1)\n",
    "        interp_pts_16 = tf.concat([v_r_2, u_r_2, y_r_2, x_r_2], 1)\n",
    "        \n",
    "        lf_r_1 = gather_4D(tf.squeeze(lf), interp_pts_1)\n",
    "        lf_r_2 = gather_4D(tf.squeeze(lf), interp_pts_2)\n",
    "        lf_r_3 = gather_4D(tf.squeeze(lf), interp_pts_3)\n",
    "        lf_r_4 = gather_4D(tf.squeeze(lf), interp_pts_4)\n",
    "        lf_r_5 = gather_4D(tf.squeeze(lf), interp_pts_5)\n",
    "        lf_r_6 = gather_4D(tf.squeeze(lf), interp_pts_6)\n",
    "        lf_r_7 = gather_4D(tf.squeeze(lf), interp_pts_7)\n",
    "        lf_r_8 = gather_4D(tf.squeeze(lf), interp_pts_8)\n",
    "        lf_r_9 = gather_4D(tf.squeeze(lf), interp_pts_9)\n",
    "        lf_r_10 = gather_4D(tf.squeeze(lf), interp_pts_10)\n",
    "        lf_r_11 = gather_4D(tf.squeeze(lf), interp_pts_11)\n",
    "        lf_r_12 = gather_4D(tf.squeeze(lf), interp_pts_12)\n",
    "        lf_r_13 = gather_4D(tf.squeeze(lf), interp_pts_13)\n",
    "        lf_r_14 = gather_4D(tf.squeeze(lf), interp_pts_14)\n",
    "        lf_r_15 = gather_4D(tf.squeeze(lf), interp_pts_15)\n",
    "        lf_r_16 = gather_4D(tf.squeeze(lf), interp_pts_16)\n",
    "        \n",
    "        v_r_1_f = tf.cast(v_r_1, tf.float32)\n",
    "        v_r_2_f = tf.cast(v_r_2, tf.float32)\n",
    "        u_r_1_f = tf.cast(u_r_1, tf.float32)\n",
    "        u_r_2_f = tf.cast(u_r_2, tf.float32)\n",
    "        y_r_1_f = tf.cast(y_r_1, tf.float32)\n",
    "        y_r_2_f = tf.cast(y_r_2, tf.float32)\n",
    "        x_r_1_f = tf.cast(x_r_1, tf.float32)\n",
    "        x_r_2_f = tf.cast(x_r_2, tf.float32)\n",
    "        \n",
    "        d_v_1 = 1.0 - (v_r - v_r_1_f)\n",
    "        d_v_2 = 1.0 - d_v_1\n",
    "        d_u_1 = 1.0 - (u_r - u_r_1_f)\n",
    "        d_u_2 = 1.0 - d_u_1\n",
    "        d_y_1 = 1.0 - (y_r - y_r_1_f)\n",
    "        d_y_2 = 1.0 - d_y_1\n",
    "        d_x_1 = 1.0 - (x_r - x_r_1_f)\n",
    "        d_x_2 = 1.0 - d_x_1\n",
    "        \n",
    "        w1 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_1), d_y_1), d_x_1)\n",
    "        w2 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_1), d_y_1), d_x_1)\n",
    "        w3 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_2), d_y_1), d_x_1)\n",
    "        w4 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_1), d_y_2), d_x_1)\n",
    "        w5 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_1), d_y_1), d_x_2)\n",
    "        w6 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_2), d_y_1), d_x_1)\n",
    "        w7 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_1), d_y_2), d_x_1)\n",
    "        w8 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_1), d_y_1), d_x_2)\n",
    "        w9 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_2), d_y_2), d_x_1)\n",
    "        w10 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_2), d_y_1), d_x_2)\n",
    "        w11 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_1), d_y_2), d_x_2)\n",
    "        w12 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_2), d_y_2), d_x_1)\n",
    "        w13 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_2), d_y_1), d_x_2)\n",
    "        w14 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_1), d_y_2), d_x_2)\n",
    "        w15 = tf.multiply(tf.multiply(tf.multiply(d_v_1, d_u_2), d_y_2), d_x_2)\n",
    "        w16 = tf.multiply(tf.multiply(tf.multiply(d_v_2, d_u_2), d_y_2), d_x_2)\n",
    "        \n",
    "        \n",
    "        lf_r = tf.add_n([w1*lf_r_1, w2*lf_r_2, w3*lf_r_3, w4*lf_r_4, w5*lf_r_5, w6*lf_r_6, w7*lf_r_7, w8*lf_r_8, \n",
    "                         w9*lf_r_9, w10*lf_r_10, w11*lf_r_11, w12*lf_r_12, w13*lf_r_13, w14*lf_r_14, w15*lf_r_15, w16*lf_r_16])\n",
    "        \n",
    "        lf_r = tf.reshape(lf_r, lfsize)\n",
    "    \n",
    "        return lf_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def lf_blur_forward_tf(lf, pts_1, pts_2, pts_3, order):\n",
    "    #tensorflow forward model to blur light field along Bezier curve motion path, with given control points\n",
    "    #assumes first control point is origin\n",
    "    #currently only implemented for orders 1 and 2\n",
    "    b = tf.zeros(lfsize)\n",
    "    for i in range(num_exp_pts):\n",
    "        t = np.true_divide(i, num_exp_pts).astype(np.float32)\n",
    "        coord_2 = lambda: tf.multiply((2.0*t*(1.0-t)).astype(np.float32), pts_1) + tf.multiply(tf.square(t), pts_2) #quadratic\n",
    "        coord_3 = lambda: tf.multiply((3.0*t*(1.0-t)*(1.0-t)).astype(np.float32), pts_1) +tf.multiply((3.0*t*t*(1.0-t)).astype(np.float32), pts_2) + tf.multiply(tf.pow(t, 3), pts_3) #cubic\n",
    "        coord = tf.cond(tf.equal(order, tf.constant(2)), coord_2, coord_3)\n",
    "        b = tf.add(b, lf_reparam_tf(lf, tf.squeeze(coord)))\n",
    "    b = tf.divide(b, tf.cast(num_exp_pts, tf.float32))\n",
    "    \n",
    "    return b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def data_loss(lf_observed, lf_sharp, pts_1, pts_2, pts_3, order):\n",
    "    #data term (l2 norm of difference between observed blurred light field and forward model predicted light field)\n",
    "    return tf.reduce_mean(tf.squared_difference(lf_blur_forward_tf(lf_sharp, pts_1, pts_2, pts_3, order), lf_observed))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tv_loss_s(x):\n",
    "    #spatial total variation loss (l1 norm of spatial derivatives)\n",
    "    temp = x[:,:,0:lfsize[2]-1,0:lfsize[3]-1]\n",
    "    dy = (x[:,:,1:lfsize[2],0:lfsize[3]-1] - temp)\n",
    "    dx = (x[:,:,0:lfsize[2]-1,1:lfsize[3]] - temp)\n",
    "    l_1 = tf.reduce_mean(tf.abs(dy)+tf.abs(dx))\n",
    "    return l_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def tv_loss_a(x):\n",
    "    #angular total variation loss (l1 norm of angular derivatives)\n",
    "    temp = x[0:lfsize[0]-1,0:lfsize[1]-1,:,:]\n",
    "    dv = (x[1:lfsize[0],0:lfsize[1]-1,:,:] - temp)\n",
    "    du = (x[0:lfsize[0]-1,1:lfsize[1],:,:] - temp)\n",
    "    l_1 = tf.reduce_mean(tf.abs(dv)+tf.abs(du))\n",
    "    return l_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sp_loss_s(x, eps):\n",
    "    #gradual non-convex approximation of l0 norm of spatial derivatives\n",
    "    temp = x[:,:,0:lfsize[2]-1,0:lfsize[3]-1]\n",
    "    dy = tf.abs(x[:,:,1:lfsize[2],0:lfsize[3]-1] - temp)\n",
    "    dx = tf.abs(x[:,:,0:lfsize[2]-1,1:lfsize[3]] - temp)\n",
    "    dy_c = tf.clip_by_value((1/(tf.square(eps)))*tf.square(dy), 0.0, 1.0)\n",
    "    dx_c = tf.clip_by_value((1/(tf.square(eps)))*tf.square(dx), 0.0, 1.0)\n",
    "    return tf.reduce_mean(dy_c+dx_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sp_loss_a(x, eps):\n",
    "    #gradual non-convex approximation of l0 norm of angular derivatives\n",
    "    temp = x[0:lfsize[0]-1,0:lfsize[1]-1,:,:]\n",
    "    dv = tf.abs(x[1:lfsize[0],0:lfsize[1]-1,:,:] - temp)\n",
    "    du = tf.abs(x[0:lfsize[0]-1,1:lfsize[1],:,:] - temp)\n",
    "    dv_c = tf.clip_by_value((1/(tf.square(eps)))*tf.square(dv), 0.0, 1.0)\n",
    "    du_c = tf.clip_by_value((1/(tf.square(eps)))*tf.square(du), 0.0, 1.0)\n",
    "    return tf.reduce_mean(dv_c+du_c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lf_subproblem(pts_1, pts_2, pts_3, lf_blur, num_iters, lam_lf, eta_lf, order, record):\n",
    "    #solve for latent sharp light field, given observed blurred light field and control points\n",
    "    lam = tf.constant(lam_lf, dtype=tf.float32)\n",
    "    lf_blur_placeholder = tf.placeholder(tf.float32, shape=[lfsize[0], lfsize[1], lfsize[2], lfsize[3]])\n",
    "    pts_1_placeholder = tf.placeholder(tf.float32, shape=[1,3])\n",
    "    pts_2_placeholder = tf.placeholder(tf.float32, shape=[1,3])\n",
    "    pts_3_placeholder = tf.placeholder(tf.float32, shape=[1,3])\n",
    "    lf_var = tf.Variable(tf.constant(lf_blur, dtype=tf.float32), name='lf_var')\n",
    "    data_term = data_loss(lf_blur_placeholder, lf_var, pts_1_placeholder, pts_2_placeholder, pts_3_placeholder, order)\n",
    "    prior_loss = (lam*tv_loss_s(lf_var)) + (lam*tv_loss_a(lf_var))\n",
    "    full_loss = data_term + prior_loss\n",
    "    train_step_lf = tf.train.AdamOptimizer(learning_rate=eta_lf).minimize(full_loss, var_list=[lf_var])\n",
    "    \n",
    "    iter_full_loss = np.zeros((num_iters))\n",
    "\n",
    "    with tf.Session() as sess:\n",
    "        sess.run(tf.global_variables_initializer())\n",
    "        for i in range(num_iters):   \n",
    "            print ('lf subproblem iteration %i'%(i))\n",
    "            if (record):\n",
    "                curr_full_loss = full_loss.eval(feed_dict={lf_blur_placeholder:lf_blur, pts_1_placeholder:pts_1, pts_2_placeholder:pts_2, pts_3_placeholder:pts_3})\n",
    "                iter_full_loss[i] = curr_full_loss        \n",
    "            sess.run([train_step_lf], feed_dict={lf_blur_placeholder:lf_blur, pts_1_placeholder:pts_1, pts_2_placeholder:pts_2, pts_3_placeholder:pts_3})\n",
    "            lf_var = tf.clip_by_value(lf_var, 0.0, 1.0)\n",
    "        return lf_var.eval(), iter_full_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate Motion Blur"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#load light field\n",
    "temp = sp.io.loadmat('Data/Synthetic/lf_12.mat')\n",
    "lf = normalize_lf(np.sum(np.array(temp['lf']), axis=4)).astype(np.float32)\n",
    "lf_r = normalize_lf(np.array(temp['lf'])[:,:,:,:,0]).astype(np.float32)\n",
    "lf_g = normalize_lf(np.array(temp['lf'])[:,:,:,:,1]).astype(np.float32)\n",
    "lf_b = normalize_lf(np.array(temp['lf'])[:,:,:,:,2]).astype(np.float32)\n",
    "#simulate motion blurred light field\n",
    "pts_1_blur = np.array([-2.0, 2.0, 1.0]).astype(np.float32)\n",
    "pts_2_blur = np.array([-4.0, -4.0, 2.0]).astype(np.float32)\n",
    "pts_3_blur = np.array([-6.0, 6.0, 3.0]).astype(np.float32)\n",
    "with tf.Session() as sess:\n",
    "    lf_blur = lf_blur_forward_tf(lf, pts_1_blur, pts_2_blur, pts_3_blur, tf.constant(path_order)).eval()\n",
    "    lf_blur_r = lf_blur_forward_tf(lf_r, pts_1_blur, pts_2_blur, pts_3_blur, tf.constant(path_order)).eval()\n",
    "    lf_blur_g = lf_blur_forward_tf(lf_g, pts_1_blur, pts_2_blur, pts_3_blur, tf.constant(path_order)).eval()\n",
    "    lf_blur_b = lf_blur_forward_tf(lf_b, pts_1_blur, pts_2_blur, pts_3_blur, tf.constant(path_order)).eval()\n",
    "tf.reset_default_graph()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Real Motion Blur (Illum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#load light field\n",
    "temp = sp.io.loadmat('Data/Real/illum_73.mat')\n",
    "lf = normalize_lf(np.sum(np.array(temp['lf']), axis=4)).astype(np.float32)\n",
    "lf_r = normalize_lf(np.array(temp['lf'])[:,:,:,:,0]).astype(np.float32)\n",
    "lf_g = normalize_lf(np.array(temp['lf'])[:,:,:,:,1]).astype(np.float32)\n",
    "lf_b = normalize_lf(np.array(temp['lf'])[:,:,:,:,2]).astype(np.float32)\n",
    "\n",
    "lf_blur = np.copy(lf)\n",
    "lf_blur_r = np.copy(lf_r)\n",
    "lf_blur_g = np.copy(lf_g)\n",
    "lf_blur_b = np.copy(lf_b)\n",
    "\n",
    "#placeholder blur control points (to avoid error when saving results)\n",
    "pts_1_blur = np.array([0, 0, 0]).astype(np.float32) \n",
    "pts_2_blur = np.array([0, 0, 0]).astype(np.float32) \n",
    "pts_3_blur = np.array([0, 0, 0]).astype(np.float32) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blind Motion Deblurring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#set up tensorflow graph\n",
    "lam = tf.placeholder(tf.float32)\n",
    "eps = tf.placeholder(tf.float32)\n",
    "order = tf.placeholder(tf.int32)\n",
    "lf_blur_placeholder = tf.placeholder(tf.float32, shape=[lfsize[0], lfsize[1], lfsize[2], lfsize[3]])\n",
    "lf_var = tf.Variable(tf.constant(lf_blur), name='lf_var')\n",
    "pts_1_var = tf.Variable(tf.zeros([1,3]), name='pts_1_var') #first non-origin control point\n",
    "pts_2_var = tf.Variable(tf.zeros([1,3]), name='pts_2_var') #second non-origin control point\n",
    "pts_3_var = tf.Variable(tf.zeros([1,3]), name='pts_3_var') #third non-origin control point\n",
    "with tf.variable_scope('loss'):\n",
    "    data_term = data_loss(lf_blur_placeholder, lf_var, pts_1_var, pts_2_var, pts_3_var, order)\n",
    "    prior_loss = (lam*sp_loss_s(lf_var, eps)) + (lam*sp_loss_a(lf_var, eps))\n",
    "    full_loss = data_term + prior_loss\n",
    "with tf.variable_scope('train'):\n",
    "    train_step_lf = tf.train.AdamOptimizer(learning_rate=eta_lf).minimize(full_loss, var_list=[lf_var])\n",
    "    train_step_pts = tf.train.AdamOptimizer(learning_rate=eta_pts).minimize(full_loss, var_list=[pts_1_var, pts_2_var, pts_3_var])\n",
    "    train_step = tf.group(train_step_lf, train_step_pts)\n",
    "    \n",
    "#losses to record\n",
    "iter_data_loss = np.zeros((num_iters))\n",
    "iter_prior_loss = np.zeros((num_iters))\n",
    "iter_full_loss = np.zeros((num_iters))\n",
    "lf_mse = np.zeros((num_iters)) #only useful for synthetic examples\n",
    "\n",
    "with tf.Session() as sess:\n",
    "    sess.run(tf.global_variables_initializer())\n",
    "    lam_curr = lam_init\n",
    "    eps_curr = eps_init\n",
    "    for i in range(num_iters):\n",
    "        #print interation information\n",
    "        print ('iteration %i'%(i))\n",
    "        print ('current points')\n",
    "        print pts_1_var.eval()\n",
    "        print pts_2_var.eval()\n",
    "        print pts_3_var.eval()\n",
    "        #calculate losses\n",
    "        curr_data_loss = data_term.eval(feed_dict={lf_blur_placeholder:lf_blur, lam:lam_curr, eps:eps_curr, order:path_order})\n",
    "        curr_prior_loss = prior_loss.eval(feed_dict={lf_blur_placeholder:lf_blur, lam:lam_curr, eps:eps_curr, order:path_order})\n",
    "        curr_full_loss = full_loss.eval(feed_dict={lf_blur_placeholder:lf_blur, lam:lam_curr, eps:eps_curr, order:path_order})\n",
    "        iter_data_loss[i] = curr_data_loss\n",
    "        iter_prior_loss[i] = curr_prior_loss\n",
    "        iter_full_loss[i] = curr_full_loss\n",
    "        lf_mse[i] = np.mean(np.square(lf - lf_var.eval()))\n",
    "        #run tensorflow session (optimization step)\n",
    "        sess.run([train_step], feed_dict={lf_blur_placeholder:lf_blur, lam:lam_curr, eps:eps_curr, order:path_order})\n",
    "        #project latent sharp light field to non-negative values    \n",
    "        lf_var = tf.clip_by_value(lf_var, 0.0, 1.0)\n",
    "        #decay regularization weight and increase regularization non-convexity\n",
    "        lam_curr = np.clip(lam_curr*lam_decay, lam_min, lam_init)\n",
    "        eps_curr = np.clip(eps_curr*eps_decay, eps_min, eps_init)\n",
    "    #evaluate final latent sharp light field and control points\n",
    "    lf_final = lf_var.eval()\n",
    "    pts_1_final = pts_1_var.eval()\n",
    "    pts_2_final = pts_2_var.eval()\n",
    "    pts_3_final = pts_3_var.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#plot losses\n",
    "plt.figure()\n",
    "plt.plot(iter_full_loss[0,:])\n",
    "plt.figure()\n",
    "plt.plot(iter_full_loss[1,:])\n",
    "plt.figure()\n",
    "plt.plot(iter_data_loss[0,:])\n",
    "plt.figure()\n",
    "plt.plot(iter_data_loss[1,:])\n",
    "plt.figure()\n",
    "plt.plot(iter_prior_loss[0,:])\n",
    "plt.figure()\n",
    "plt.plot(iter_prior_loss[1,:])\n",
    "plt.figure()\n",
    "plt.plot(lf_mse[0,:])\n",
    "plt.figure()\n",
    "plt.plot(lf_mse[1,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#solve LF subproblem for each color channel for final light field estimate\n",
    "tf.reset_default_graph()\n",
    "lf_deblur_r, loss_deblur_r = lf_subproblem(pts_1_final, pts_2_final, pts_3_final, lf_blur_r, iters_lf, lam_lf, eta_lf, path_order, False)\n",
    "tf.reset_default_graph()\n",
    "lf_deblur_g, loss_deblur_g = lf_subproblem(pts_1_final, pts_2_final, pts_3_final, lf_blur_g, iters_lf, lam_lf, eta_lf, path_order, False)\n",
    "tf.reset_default_graph()\n",
    "lf_deblur_b, loss_deblur_b = lf_subproblem(pts_1_final, pts_2_final, pts_3_final, lf_blur_b, iters_lf, lam_lf, eta_lf, path_order, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#save inputs and results\n",
    "deblur_dict = {'lf_in_r':lf_blur_r, 'lf_in_g':lf_blur_g, 'lf_in_b':lf_blur_b, \n",
    "             'lf_out_r':lf_deblur_r, 'lf_out_g':lf_deblur_g, 'lf_out_b':lf_deblur_b,\n",
    "             'pts_1_in':pts_1_blur, 'pts_2_in':pts_2_blur, 'pts_3_in':pts_3_blur,\n",
    "             'pts_1_out':pts_1_final, 'pts_2_out':pts_2_final, 'pts_3_out':pts_3_final}\n",
    "sp.io.savemat('deblur_results.mat', deblur_dict)"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}