TVDirectionalTest.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TV Minimization\n",
    "\n",
    "\n",
    "In this example $Y=AX+N$ simulation enviroment will be created.\n",
    "In this equation:\n",
    "- Y is the observation\n",
    "- A is system geometry parameters\n",
    "- X is varible (the thing which is observed )\n",
    "- N is randomly distrubuted zero mean gaussian noise\n",
    "\n",
    "Size of this vector and matrices are:\n",
    "$$ Y_{mx1}=A_{mxn}X_{nx1} + N_{mx1} $$\n",
    "\n",
    "\n",
    "now contiune with code.\n",
    "\n",
    "\n",
    "## 1. Configuration\n",
    "\n",
    "### Import libraries\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "from math import sqrt\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image\n",
    "%matplotlib inline \n",
    "\n",
    "plt.rcParams['figure.figsize'] = [16, 10]\n",
    "plt.gray()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def normalize(a):\n",
    "    return a/255\n",
    "\n",
    "def denormalize(a):\n",
    "    return a*255\n",
    "\n",
    "# image to matrix\n",
    "def i2m(a):\n",
    "    return np.array(a)/255\n",
    "\n",
    "# matrix to image\n",
    "def m2i(a):\n",
    "    return Image.fromarray(a*255).convert(\"L\")\n",
    "\n",
    "def rmse(X1, X2):\n",
    "    return np.sqrt(np.mean((X1-X2)**2))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPS = 0.0000000000000001 # epsilon\n",
    "BETA =  0.2 #\n",
    "LEARNING_RATE = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Init\n",
    "\n",
    "To simulate the problem the following equation is defined\n",
    "\n",
    "$ Y = X + N $\n",
    "\n",
    "where:\n",
    "- X is the real & non-noisy image ( It is the varible which should be estimated )\n",
    "- N is random gaussion noise\n",
    "- Y is the observation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Xreal = i2m(  Image.open('test/zebra.jpg').convert('L'))#.resize((256,256), Image.ANTIALIAS))#.rotate(-90) )\n",
    "Xreal = i2m(  Image.open('test/lines.jpg').convert('L'))\n",
    "pixelSizeX, pixelSizeY = np.shape(Xreal)\n",
    "noise = np.random.rand(pixelSizeX, pixelSizeY) * 0.2\n",
    "\n",
    "Y = Xreal + noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m2i(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. TV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define TV\n",
    "\n",
    "Tv\n",
    "\n",
    "$\n",
    "T V_{2 D}(X)=\\sum_{i}^{K} \\sum_{j}^{L} \\sqrt{\\left(X_{i, j}-X_{i-1, j}\\right)^{2}+\\left(X_{i, j}-X_{i, j-1}\\right)^{2}}\n",
    "$\n",
    "\n",
    "Gradient of tv\n",
    "\n",
    "$Grad(X)=\\frac{2\\left(X_{i, j}-X_{i-1, j}\\right)+2\\left(X_{i j}-X_{i j-1}\\right)}{\\sqrt{\\left(X_{i, j}-X_{i-1, j}\\right)^{2}+\\left(X_{i, j}-X_{i, j-1}\\right)^{2}+\\varepsilon}}\n",
    "-\\frac{2\\left(X_{i+1, j}-X_{i, j}\\right)}{\\sqrt{\\left(X_{i+1, j}-X_{i, j}\\right)^{2}+\\left(X_{i+1, j}-X_{i+1, j-1}\\right)^{2}+\\varepsilon}}\n",
    "-\\frac{2\\left(X_{i, j+1}-X_{i, j}\\right)}{\\sqrt{\\left(X_{i, j+1}-X_{i, j}\\right)^{2}+\\left(X_{i, j+1}-X_{i-1, j+1}\\right)^{2}+\\varepsilon}}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classic TV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tvNorm(x):\n",
    "    \"\"\"Computes the total variation norm and its gradient. From jcjohnson/cnn-vis.\"\"\"\n",
    "    x_diff = x - np.roll(x, -1, axis=1)\n",
    "    y_diff = x - np.roll(x, -1, axis=0)\n",
    "    grad_norm2 = x_diff**2 + y_diff**2 + EPS\n",
    "    norm = np.sum(np.sqrt(grad_norm2))\n",
    "    return norm\n",
    "\n",
    "def tvGrad(x):\n",
    "    \"\"\"Computes the total variation norm and its gradient. From jcjohnson/cnn-vis.\"\"\"\n",
    "    x_diff = x - np.roll(x, -1, axis=1)\n",
    "    y_diff = x - np.roll(x, -1, axis=0)\n",
    "    grad_norm2 = x_diff**2 + y_diff**2 + EPS\n",
    "    dgrad_norm = 0.5 / np.sqrt(grad_norm2)\n",
    "    dx_diff = 2 * x_diff * dgrad_norm\n",
    "    dy_diff = 2 * y_diff * dgrad_norm\n",
    "    grad = dx_diff + dy_diff\n",
    "    grad[:, 1:] -= dx_diff[:, :-1]\n",
    "    grad[1:, :] -= dy_diff[:-1, :]\n",
    "    return grad\n",
    "\n",
    "\n",
    "def l2Norm(x):\n",
    "    \"\"\"Computes 1/2 the square of the L2-norm and its gradient.\"\"\"\n",
    "    return np.sum(x**2)\n",
    "\n",
    "def l2NormGrad(x):\n",
    "    return -2 * x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Directional TV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def directionalTVNormM(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    tv = 0\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            payda1 = (XM[i,j] - XM[i+1,j])**2\n",
    "            payda2 = (XM[i,j] - XM[i,j+1])**2\n",
    "            payda3 = (XM[i,j-1] - XM[i,j+2])**2\n",
    "            tv += sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "    return tv\n",
    "\n",
    "def directionalTVNorm(x):\n",
    "    payda1 = x - np.roll(x, +1, axis=0)\n",
    "    payda2 = x - np.roll(x, +1, axis=1)\n",
    "    payda3 = np.roll(x, -1, axis=1) - np.roll(x, +2, axis=1)\n",
    "    return np.sum( np.sqrt( payda1**2 + payda2**2 + payda3**2) )\n",
    "\n",
    "\n",
    "def gradMA(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    gradX = np.zeros((M,N))\n",
    "    # XM(4:end-3,4:end-3)=X\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            pay = 0.5 * -2 * ( XM[i,j-3] - XM[i,j] )\n",
    "            payda1 = (XM[i,j-2] - XM[i+1,j-2])**2\n",
    "            payda2 = (XM[i,j-2] - XM[i,j-1])**2\n",
    "            payda3 = (XM[i,j-3] - XM[i,j])**2\n",
    "            payda = sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "            gradX[i,j] = pay / payda\n",
    "    return gradX\n",
    "\n",
    "def gradA(x):\n",
    "    pay = -1 * ( x - np.roll(x, -3, axis=1) )\n",
    "    payda1 = np.roll(x, -2, axis=1) - np.roll(np.roll(x, 1, axis=0), -2, axis=1 )\n",
    "    payda2 = np.roll(x, -2, axis=1) - np.roll(x, -1, axis=1)\n",
    "    payda3 = np.roll(x, -3, axis=1) - x\n",
    "    payda = np.sqrt( payda1**2 + payda2**2 + payda3**2 + EPS )\n",
    "    gradX = pay / payda\n",
    "    return gradX\n",
    "    \n",
    "\n",
    "def gradMB(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    gradX = np.zeros((M,N))\n",
    "    # XM(4:end-3,4:end-3)=X\n",
    "\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            pay = 0.5 * -2 * ( XM[i,j-1] - XM[i,j] )\n",
    "            payda1 = (XM[i,j-1] - XM[i+1,j-1])**2\n",
    "            payda2 = (XM[i,j-1] - XM[i,j])**2\n",
    "            payda3 = (XM[i,j-2] - XM[i,j+1])**2\n",
    "            payda = sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "            gradX[i,j] = pay / payda\n",
    "    return gradX\n",
    "\n",
    "def gradB(x):\n",
    "    pay = -1 * ( np.roll(x, -1, axis=1) -  x )\n",
    "    payda1 = np.roll(x, -1, axis=1) - np.roll(np.roll(x, 1, axis=0), 1, axis=1 )\n",
    "    payda2 = np.roll(x, -1, axis=1) - x\n",
    "    payda3 = np.roll(x, -2, axis=1) - np.roll(x, +1, axis=1)\n",
    "    payda = np.sqrt( payda1**2 + payda2**2 + payda3**2 + EPS )\n",
    "    gradX = pay / payda\n",
    "    return gradX\n",
    "\n",
    "\n",
    "def gradMC(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    gradX = np.zeros((M,N))\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            pay = 0.5 * 2 * ( XM[i,j] - XM[i+1,j] + XM[i,j] - XM[i,j+1] )\n",
    "            payda1 = (XM[i,j] - XM[i+1,j])**2\n",
    "            payda2 = (XM[i,j] - XM[i,j+1])**2\n",
    "            payda3 = (XM[i,j-1] - XM[i,j+2])**2\n",
    "            payda = sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "            gradX[i,j] = pay / payda\n",
    "    return gradX\n",
    "\n",
    "def gradC(x):\n",
    "    pay =  2 * x - np.roll(x, +1, axis=0) - np.roll(x, +1, axis=1)\n",
    "    payda1 = x - np.roll(x, +1, axis=0)\n",
    "    payda2 = x - np.roll(x, +1, axis=1)\n",
    "    payda3 = np.roll(x, -1, axis=1) - np.roll(x, +2, axis=1)\n",
    "    payda = np.sqrt( payda1**2 + payda2**2 + payda3**2 + EPS )\n",
    "    gradX = pay / payda\n",
    "    return gradX\n",
    "\n",
    "\n",
    "\n",
    "def gradMD(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    gradX = np.zeros((M,N))\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            pay = 0.5 * 2 * ( XM[i,j] - XM[i,j+3] )\n",
    "            payda1 = (XM[i,j+1] - XM[i+1,j+1])**2\n",
    "            payda2 = (XM[i,j+1] - XM[i,j+2])**2\n",
    "            payda3 = (XM[i,j] - XM[i,j+3])**2\n",
    "            payda = sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "            gradX[i,j] = pay / payda\n",
    "    return gradX\n",
    "\n",
    "\n",
    "def gradD(x):\n",
    "    pay = x - np.roll(x, +3, axis=1)\n",
    "    payda1 = np.roll(x, +1, axis=1) - np.roll(np.roll(x, +1, axis=0), +1, axis=1 )\n",
    "    payda2 = np.roll(x, +1, axis=1) - np.roll(x, +2, axis=1)\n",
    "    payda3 = x - np.roll(x, +3, axis=1)\n",
    "    payda = np.sqrt( payda1**2 + payda2**2 + payda3**2 + EPS )\n",
    "    gradX = pay / payda\n",
    "    return gradX\n",
    "\n",
    "def gradME(X):\n",
    "    M,N = np.shape(X)\n",
    "    XM = X\n",
    "    gradX = np.zeros((M,N))\n",
    "    # XM(4:end-3,4:end-3)=X\n",
    "\n",
    "    for i in range(3, M-3):\n",
    "        for j in range(3, N-3):\n",
    "            pay = 0.5 * -2 * ( XM[i-1,j] - XM[i,j] )\n",
    "            payda1 = (XM[i-1,j] - XM[i,j])**2\n",
    "            payda2 = (XM[i-1,j] - XM[i-1,j+1])**2\n",
    "            payda3 = (XM[i-1,j-1] - XM[i-1,j+2])**2\n",
    "            payda = sqrt( payda1 + payda2 + payda3 + EPS)\n",
    "            gradX[i,j] = pay / payda\n",
    "    return gradX\n",
    "\n",
    "def gradE(x):\n",
    "    pay = -1 * ( np.roll(x, -1, axis=0) - x )\n",
    "    payda1 = np.roll(x, -1, axis=0) - x\n",
    "    payda2 = np.roll(x, -1, axis=0) - np.roll(np.roll(x, -1, axis=0), +1, axis=1 )\n",
    "    payda3 = np.roll(np.roll(x, -1, axis=0), -1, axis=1 ) - np.roll(np.roll(x, -1, axis=0), +2, axis=1 )\n",
    "    payda = np.sqrt( payda1**2 + payda2**2 + payda3**2 + EPS )\n",
    "    gradX = pay / payda\n",
    "    return gradX\n",
    "\n",
    "\n",
    "def directionalTVGradM(X):\n",
    "    return gradMA(X)+gradMB(X)+gradMC(X)+gradMD(X)+gradME(X)\n",
    "\n",
    "def directionalTVGrad(X):\n",
    "    return gradA(X)+gradB(X)+gradC(X)+gradD(X)+gradE(X)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.imshow( directionalTVGrad(Xreal))\n",
    "plt.figure()\n",
    "#plt.imshow( gradMA(Xreal) )\n",
    "#(gradA(Xreal) - gradMA(Xreal)).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit directionalTVGrad(Xreal)\n",
    "%timeit directionalTVNorm(Xreal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Cost functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classic TV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Cost Function\n",
    "\n",
    "$\\hat{X}=\\underset{X}{\\arg \\min }\\left[\\|Y- X\\|_{2}+\\beta T V(X)\\right]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tvCost(Y, X, beta):\n",
    "    return l2Norm(Y-X) + beta * tvNorm(X)\n",
    "\n",
    "def tvCostGrad(Y, X, beta):\n",
    "    return l2NormGrad(Y-X) + beta * tvGrad(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Directional TV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def directionalTVCost(Y, X, beta):\n",
    "    return l2Norm(Y-X) + beta * directionalTVNorm(X)\n",
    "\n",
    "def directionalTVCostGrad(Y, X, beta):\n",
    "    return l2NormGrad(Y-X) + beta * directionalTVGrad(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [16, 10]\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(121)\n",
    "ax.title.set_text('Original Image Variation on X Axis')\n",
    "ax.imshow(m2i(Xreal-np.roll(Xreal, -1, axis=1)))\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "ax.title.set_text('Noisy Image Variation on X Axis')\n",
    "ax.imshow(m2i(Y-np.roll(Y, -1, axis=1)))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [16, 10]\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(121)\n",
    "ax.title.set_text('Original Image')\n",
    "ax.imshow(m2i(directionalTVGrad(Xreal)))\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "ax.title.set_text('Noisy Image')\n",
    "ax.imshow(m2i(directionalTVGrad(Y)))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Minimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient Descent\n",
    "\n",
    "Generic form of gradient descent.\n",
    "\n",
    "$X^{k+1} = X^k - \\alpha \\nabla Cost$\n",
    "\n",
    "Gradient descent method for tv\n",
    "\n",
    "$X^{k+1} = X^k - \\alpha \\nabla( |Y- X\\|_{2}+\\beta T V(X) ) $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sum all up for TV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimize(costFunction=None, costFunctionGrad=None, \\\n",
    "             real=None, initial=None, \\\n",
    "             beta=0.3, learningrate=0.001, maxIter=1000, finishCondition=10):\n",
    "    costHistory = []\n",
    "    rmseHistory = [] \n",
    "    Xk = initial\n",
    "    for i in range(maxIter):\n",
    "        if i%30 == 0:\n",
    "            print(i, costFunction(initial, Xk, beta))\n",
    "        Xnext = Xk - learningrate * costFunctionGrad(initial, Xk, beta)\n",
    "        costHistory.append( costFunction(initial, Xk, beta) )\n",
    "        rmseHistory.append( rmse(real, Xk) )\n",
    "        # init for next iteration\n",
    "        if  costFunction(initial, Xk, beta) - costFunction(initial, Xnext, beta) < finishCondition :\n",
    "            break\n",
    "        Xk = Xnext\n",
    "    return Xnext, costHistory, rmseHistory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def minimizeFast(costFunction=None, costFunctionGrad=None, \\\n",
    "                 real=None, initial=None, \\\n",
    "                 beta=0.3, learningrate=0.001, maxIter=1000, finishCondition=10):\n",
    "    e = 0.1\n",
    "    change = 0\n",
    "    costHistory = []\n",
    "    changeHistory = []\n",
    "    rmseHistory = []\n",
    "    Xk = np.copy(initial)\n",
    "    Xnext = np.copy(initial)\n",
    "    Xold = np.copy(initial)\n",
    "    innerLoopCount = 0\n",
    "    outerLoopCount = 0\n",
    "    while True:\n",
    "        costK = EPS\n",
    "        costNext = 0\n",
    "        n = 0\n",
    "        Xk = Xold # discard last iteration\n",
    "        while ( costNext< costK ):\n",
    "            innerLoopCount += 1\n",
    "            Xnext = Xk - (2**n) * learningrate * costFunctionGrad(initial, Xk, beta)\n",
    "            costK = costFunction(initial, Xk, beta)\n",
    "            costNext = costFunction(initial, Xnext, beta)\n",
    "            print(n, costK, costNext, innerLoopCount)\n",
    "            Xold = np.copy(Xk)\n",
    "            Xk = np.copy(Xnext)\n",
    "            n = n + 1\n",
    "            \n",
    "        costHistory.append( costFunction(initial, Xold, beta) )\n",
    "        rmseHistory.append( rmse(real, Xold) )\n",
    "        change = np.linalg.norm((Xk-Xold),2);\n",
    "        #changeHistory.append(change)\n",
    "        outerLoopCount += 1\n",
    "        print (outerLoopCount, \"change\", change)\n",
    "        if ( change < 0.01 or innerLoopCount > maxIter or abs(costK-costNext) < finishCondition ): \n",
    "            break\n",
    "    return Xnext, costHistory, rmseHistory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "imTV, chistoryTV, rhistoryTV = minimize(\n",
    "    costFunction=tvCost,\n",
    "    costFunctionGrad=tvCostGrad,\n",
    "    real=Xreal,\n",
    "    learningrate=0.001,\n",
    "    initial=Y,\n",
    "    beta=0.5,\n",
    "    finishCondition=1\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(imTV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "imTVD, chistoryTVD, rhistoryTVD = minimize(\n",
    "    costFunction=directionalTVCost,\n",
    "    costFunctionGrad=directionalTVCostGrad,\n",
    "    real=Xreal,\n",
    "    learningrate=0.001,\n",
    "    initial=Y,\n",
    "    beta=0.5,\n",
    "    finishCondition=1\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(imTVD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "plt.figure()\n",
    "plt.title(\"Cost Over Iterations\")\n",
    "plt.plot(chistoryTV)\n",
    "plt.plot(chistoryTVD)\n",
    "plt.figure()\n",
    "plt.title(\"RMSE Over Iterations\")\n",
    "plt.plot(rhistoryTV)\n",
    "plt.plot(rhistoryTVD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(221)\n",
    "ax.title.set_text('Original Image')\n",
    "ax.imshow(m2i(Xreal))\n",
    "\n",
    "ax = fig.add_subplot(222)\n",
    "ax.title.set_text('Noisy Image')\n",
    "ax.imshow(m2i(Y))\n",
    "\n",
    "ax = fig.add_subplot(223)\n",
    "ax.title.set_text('TV')\n",
    "ax.imshow(m2i(imTV))\n",
    "\n",
    "ax = fig.add_subplot(224)\n",
    "ax.title.set_text('TV Directional')\n",
    "ax.imshow(m2i(imTVD))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "\n",
    "ax = fig.add_subplot(231)\n",
    "ax.title.set_text('Original Image')\n",
    "ax.imshow(m2i(Xreal))\n",
    "\n",
    "ax = fig.add_subplot(232)\n",
    "ax.title.set_text('Noisy Image')\n",
    "ax.imshow(m2i(Y))\n",
    "\n",
    "ax = fig.add_subplot(233)\n",
    "ax.title.set_text('TV B=0.3')\n",
    "ax.imshow(m2i(denoiseTV(Y, 0.3)))\n",
    "\n",
    "ax = fig.add_subplot(234)\n",
    "ax.title.set_text('TV B=0.7')\n",
    "ax.imshow(m2i(denoiseTV(Y, 0.7)))\n",
    "\n",
    "ax = fig.add_subplot(235)\n",
    "ax.title.set_text('TV B=0.9')\n",
    "ax.imshow(m2i(denoiseTV(Y, 0.9)))\n",
    "\n",
    "ax = fig.add_subplot(236)\n",
    "ax.title.set_text('TV B=1.1')\n",
    "ax.imshow(m2i(denoiseTV(Y, 1.1)))\n",
    "\n",
    "plt.show()\n",
    "plt.rcParams['figure.figsize'] = [16, 10]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}