Linear Programming | Image Denoising

Linear Programming | Image Denoising

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Programming | Image Denoising"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import cv2\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import time\n",
    "from scipy import optimize\n",
    "from __future__ import division"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to display the image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def img_display(img, name):\n",
    "    plt.figure()\n",
    "    plt.imshow(img, cmap='gray', interpolation='bicubic')\n",
    "    plt.title(name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Function to create noisy gaussian image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def noisy(image):\n",
    "    row,col= image.shape\n",
    "    mean = 0\n",
    "    var = 0.1\n",
    "    sigma = var**0.5\n",
    "    gauss = np.random.normal(mean,sigma,(row,col))\n",
    "    gauss = gauss.reshape(row,col)\n",
    "    noisy = image + gauss\n",
    "    return noisy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://glimpglobe.com/proj/lp/imgtrue.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img1 = cv2.imread('image.png', cv2.IMREAD_GRAYSCALE)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "img1 = np.float32(img1)\n",
    "img_noisy = noisy(img1)\n",
    "\n",
    "x = img1.shape[0]\n",
    "y = img1.shape[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noisy Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://glimpglobe.com/proj/lp/imgnoisy.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Just some matrix operations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "alpha = np.zeros((20, 20))\n",
    "alpha.fill(1/(x*y))\n",
    "alphab = np.zeros((20, 20))\n",
    "alphab.fill(1/(x*y))\n",
    "\n",
    "b1 = img_noisy.ravel() \n",
    "b2 = img_noisy.ravel()\n",
    "b = -1 * np.concatenate((b1, b2), axis=0)\n",
    "\n",
    "c = np.concatenate((alpha, alphab), axis=0)\n",
    "c = c.ravel()\n",
    "c = np.append(c, 0)\n",
    "c = np.append(c, 0.4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "A=np.zeros((2*x*y, (2*x*y)+1))\n",
    "A = np.zeros((2*x*y, 2*(x*y+1)))\n",
    "A[:x*y, 2*x*y] = -img1.ravel()\n",
    "\n",
    "A[x*y-1:-1, 2*(x*y)] = img1.ravel()\n",
    "\n",
    "n = 1\n",
    "N = x*y+1\n",
    "\n",
    "alphar = alpha.ravel()\n",
    "alphabr = alphab.ravel()\n",
    "\n",
    "for i in range(0,x*y):\n",
    "    A[i][n]=alphar[i]\n",
    "    A[i][N]=-1*alphabr[i]\n",
    "    n = n+1\n",
    "    N = N+1\n",
    "\n",
    "n = 1\n",
    "N = x*y+1\n",
    "\n",
    "for i in range(0,x*y):\n",
    "    A[i][n]=-1*alphar[i]\n",
    "    A[i][N]=alphabr[i]\n",
    "    n = n+1\n",
    "    N = N+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Python Scipy library includes linprog function which takes in Coefficients of the linear objective function to be minimized, which here is \"c\". A_ub which is a matrix that after multiplied gets the upper-bound inequality, b_ub is an array for the upper-bound in 1-D. A_eq gets the values of the equality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "res = optimize.linprog(c, A_ub=A, b_ub=b, bounds=(0,800),options=dict(bland=True, tol=1e-15))\n",
    "\n",
    "alpha=res[0:(x*y*2)];\n",
    "alpha=alpha[0:x*y]-alpha[(x*y):x*y*2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ Xn+\\alpha_{n} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "alpha=alpha+img_noisy\n",
    "# img_noisy = Xn\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img_est = np.reshape(alpha, (20, 20))\n",
    "img_display(img_est, \"Image Estimate\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final Result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://glimpglobe.com/proj/lp/imestt.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "img_diff = img1 - img_est\n",
    "img_display(img_est, \"Image Result\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error Estimate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://glimpglobe.com/proj/lp/imgest.png\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}