python/opengm/01_image_denoising.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import opengm\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from PIL import Image\n",
    "\n",
    "# images source: https://www.irisa.fr/vista/Themes/Demos/Debruitage/ImageDenoising.html\n",
    "# input_img = './images/Brain_9pct_u881_z.png'\n",
    "# image = Image.open(input_img).convert(\"L\")\n",
    "# img = np.asarray(image).astype(float)/255\n",
    "# plt.imshow(img, cmap = cm.Greys_r)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate random binary image\n",
    "# source: http://www.shogun-toolbox.org/static/notebook/current/Binary_Denoising.html\n",
    "import scipy.ndimage\n",
    "im_size = np.array([40, 60], np.int32)\n",
    "im_rand = np.random.random_sample(im_size)\n",
    "# print im_rand[0:10, 0:10]\n",
    "im_bw = np.round(scipy.ndimage.gaussian_filter(im_rand, sigma=3))\n",
    "# print im_bw[0:10, 0:10]\n",
    "noise_level = 5\n",
    "t = np.random.random_sample(im_bw.shape)\n",
    "im_noise = im_bw*(1-t**noise_level) + (1-im_bw)*(t**noise_level)\n",
    "im_noise = np.round(im_noise)  # image binarization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAD2CAYAAADF7Q/tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFLVJREFUeJzt3XmsLFldB/DvD0ZFRAEBlVGcUcGFaMS8O5K4RBMXokjA\nRMWFRcUFiIiJBFEjkHFJjFtcILgkiCIK7iAuaBQS3OC9GFGDGzDICAozwrDJfvyj+5n7Ln2H7lvV\nVXX6fj7Jy7x3e6lTp6pP/abu+fap1loAAKBHt5m7AQAAcFaKWQAAuqWYBQCgW4pZAAC6pZgFAKBb\nilkAALqlmOUKVfW0qvqBBbTj46vqrVV127nbAnAIquqPqurhS39P2FX5ntnzpapuSPLRSd6b5N1J\n/irJI1trr5mzXQDcuvX4ffskn9Bae9v6Z9+S5CGttS+csWkwK3dmz6cHtNbukOTuSf47yc/ue4NV\nddW+twFwDtw2yWPnbgQsiWL2HGutvSPJbyW59+WfVdUvV9UPrf/+hVV1Y1V9d1W9vqpeV1XfdOy5\n96+qv6uqN1fVa6rqycceu7aqWlU9oqr+I8mfV9Xzq+oxx9tQVS+rqq882bZjr79q/e8XVtUPVdVf\nracfPK+q7lJVv7be/kur6tpjr//pdZveXFWXqurzjz32oVX1jKp6Y1W9vKoeX1U3Hnv86qr67ap6\nQ1W9qqq+c0g/A4zox5I8rqrutOnBqvqc9Xh4y/q/n3PssReu7+Smqu5ZVS9aP++mqnr2sed9alX9\naVX9T1X9S1V9zWmNOfGe31hVf1lVP1VVb6qqV67b843r8fj1x6ck3No1ZP34w6rq1VV1c1X9QFXd\nUFVfvH7sNlX1hKp6xfrx51TVR56tS+mdYvYcq6rbJ3lwkr+5lad9TJI7JvnYJI9I8pSquvP6sbcl\neViSOyW5f5JHVdWDTrz+C5J8WpL7JXlGkocc2/5nrt/3+Vs2+WuTPHT9mk9K8tdJnp7kI5O8PMmT\njj33pUnus37sWUl+s6put37sSUmuTfKJSb7kRJtuk+R5Sf5+vZ0vSvJdVXW/LdsIsE8Xk7wwyeNO\nPrAu5p6f5GeS3CXJTyZ5flXdZcP7/GCSFyS5c5KPy/o3dFX1YUn+NKtx86OyGnefWlX33vAem9w3\nycvW239Wkt9Icl2Se2Y11v5cVd1h/dxTryHr7T01yTdk9VvEy9ehyx6T5EFZXWOuTvLGJE/Zso0c\nGMXs+fR7VfWmJLdkVcz92K08991Jrm+tvbu19odJ3prkU5KktfbC1to/tNbe11p7WZJfz2pgOe7J\nrbW3tdb+N8lzk3xyVd1r/dhDkzy7tfauLdv99NbaK1prtyT5oySvaK39WWvtPUl+M8lnXX5ia+2Z\nrbWbW2vvaa39RJIPudzuJF+T5Edaa29srd2Y1cB/2XVJ7tZau7619q7W2iuT/GJWAzrAEjwxyWOq\n6m4nfn7/JP/WWvvV9dj360n+OckDNrzHu5Nck+Tq1to7WmsvXv/8K5Lc0Fp7+vo9/i7Jbyf56i3b\n9qr1a9+b5NlJ7pHVNeSdrbUXJHlXVoXtB7qGfFWS57XWXry+RjwxyfGQzyOTfH9r7cbW2juTPDnJ\nV5nSdj4pZs+nB7XW7pTkdkm+I8mLqupjTnnuzeti8bK3J7lDklTVfavqL9a/jr8lq8Hlride///B\nsvW0hmcnecj6DujXJfnVHdr938f+/r8b/n35//ZTVY9bTyG4ZV243/FY264+3q4Tf78mydXrX5G9\naf3a78sqNAcwu9baPyb5gyRPOPHQ1UlefeJnr86VdzQve3ySSvKSqvqnqvrm9c+vSXLfE2PgN2T1\nW7ptnByX01rbOFZ/gGvIFeN0a+3tSW4+9j7XJPndY218eVbBZmP1OaSYPcdaa+9trf1OVgPA553h\nLZ6V1d3We7TW7pjkaVkNjlds5sS/n5HVwPhFSd7eWvvrM2z3Vq3nxz4+qzuwd14X7rcca9vrsvq1\n2mX3OPb312R1Z+FOx/58eGvty8duJ8AAT0ryrbmyUH1tVkXecR+f5D9Pvri19l+ttW9trV2d5Nuz\nmkpwz6zGwBedGAPv0Fp71B724dauIVeM01X1oVlNXbjsNUm+7EQ7b9dae7995fApZs+xWnlgVnOm\nXn6Gt/jwJP/TWntHVX12kq//QC9YF6/vS/IT2e2u7K7tek+SNyS5qqqemOQjjj3+nCTfW1V3rqqP\nzeru9GUvSfKWqvqedVDstlX16VV13Z7aCrCz1tq/Z/WbruMB1T/MairX11fVVVX14KwCvn9w8vVV\n9dVVdblYfGNWNx7et37uJ1fVQ6vqg9Z/rquqT9vDbtzaNeS3kjxgHSD74KymERy/WfK0JD9cVdes\n9+du6+sZ55Bi9nx6XlW9Ncmbk/xwkoe31v7pDO/z6CTXV9VbsprP9JwtX/crST4jyTPPsM1t/EmS\nP07yr1n9iu0duXIqwfVJbkzyqiR/ltWg+c5kdbc6qzlj91k/flOSX8pqmgLAklyf5MMu/6O1dnNW\n49d3Z/Ur+ccn+YrW2k0bXntdkr9dXwuem+SxrbVXttbekuRLs8oJvDbJfyX50axyB2M79RqyviY9\nJqsA2euyymu8PuuxOslPr9v9gvXr/yar8BnnkEUTmFxVPSzJt7XWzjK1YXRV9agkX9taOxleA2AB\n1t+A8KYk92qtvWru9rAs7swyqfXXgT06yS/M2Ia7V9Xnrr+n8FOyuovxu3O1B4D3V1UPqKrbr78u\n7MeT/EOSG+ZtFUukmGUy6+9qfUNWaddnzdiUD07y80nekuTPk/x+Vt9nCMByPDCrqQ6vTXKvrH6D\n5tfJvB/TDAAA6JY7swAAdEsxCwBAtxSzAAB0SzELAEC3FLMAAHRLMQsAQLcUswAAdEsxCwBAtxSz\nAAB0SzELAEC3FLMAAHTrql2eXFVtzI1fuHBhzLdjjy5dujR3E65w6OfO0vr7pB77/4YbbshNN91U\nc7djSne9613btddee8XPhpxb2x7307bR43nTq6WNIYd+7Df196Z9nuu49Nj/u4zZ1dr29enYxewu\n22ZeVcuqAQ793Flaf5/UY/8fHR3l4sWLy+7YkR0dHbWLFy9e8bMh59a2x/20bfR43vRqaWPIoR/7\nTf29aZ/nOi499v8uY7ZpBgAAdEsxCwBAtxSzAAB0a6cAGOfD0uZaAWdz6dKlvX+et50reGjbPmmX\necJjt3vbYzykLa4Luxt7fvrYx2DosV/SPFx3ZgEA6JZiFgCAbilmAQDolmIWAIBuCYDt0ZLCCQD7\ncEhj2pAxe5d+WHqf9RD2cn0dbh/n9lzcmQUAoFuKWQAAuqWYBQCgW4pZAAC6NVkArIcJxFMYslLL\nPvQw0Z/lEb7ow4ULF3Lx4sUrfjbFalG7mCp0ddb3m3PMnmLFp14tbbwZe8WuQzpWm4w9lrgzCwBA\ntxSzAAB0SzELAEC3FLMAAHRrpwDYpjABK3NO1N920vQhTSgXQIKzGfI52UcobK7P7VTj4bb9M1cw\nb0g/nLaNKa6HSxvvxz5+Q0wV3lwSd2YBAOiWYhYAgG4pZgEA6JZiFgCAbk22Ahj7c0jBLlamCAkc\nEn0znzlXK1xSgGWXc3CKlaGm+EzsY5+XdEx3se2YvbTjMpexr3HuzAIA0C3FLAAA3VLMAgDQLcUs\nAADdEgCDTiw9FLaPgM6S9u+QjX3shh63TdvuNRg0xNI/80Nte0wd+5WxVxk7pDHbnVkAALqlmAUA\noFuKWQAAuqWYBQCgW4pZAAC65dsMOBhLX+5yH5aedl5SWzjd2J+TXd5v6efI0M/YkNfP1TdT7XOP\nY/Zp/TD2N3BMseTx2K+dkzuzAAB0SzELAEC3FLMAAHRLMQsAQLcEwM65uSbb9zrJHBhXr4GfuSy9\nv5Ll9dmY9tH/Y/fX0oPB++DOLAAA3VLMAgDQLcUsAADdUswCANAtATBmMdUKQT2uMAOwi23HuSHj\noTG7H4ce9trEnVkAALqlmAUAoFuKWQAAuqWYBQCgWwJg55zJ9n07jxP92d6lS5fe7xzp9fM9ZKya\nc4Wlsbc9VYhrCkOO3xT7ssuxm+tcZMWdWQAAuqWYBQCgW4pZAAC6pZgFAKBbOwXANoUJ5rS0yezs\nxy5hi0Nx6Ps3tm3GgqOjowlasnxznltDxuyxXztk1ayp+nDI2Dd2YG6X/h+7H8e+1m+7f1OtejbX\nNW7Oa+vYY7Y7swAAdEsxCwBAtxSzAAB0SzELAEC3ul4BbOwJ16xYFWw/hLrG4Vycxj7CIVOEZIas\nhrW0cNxc/TX0MzZXP4597Rq6H0sLFJ51u6f14VwBvk3cmQUAoFuKWQAAuqWYBQCgW4pZAAC61XUA\nbIixV0c5dEsLhS09BHQez5F9WPpxPhQ9hKaGrNi1yRQriu3y3CH9OGSf9zG2z/W5Xdp5PPZ7TtHu\npQUjt+XOLAAA3VLMAgDQLcUsAADdUswCANCtcxsA26SHSc5LsrRQ2BScI/tz6OfOki1pJZ9dDAmF\nTRG4GrqdbQ3Z5yHB56WdD5tMNWZPEZyaYl922caSjr87swAAdEsxCwBAtxSzAAB0SzELAEC3BMAA\n2EoPAcghbdzHimJDQmpThNm2tUvYZ+x9nssu58OSgn69GrIv7swCANAtxSwAAN1SzAIA0C3FLAAA\n3RIA48yWNFGfZXKO0JOlrQo2drhnaat9LX182Ee4auzw39jmPCYnt310dLT1a92ZBQCgW4pZAAC6\npZgFAKBbilkAALolAPYBDJkM3evKHEuflM+tO4/BAZZll9Wwhr7nNtsY8toejP3ZG/p+Y7dn6Cps\nYxp6bgt77Yc7swAAdEsxCwBAtxSzAAB0SzELAEC3BMCOmWsS/VShg0Oa7M3KtsfUseeyKcIm+xjT\ntg3ODFlhyZi9TFNcm+dale00c4Xoplr9bWzuzAIA0C3FLAAA3VLMAgDQLcUsAADdEgBbgB4mV7My\n10otMJbzeL4OCYVt+35TmSKg02sIaJOljdlz9WMP5/YQ7swCANAtxSwAAN1SzAIA0C3FLAAA3RIA\ng070OjF/Cucx1LSNCxcu5OLFi1f8TF+t9Pp5GjvMtu02DslUobAh/ThX0G+TfZwPY++fO7MAAHRL\nMQsAQLcUswAAdEsxCwBAtxSzAAB0q3ZJj1XVQUccDz3ByX4sLQV76Ib0d2vtXEX5N43ZS1vec5Mp\nEvtzmmtftt3unCn3bQ3pw7H74TS9np8nndYPQ/p7W9uO2e7MAgDQLcUsAADdUswCANAtxSwAAN0S\nADuDQ5nUze6mCso4x1bG7m8BsOWZIpDm8zRcrwG8pZ1LvfbjJlNcDwXAAAA4eIpZAAC6pZgFAKBb\nilkAALp11Zwb38dE6ikc0gRuWIqlrULFdnYMEW/1syEOffWwKczZXz2sPjZEr+0e29jBT3dmAQDo\nlmIWAIBuKWYBAOiWYhYAgG5NFgDbx6TnKVaO2dahT1qHsQh69WGqMXtbU5w3Q4JPp7VvrkDakG3M\neZ1a0nV9k13acijX+330/9jv6c4sAADdUswCANAtxSwAAN1SzAIA0K1ZVwA7j6xE0485QweHcp4s\nKbjBMo19jmz7fkOet+mzuMvnc+wg3ND2nHUbUxkawlvKNoYa+7gM3ZchYT0rgAEAwJpiFgCAbilm\nAQDolmIWAIBu7SUAtrSJ4kIoMC6fqcMyVYhk26DStkGXucImU5kriLWPbQw5pnPVFFOdI3OF+rY1\n1fkwhDuzAAB0SzELAEC3FLMAAHRLMQsAQLesAAan6DU0ArsaO3QzNDAydihsis/tadsYux+XFJDa\nxbZtPPSw1yZLP349XPfcmQUAoFuKWQAAuqWYBQCgW4pZAAC6dS4CYEsK8ix9ojfL4DxhDBcuXMjF\nixev+NmQANHSwkdTjONDg0tLC9cdsjmv9WOvzjX2cV5SHbQP7swCANAtxSwAAN1SzAIA0C3FLAAA\n3ToXAbC5mKgP4zvr5+ro6GjklizfpUuXRg15HPqYto8Vt5beZ/sIH80VMuw10DTFSm899M3J/dtl\nzHZnFgCAbilmAQDolmIWAIBuKWYBAOjWXgJgS1slZgqHvn9w3FSryfhcHZYlXRuGruw1trn6Zh/b\nmGuFsx6CYtse57EDc0PG7B7GYXdmAQDolmIWAIBuKWYBAOiWYhYAgG5NtgLYkib+z73tsQ2ZzH5I\n/TCFsYNPh9T/U4XCGN9U4/NcwaAhhvbNkM/A0q6bJ522b0tq4yZTrEa2y7a33c7Y7dl2G7ucw3Ot\nUubOLAAA3VLMAgDQLcUsAADdUswCANCtyQJgmyx9cvsSjR2oOaSAztJXEjqP9E2/ehifp1jBqId9\n7iF8uaTzaexVuKYyxSpeSwrmHx0dbf1ad2YBAOiWYhYAgG4pZgEA6JZiFgCAbs0aANtkHytNLG0i\nPMCS7GMloG23s8lcQdderylztee0/tq2PUsKhR2SsVfh6uGYuDMLAEC3FLMAAHRLMQsAQLcUswAA\ndEsxCwBAtxb3bQa7WFqilHn1msJkOGPBMFP135BvFZiijXOeRz0sSXsolnZdGHu55SHXwl3Ow7G/\nNWEId2YBAOiWYhYAgG4pZgEA6JZiFgCAbnUdAIPjljapn+kIzwwzJESyj34+j8eux33eJRg05Hya\na2xfUltOM0W/Dt3nKcYSd2YBAOiWYhYAgG4pZgEA6JZiFgCAbgmAAQfpZMDg6OhoppYs37bBi6Eh\npSGhj6WH/E4LyYzdxqWF9TYZeztDwktzBp+mMOe+bNu3U/SjO7MAAHRLMQsAQLcUswAAdEsxCwBA\ntwTAADizXYJGQ4JBcwWplmaqsN5cxg4LLS3ENSTweOiGhHbdmQUAoFuKWQAAuqWYBQCgW4pZAAC6\nJQBGl87j5HgYw9jBp30EjaZYKWyKVap6NfQcmWJ8nmtVsKHmunbNuc/bGvIZcmcWAIBuKWYBAOiW\nYhYAgG4pZgEA6JYAGMA50kNQaYpA2hQhs329/qyGBKR2eb8pwkZD3m/O4NO2fTN2H+5jn8cO3FkB\nDACAc0kxCwBAtxSzAAB0SzELAEC3dg2A3ZTk1ftoCOyihxALi3PN3A2YgTH7FHOuejbX+DXnqmfG\n7NNt2zeH1Idb7svWY3YtbTkzAADYlmkGAAB0SzELAEC3FLMAAHRLMQsAQLcUswAAdEsxCwBAtxSz\nAAB0SzELAEC3FLMAAHTr/wDUpDZhEXUvlAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f57368382d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# black is 0, white is 1\n",
    "fig, plots = plt.subplots(1, 2, figsize=(12, 4))\n",
    "plots[0].matshow(im_bw, cmap=cm.Greys_r)\n",
    "plots[0].set_title('Binary image')\n",
    "plots[1].matshow(im_noise, cmap=cm.Greys_r)\n",
    "plots[1].set_title('Noise image')\n",
    "# print im_bw[0:10, 30:40]\n",
    "for p in plots:\n",
    "    p.set_xticks(())\n",
    "    p.set_yticks(())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dimy, dimx = im_noise.shape\n",
    "nvar = dimx * dimy\n",
    "num_labels = 2\n",
    "nstates = np.ones(nvar, dtype=opengm.index_type) * num_labels\n",
    "gm = opengm.graphicalModel(nstates)\n",
    "\n",
    "beta = 0.3\n",
    "nu = 0.3\n",
    "\n",
    "# add unary factors\n",
    "for y in xrange(dimy):\n",
    "    for x in xrange(dimx):\n",
    "        f = np.ones(2, dtype=np.float32)\n",
    "        if im_noise[y,x] == 0:\n",
    "            f[0] = -nu\n",
    "            f[1] = nu\n",
    "        else:\n",
    "            f[1] = -nu\n",
    "            f[0] = nu\n",
    "        fid = gm.addFunction(f)\n",
    "        gm.addFactor(fid, (y * dimx + x,))\n",
    "\n",
    "# x_i and x_j factor\n",
    "f1=np.ones((2,2), dtype=np.float32)\n",
    "f1[0,0] = -beta\n",
    "f1[1,1] = -beta\n",
    "f1[0,1] = beta\n",
    "f1[1,0] = beta\n",
    "\n",
    "\n",
    "fid1 = gm.addFunction(f1)\n",
    "\n",
    "for y in xrange(dimy):\n",
    "    for x in xrange(dimx):\n",
    "        if x + 1 < dimx:\n",
    "            gm.addFactor(fid1, np.array([y * dimx + x, y * dimx + x + 1]))\n",
    "        if y + 1 < dimy:\n",
    "            gm.addFactor(fid1, np.array([y * dimx + x, (y+1) * dimx + x]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# opengm.visualizeGm(gm,plotFunctions=False, layout='neato', plotUnaries=True,\n",
    "#                     plotNonShared=True,relNodeSize=1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beginning of inference\n",
      "energy   -1362.60005414\n",
      "energy   -1672.80006647\n",
      "energy   -1743.00006926\n",
      "energy   -1750.80006957\n",
      "energy   -1753.80006969\n",
      "energy   -1750.80006957\n",
      "energy   -1749.0000695\n",
      "energy   -1753.20006967\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "energy   -1755.00006974\n",
      "end of inference\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "opengm.inference._inference.InferenceTermination.NORMAL"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "imgplot=[]\n",
    "\n",
    "class PyCallback(object):\n",
    "    def appendLabelVector(self,labelVector):\n",
    "        #save the labels at each iteration, to examine later.\n",
    "        labelVector=labelVector.reshape(self.shape)\n",
    "        imgplot.append([labelVector])\n",
    "    def __init__(self,shape,numLabels):\n",
    "        self.shape=shape\n",
    "        self.numLabels=numLabels\n",
    "        #matplotlib.interactive(True)\n",
    "        plt.interactive(True)\n",
    "    def checkEnergy(self,inference):\n",
    "        gm=inference.gm()\n",
    "        #the arg method returns the (class) labeling at each pixel.\n",
    "        labelVector=inference.arg()\n",
    "        #evaluate the energy of the graph given the current labeling.\n",
    "        print \"energy  \",gm.evaluate(labelVector)\n",
    "        self.appendLabelVector(labelVector)\n",
    "    def begin(self,inference):\n",
    "        print \"beginning of inference\"\n",
    "        self.checkEnergy(inference)\n",
    "    def end(self,inference):\n",
    "        print \"end of inference\"\n",
    "    def visit(self,inference):\n",
    "        self.checkEnergy(inference)\n",
    "\n",
    "inf=opengm.inference.BeliefPropagation(gm,parameter=opengm.InfParam(damping=0.00))\n",
    "callback=PyCallback(im_noise.shape,num_labels)\n",
    "visitor=inf.pythonVisitor(callback,visitNth=1)\n",
    "\n",
    "inf.infer(visitor)\n",
    "# inf.infer()\n",
    "# result = inf.arg()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAFhCAYAAAAr/EDTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X20JVV55/HfAzQ2NNC0CAqiMIKoRBbO9AvjSyKMJo62\nLJg14isCviS+rJA4IxFGCbIYND1LE2N8GWKSAdsWFYw6IjGjLsVZRBngLhb2GJgkyFtLK9I2DbS2\nNrDnj6p72Gdza999au86p865389avbrurbddVafOvvU8tfc255wAAEB7e0y6AAAATDsqUwAAMlGZ\nAgCQicoUAIBMVKYAAGSiMgUAIBOVaQFmdomZ/XEPyvF0M3vIzPacdFmAWWdmXzezM/u+TYyH0c50\ncWZ2h6QnS3pE0m5J35P0dufc3ZMsF4D26vt6X0n/yjm3s/7dWyWd7pw7cYJFwxTiyTTdyc65/SQd\nKumnkj7W9Q7NbK+u9wEscXtK+sNJFwLTj8p0RM65XZK+KOnY+d+Z2WVmdnE9faKZbTGzd5vZvWa2\n1cze5C273sxuMrMHzOxuM7vQm3ekmTkze4uZ3SXp22Z2tZmd7ZfBzH5gZv8hLJu3/l71z9eY2cVm\n9r06/HuVmR1kZp+t93+DmR3prf/RukwPmNmcmf2mN28fM/u0mW03s1vM7D1mtsWbf5iZ/a2Z/czM\nbjezP8g5z8CYfEjSOWZ24EIzzewF9X2yo/7/Bd68a+onWZnZ0Wb23Xq5+8zsC95yzzazb5rZz83s\n/5nZq5sKE2zzLDP7BzP7iJndb2Y/qstzVn2f3uuHhGPfLfX8M8zsTjPbZmZ/bGZ3mNlL63l7mNl5\nZnZbPf8KM3tiu1O6NFGZjsjM9pX0GknXRRZ7iqSVkp4q6S2SPmFmq+p5OyWdIelASeslvcPMTg3W\nf7Gk50h6maRPSzrd2//x9XavTizyayW9sV7nKEnfl3SppCdKukXS+71lb5D0vHre5ZKuNLPl9bz3\nSzpS0jMk/XZQpj0kXSXp5no/L5H0LjN7WWIZgUm5UdI1ks4JZ9SVydWS/kLSQZL+TNLVZnbQAtv5\nr5K+IWmVpMNVR67MbIWkb6q6nw5RdT9+0syOXWAbCzlB0g/q/V8u6fOS1ko6WtU9+HEz269etvG7\npd7fJyW9QVV0bf77ad7Zkk5V9d1zmKTtkj6RWEaIynQUXzGz+yXtUFWZfCiy7G5JFznndjvn/k7S\nQ5KeJUnOuWucc5udc486534g6XOqPsC+C51zO51zv5T0VUnHmNkz63lvlPQF59yvE8t9qXPuNufc\nDklfl3Sbc+5bzrmHJV0p6V/PL+ic2+Sc2+ace9g596eSnjBfbkmvlvRB59x259wWVV8w89ZKOtg5\nd5Fz7tfOuR9J+itVXxxA310g6WwzOzj4/XpJ/+yc+0x9T3xO0q2STl5gG7slHSHpMOfcLufctfXv\nXynpDufcpfU2bpL0t5JOSyzb7fW6j0j6gqSnqfpu+ZVz7huSfq2qYl3su+VVkq5yzl1bf3dcIMl/\nYebtkt7nnNvinPuVpAslvYpUUzoq03SnOucOlLRc0u9L+q6ZPaVh2W11ZTXvF5L2kyQzO8HMvlOH\nQ3eo+hA/KVh/8GJTHVb+gqTT6yfA10n6zAjl/qk3/csFfp7/q1Zmdk4dwt1R/+Gw0ivbYX65gukj\nJB1Wh6Lur9d9r6qXtoBec879X0lfk3ReMOswSXcGv7tTw090894jySRdb2Y/NLM3178/QtIJwb3x\nBlXRqxTh/Srn3IL38CLfLUP3r3PuF5K2eds5QtKXvTLeouqFS+7hRFSmI3LOPeKc+5KqD9qLWmzi\nclVPm09zzq2UdImqm3BoN8HPn1Z1A75E0i+cc99vsd+oOj/6HlVPoKvqPxx2eGXbqip8Ne9p3vTd\nqv6CPtD7t79z7hWlywl05P2SflfDFeU9qioZ39Ml/Thc2Tn3E+fc7zrnDpP0NlWh3KNV3RvfDe6N\n/Zxz7+jgGGLfLUP3r5ntoyp0PO9uSS8PyrncOfe4Y8XCqExHZJVTVOVGbmmxif0l/dw5t8vM1kl6\n/WIr1JXno5L+VKM9lY5arocl/UzSXmZ2gaQDvPlXSPovZrbKzJ6q6ul83vWSHjSzc+sXlfY0s+ea\n2dqOygoU5Zz7F1URIP/Fub9TlWJ5vZntZWavUfXi4dfC9c3sNDObr6y2q/qD+NF62WPM7I1mtqz+\nt9bMntPBYcS+W74o6eT6Baa9VYVx/T/iL5H0ATM7oj6eg+vvOSSiMk13lZk9JOkBSR+QdKZz7oct\ntvNOSReZ2YOq8hZXJK63UdJxkja12GeK/yXp7yX9k6pQ1i4Nh3IvkrRF0u2SvqXq5vyVVD2tq8oN\nPa+ef5+kv1YVJgamxUWSVsz/4Jzbpupz/W5VIdH3SHqlc+6+BdZdK+n/1N8RX5X0h865HznnHpT0\nO6reH7hH0k8k/TdV7yOU1vjdUn9Xna3qBaatqt7juFf1PSzpo3W5v1Gvf52ql5+QiE4bpoSZnSHp\n95xzbULLxZnZOyS91jkXvjwFoOfqN4Dvl/RM59ztky7PLODJdArUzXHeKelTEyzDoWb2wro92rNU\n/bX+5UmVB8BozOxkM9u3bq7zYUmbJd0x2VLNDirTnqvbav5M1Vt9l0+wKHtL+ktJD0r6tqT/qard\nGoDpcIqqUPM9kp6pKrJEaLIQwrwAAGTiyRQAgExUpgAAZKIyBQAgE5UpAACZqEwBAMhEZQoAQCYq\nUwAAMlGZAgCQaaSBX82sVQ8Pq1evbrPakjA3N9fZtvt43rs83oVM6hzMzc3d55wLB5vujSc96Unu\nyCOPlDTaNWk6n+E2+vjZ69pSv5f9MnZxLvp+L4/UA1LbypRelpqZhUOZltPH897l8S5kUufAzOac\nc2smsvMEa9ascTfeeKOk0a5J0/kMt9HHz17Xlvq97Jexi3PR93t5pCdT5Bt3ZQIsZG5uLvuz6K9f\n4ouu6wo5Vt42xxI7f7GKpetKZ1JSj6XE8Ze+liWQMwUAIBOVKQAAmahMAQDItGRzppOKqwOzoi/3\nTeq93HZeaePMk87q91xfrqWPJ1MAADJRmQIAkKmTytQ5N/Sv78ys8V/p7WO8OPcLW716ddL9Gd7L\nbc5n6jptvzeaylfiXm67vb587rr8Hm57vfpybnwlPis8mQIAkInKFACATFSmAABkGqlpzOrVqzXf\nn+e0KdFtVSi1n9JJob9UpBg1PzmvRJOUNkrcX227+EvtJjC2r5Rtj7JebJ0u7/lRvhvbHEtfuglM\nxZMpAACZqEwBAMi0ZHtAKqEv4dxpUiL80xfTXPYSSo8U03XoLna9Uq9liW2U3vakQp6xe7mLUHyX\nSnwv8WQKAEAmKlMAADIR5sVEjXOg5NIDQM+qLs9TuL0+vpXp62taoum89f18Su3enA6X7eO9zJMp\nAACZqEwBAMhEZQoAQCZypktI33sQGWd+qi+5rz5q89mIrdPHJg6py3VR9hK9LU3qXo7tt3TusvRy\nXePJFACATFSmAABkIsw7Bn3s6BsYl0mFIbswS8fSRtvjL93RfR/xZAoAQCYqUwAAMlGZAgCQiZzp\nFCrRDKHvzWSAvksdYLztwOlL/V7ue440xJMpAACZqEwBAMhEmHcMZjUMU9q0hXWm2dzc3OB8T/Iz\nmXpvtP1sNIVe226vRPi2hKbtx0LPbaX2IpWyzizjyRQAgExUpgAAZBopzOuHhrpACDTfOAfbbqOP\nZWprVsJcXXfmXnq51LdoF/p5VKP0wtMmZJ0aKk49rhI9FKV2YD/KuS39vdTl9kKp2+fJFACATFSm\nAABkojIFACBTr5rGpMampzm32sWr6300TTnDUczq9WrSduSONvdy6r3R9WgibfNxqT0glR5EO1Xb\n754S17LrY0ktR9NyJfBkCgBAJipTAAAy9SrMm6rrXlPGqcvekZb6QMYlLLWwrtS+uUYbsc9/ieYf\nqaHXNr0jtQ01lm42lKrr6zqpAcD7Mog4T6YAAGSiMgUAIBOVKQAAmaYyZ+qbpVzdNIwuM0vnu0lf\nz/24jLM5QWy/JUZ8aVq2iy7omspYomvBLoyzq8WuyrDYNsZ5TnkyBQAgE5UpAACZpj7MC2B8xhnm\nL9HzTupyXTahiWkbAh5neLj0AOtN2y61zRLalIMnUwAAMlGZAgCQiTBvjyz1t0jHiXM9O0r0WFT6\nrd+u3+Dte+9mqWHprkPlJbbJ4OAAAIwJlSkAAJmoTAEAyDRzOdNJDbybilxdumnIraDStolDidzi\ntN3LXY4a03VvQKn3ZGrToxLl6AueTAEAyERlCgBApqkP83YRaindATTSpA5ejLIm2WF5LJSb2isR\n9/JjSnwftmlC1HZfo0htvjSpa8aTKQAAmahMAQDIRGUKAECmqc+ZdmHa8iTjVPoVd0zeNFzH1BFf\nYuuVUKIrwEkOCO7r8l7uuou/rvfVBk+mAABkojIFACATlSkmyjk3+DdtzKzxX9+tXr16as+7/5kZ\nd/n9/ba95n38zJcoU+o2Stwnpe+9EtugMgUAIBOVKQAAmXibF2PVp9BWG9MQwk0xNzc3OJYSPSCl\nGue+RtHmDdu2nfQ3HXMX90abwQK6vkalO9wPlRggoc0x82QKAEAmKlMAADJRmQIAkGnqc6ZdD4aL\nYX3JcY3LUjjero+xxIgkXd/Lpbef2ivTOM99KHUUltTtxbYxzu/iSfXSxpMpAACZqEwBAMg0ljBv\n29BArr50KI3pshRCuyWkhhDbns/UwcGnTelj6WOqa9quUYnmQDyZAgCQicoUAIBMVKYAAGTqJGfa\nNl7e5SvNfcwroD/Iky6sxL0cU+K8l+4yr+tBv7tshtP153gavkfbnIMS540nUwAAMlGZAgCQaep7\nQCphkr119F3X4fZJnWvCuuOTM2DzqNuLfb5KhF5jodwS28jVdpSU0ttI1fb7oM2xjBKipmkMAAAT\nQGUKAECmYmHeaX5DDbOLz87oSoTamsKrsbBe6vbGeU27SEu02Ubqeet6gPES20gNvbctR4nz2wZP\npgAAZKIyBQAgE5UpAACZaBqDqBKvjGO6tMnBtV0utVeicfaU5Ouip6Q2xtmjUhf3+DjP1aS+o3gy\nBQAgE5UpAACZpqIydc4N/Zu27eMxnOt+WL16deM1iF0fMxv8K30d/W3nDCjeVC7/9233Nc2f3dL3\n3ijby72u4f76aCoqUwAA+ozKFACATFSmAABkWrJNY/oad0e/jXPQ6y7Nzc21KmMf75vU7gpTft+F\nEiOjtB30PFaOcSrR/eE4y9+m+0qeTAEAyERlCgBApmJh3kn1DJKqj2VC93JHG+Fz086kBsAusa/S\nZW+7ja57n/KNc9DvtiHrtgN9jwtPpgAAZKIyBQAgUydv83YxoG6X22ury8Fq+6pN2LQvx8+A82lK\nhDlLhyhTpX73pF7/rr/LJrWvUIm3hWPba9pG22OObaPp2nb95jBPpgAAZKIyBQAgE5UpAACZxtID\nUt+bzaRqm2frS35uVnPXbUxz2cdpkjnDNvvt+rqmDmbe5X4X2nfJ7ceapJQ2yjkcZ29W9IAEAMAE\nUJkCAJBp7B3dt3mM70uYFFgK2jQ7iG2jzfqjrDepzt1Lb3uUkOespM5i2g4qQQ9IAABMKSpTAAAy\nUZkCAJCpV4ODkxvt1lLIs4zTrH5eSxxXapduXe6rraXY7WTX3wcluglMbRrT5nuuxHXmyRQAgExU\npgAAZOpVmBfdIrRb1jh7wBmn1FBbm+OdhnM0qTKOcq6b5nVxj48zPZTaw1SJAeJL92bFkykAAJmo\nTAEAyESYFyhkVt4CTe15Jyb1XIzznJXYV5ch8FCJMHqJQb/7Mkh5m3KMswN/nkwBAMhEZQoAQCYq\nUwAAMpEzBdBKidFl2uQFJ9kkqUQ+ubS2ub++D/Q+SQwODgDABFCZAgCQiTDvDJuGcAomr23YtOsQ\nbVOorYsBxsepTbOhtlKbv5Tu5WiWemJKxZMpAACZqEwBAMhEZQoAQCZypsAS15emJanLlsi7tl0u\nVWp+MrZOiRxhm0G5S+hilJs2ZWybJ6ZpDAAAE0BlCgBAplHDvPdJurOLgqC8vjQFWKKOmHQBFjG1\n93KJUHHONrvaXh97Vyqlqfx9Oa5FypF0LxttEQEAyEOYFwCATFSmAABkojIFACATlSkAAJmoTAEA\nyERlCgBAJipTAAAyUZkCAJCJyhQAgExUpgAAZKIyBQAgE5UpAACZqEwBAMhEZQoAQCYqUwAAMlGZ\nAgCQicoUAIBMVKYAAGSiMgUAIBOV6QSZ2RvM7BsdbNeZ2dEN875uZmeW3icwDczsMjO7OHMbJ5rZ\nllJl6pKZvdfM/rrQthq/VyDtNekCTJqZ3SHpyZIelvSIpH+UtFHSp5xzj3a5b+fcZyV9tst9LLDP\nl49zf8A4effzI5J2S/qepLc75+6eZLkmxTn3wUmXYangybRysnNuf0lHSNog6VxJfzPZIgFo6WTn\n3H6SDpX0U0kfm3B5OmFmS/5hqE+oTD3OuR3Oua9Keo2kM83suZJkZk8wsw+b2V1m9lMzu8TM9qnn\nnWhmW8zs3WZ2r5ltNbM3zW/TzFaa2UYz+5mZ3Wlm55vZHvW8s8zs2nrazOwj9TYeMLPNKfuv5/9R\nvd97zOzNsWM0s2vM7K3e/v+h3u/9ZvYjM3tB/fu767Kc6a273sxuqst3t5ldGGz7jPoYt5nZH5vZ\nHWb20nreHmZ2npndVs+/wsyemHG5gCjn3C5JX5R07ELzzWyVmX2tvje319OHe/OfaGaX1vfVdjP7\nSsN2/sDM/tFf15tX7B4zsyPrUOtbzOwuSd+ufx+77y40s03B+mfW3yX3mdn7vO2vM7Pv1+XcamYf\nN7O9G475FfUxP2hmPzazc2LXYimgMl2Ac+56SVsk/Wb9qw2SjpH0PElHS3qqpAu8VZ4iaWX9+7dI\n+oSZrarnfaye9wxJL5Z0hqQ36fF+R9Jv1ftZKenVkrYttn8z+/eSzpH025KeKemlIx7uCZJ+IOkg\nSZdL+ryktfV+Tpf0cTPbr152Z13+AyWtl/QOMzu1Lsexkj4p6Q2qngjmz8e8syWdWp+DwyRtl/SJ\nEcsKJDOzfVX9YXxdwyJ7SLpUVUTq6ZJ+Kenj3vzPSNpX0m9IOkTSRxbYxwWSzpL0YudcUx61yD3m\nebGk50h6WcJ9t5AXSXqWpJdIusDMnlP//hFJ/0nSkyQ9v57/zoZt/I2kt9URveeqrtiXNOfckv4n\n6Q5JL13g99dJep8kU/UBP8qb93xJt9fTJ6q6Cffy5t8r6d9K2lPSryUd6817m6Rr6umzJF1bT/87\nSf9Ur7eHt/xi+/8fkjZ4846R5CQd3XC810h6q7f/f/bmHVev+2Tvd9skPa9hW38u6SP19AWSPufN\n27c+9pfWP98i6SXe/ENV5bT2Wmjb/ONfm3/1/fyQpPvrz9c9ko7z5l8m6eKGdZ8naXs9faikRyWt\nWmC5EyX9WNKfSbpW0spIeUreY0fW6z7Dm7/YfXehpE3B+od7y18v6bUN+36XpC97Pw++VyTdpeq7\n7IBJX/O+/OPJtNlTJf1c0sGqPqBzdfjjfkl/X/9+3jbn3MPez7+QtJ+qv/CWSbrTm3enFvjL0Tn3\nbVV/FX9C0r1m9ikzOyBh/4dJ8l+u8PeV4qfe9C/rsoS/20+SzOwEM/tOHRbbIent9TE+rhzOuV/o\nsSdrqfrr/8veMdyi6i/hJ49YXmAxpzrnDpS0XNLvS/qumT0lXMjM9jWzv6xDpA9I+t+SDjSzPSU9\nTdLPnXPbG/ZxoKTfk/Qnzrkdi5Sn1D02z7/fF7vvFvITb3r+u0pmdkwd6v5JfT4+uMC+5/1HSa+Q\ndKeZfdfMnr/IPmcelekCzGytqgrvWkn3qfqw/4Zz7sD630pXveCwmPtU/XV8hPe7p6v6q/ZxnHN/\n4ZxbrSrHc4ykP0rY/1ZVN76//a5cLumrkp7mnFsp6RJVT87z5fDzTfuoCmvNu1vSy71jONA5t9w5\nt+C5AHI55x5xzn1J1R9tL1pgkXerCnee4Jw7QFWaRao+03dLeqKZHdiw+e2SXinpUjN7YcFix+6x\nec6bXuy+G8V/l3SrpGfW5+O9C+y7KoBzNzjnTlEV/v6KpCta7nNmUJl6zOwAM3ulqpzGJufcZlc1\nj/krSR8xs0Pq5Z5qZi9bbHvOuUdUfcg+YGb7m9kRkv6zpE0L7Htt/VfpMlVh3V2SHk3Y/xWSzjKz\nY+sc0fvzzkLU/qr+Wt9lZuskvd6b90VJJ9cvV+ytKrzk34iXqDoPR9THcLCZndJhWbHEWeUUSatU\nRUJC+6v6Q/V+q16GG9w7zrmtkr4u6ZNWvai0zMx+y1/ZOXeNqlzll+r7oYTYPbaQxe67Uff9gKSH\nzOzZkt6x0EJmtrdVbeRXOud21+t02oxwGlCZVq4yswdV/TX6PlW5EP8loXMl/Yuk6+rwx7dU/UWb\n4mxVleOPVD3pXq4qzxk6QFWluV1VqHabpA8ttn/n3NdV5VW+XS/T5YsA75R0UX2uLpD316hz7oeq\njvXzqv5afkhV7vhX9SIfVfUX9zfq9a9T9WIGUNpVZvaQqi/5D0g6s/58hv5c0j6qoj/XqUqf+N6o\nKrJ0q6rP8rvCDTjnvinpzfU+/02BsjfeYwtJuO9GcY6qyvtBVd9FX4gs+0ZJd9TfR29X9UfFkmZ1\nMhkoqn478X5VIaPbJ10eYCngvpscnkxRjJmdXL/UsULShyVtVvV2JYCOcN/1A5UpSjpFVVOEe1S1\neX2tI/QBdI37rgcI8wIAkIknUwAAMlGZAgCQicoUAIBMVKYAAGSiMgUAIBOVKQAAmahMAQDIRGUK\nAECmvUZZ2Mxa9fCwevXqNqvNpLm5ubHtq4/nfZzHL03uHMzNzd3nnDt48SUnY5999nH777+/JGnH\njuHhOB9++LGheVesWDE078gjjxxML1u2bDC9e/fuoeX8ebPEP87wvG3f/tjQpzt37szel3/u/fMu\nTe78+sd/xx13NC5X+vil5s9eF/zOjG666aake3mkHpDaVqb0svQYs7ajI42uj+d9nMcvTe4cmNmc\nc27NRHae4JBDDnGvetWrJElXX3310Lxt2x4bW3rduuGRxS677LLB9KGHHjqY3rp169By/rxZ4h9n\neN6uvPLKwfT111+fvS//3PvnXZrc+fWP/6yzzmpcrvTxS82fvS7s2rVrMH3AAQck3csjPZkCmA07\nduwYVAZ33XVX8nqrVq0aTPtfOE9/ev6Y9A8++ODQz/NPzqX42w+3HZvXZPPmzUM/+xXIQQc9Nj53\neH79cxU797fddttg2j/v0uSeTMNy+L7zne8kbaPN8Yf79j97Ja5lCeRMAQDIRGUKAEAmKlMAADKR\nMwWWoCc84Qk66qijFpzn//7Zz3524zbGmY+K8fNusZze8uXLG+f5x+K/sRq+WLVly5bB9K233jo0\nz3+D1c+ZxsrbtS7zh+Fnw88Zl3ibN+S/LR3L0U/qc8mTKQAAmahMAQDI1EmYt4/tG2NibR/bHMu4\n21ICo1q1apVOO+20RZc7/vjjh372Q22p4bTUUGPq9sIwqR96Ddt++vxjOfzwwxuX87e3adOmoXl+\naDfWFKR0KNc/71L6ueoy5HnccccN/eyHtmNh3rbnxr8ufjvTEs2EwmZZsZRAE55MAQDIRGUKAEAm\nKlMAADKNlDNdvXq1brzxxq7K0qm2eUzyn9PNv37Tlsvv0sqVK7V+/XpJozUn8fNTpXOhMf6+NmzY\nMDTPz2OGXdD5zXzC7v+a+NsL+5iNNXlpaq6RmiOMNffw84Wx9brukjFVeCxt8qRh0y0/f+3nvEt0\nZRnmpNvgyRQAgExUpgAAZKIHpBlWuslPCeF+ZymMPk3HsmzZskHzgrZNC0qHb1O3F/Y8lNpEJXVY\nsFizji6afDSt74/nGTbR8ZulnH766Vn7HYUf9g+bTflh2dQRZGJi2/DD3mGaoqlZiz/SjDQc2g3D\n6DfffHNyOefxZAoAQCYqUwAAMhHmxZLR5s3eaQrdltIm9Bq+Rdok3F6bUHH4xm6qNp2vl3grNeSH\nb1NDymFoO+x9aF7Xb+/6KYGwF6nYoAipA6f7wnPvb8MPe6eGuWO9WcXeAk/FkykAAJmoTAEAyERl\nCgBAJnKmmKhY7rLLfOVSzIWmapN3i62Tmk8tITXHmbpc+Ps2ve2E+Tc/t+iPchPbV9isJ5af7JJ/\nLcNzcd555w2mw6Ymfr5y48aNSfuK5VP9bcRGCvJt27Zt6OdYvjpcNgVPpgAAZKIyBQAgE2HeCWPw\nccy6Es01UjsiLz0odyyUGWoz+LgvbP7S1Jwk1GbA9rZi2/cH7A6P0T+21CZKsVC8v402TZ66wJMp\nAACZqEwBAMhEZQoAQCZyplNoUs1JMDt2796trVu3ShrOdUntR5Epbffu3YPp2ODYfeHnSf2RTMI8\no3++/XVizUli/PUmeS3DUVl8bbqALJ3/7hpPpgAAZKIyBQAg00hh3rm5uUEYcVKDS4+ijyHPSZ63\nvl+zPl6vWbVjx45BzzHhIM9N4cpR+AM0x0KNsRFq5sPQUrvBmqXh5hV+2DA1hBj2hLN58+bBdDha\niR9ijR2zPy/Wo1LTyDDS8DXzlwtDrSXCvP41ijVR8sPNYYi6TY9CXfNH7wmbHrUZ6JwnUwAAMlGZ\nAgCQqfXbvF2E5PoehpwGfT+HsxTKnea3qrdv364rr7xSkrRhw4ZW2/BDYWHH634INDVsHIZem97g\njQ3cHG4j943QMPwX9lLk88PSsfBtU2g7XGf9+vWD6djbzH7otU1H/NLweQrfCPa3Hwvl+ucm7Jjf\nP48leixqCt/H+GFdSTrjjDMG07EBxl/4whcmbZ8nUwAAMlGZAgCQicoUAIBMveoBqe95JqSb1WvZ\n95x0qp07dw7yWrEcVmzkDn86zJH5AzbHcqux5h9Nwm2E+87VJh8X8tcLz2HqyC6pzZJSlwsHafeb\nL/m50FhPTG3zoqVHdoldl6brt27duqHl/Dxpm89hiCdTAAAyUZkCAJCpV2HepSAW/pzmEOKshnWX\noljziqalhU1VAAAK7UlEQVQQWiysFwsH+6G3MHzrhxT9jtLDpjGpYepY6DW1dyR/3vnnnz80zy9/\nrKlFavMVPxxcYtDvcBtNxxk2efFD9rFz0yaUG7sOJbZx0kknDabDz5ffZCs8N/4gC6l4MgUAIBOV\nKQAAmahMAQDIZKPk6cxsepN6U6av+dOlkBstce7NbM45t6ZAcTqx5557un322UfSaLmuEs1GmoTd\nvflNLfo4UHSYq/NzuaeddtpgOsyf+vm52Kg5pYVNY84999zBdKzJS9Pno0S+swt+nvTiiy8eTPs5\nUimeu/ZzpnvvvXfSvcyTKQAAmahMAQDIRNMYAMm6DOW17TUnNfQcW65EE5rYaDZNYqHdWAi4RHjY\nD+3GBsDOHWA9pm2o2E8JhD0b+U1g/NDuKCPq+CMApeLJFACATFSmAABkIszbI319g3cWca77wQ/X\nhQNxp9q2bVvScqk9G8V+nxpePP744wfTfofyUpkQbYk3f/1epWJK9ErkX6O2b2n72/QHTg/flm4b\n2vWFA6Sn4MkUAIBMVKYAAGSiMgUAIBM50zEjVzc5nPvupDZxCHs5OuOMMwbTbQdo3rx582Dab+4h\nDffmk9rUJjUvGsvV+Tm3ZcuWDS0X/tykTV407OUoxm/Kk5q7bNsDln+dw2vUJGxq1HTuw89N1z1J\nNeHJFACATFSmAABkmvowbxehu9zO3AknlsF5HI8SHZanrhOGWv2QX2xAbV/YYbnfTGLLli1D8/yB\nrpsGG5eGQ4ptQrlSevi2S6k9KkmPHyx7Xqw3pNQepvxrIg2fR/8crlq1qnF727dvH/rZv7ZtUwJd\n4skUAIBMVKYAAGSiMgUAINPU50y7QK6umX9uSg8UznmfjEkO6uznLsNcaFOXbrt27Rr62c8Thuvk\n5tbCrgD9fbfNkY5zQHBfeCx+HvP8888feXthM6dYF3/+dUg95nC5WHOjXP5g4NLjP2MpeDIFACAT\nlSkAAJkI82KiCO1OxooVK7RmzRpJ8aYQk9QUyouF+Nr2NtSmTGFTk7bhyz7wm8n4vUZJzT1Hhb1D\npfZK1DbMHTv3TfvyhaFcfwDwsElVGzyZAgCQicoUAIBMhHmBEZR+g3lSdu7cOQjnlegBKVVssOZY\nqK3LNzmldqHH1N6GwuWaQpRdhH/9fYVv8/r83qHCQdr9MO9JJ500mA57UPLfxo4dS4njTD33fi9K\nqb1jtcWTKQAAmahMAQDIRGUKAECmqc+ZxnJYNLsob1ZyhjFL4Rh9o+RI/V5vwtyab9u2bQsuF9tX\n2AvPOAeDLp2vTM0ZjjKYd245wn35OcPYIOp+nvu0005r3Jc/AkzbZkMl+HnSDRs2DKZHGTg+7N0p\nBU+mAABkojIFACDT1Id5Y/xwHSFfYHR+UwhpOPQa60R+8+bNg2k/vOaHf6Xh8FrY846/Lz8kOW0D\ndodKd3Sf2htQONi2f12aejlqa5xh3fD4/SYw/jHGevoKm2yFn9MUPJkCAJCJyhQAgExUpgAAZOok\nZ1oiP1m6eQJNaBCz1JrDxPh50osvvnhonp+f9JtChPm4448/fjB98803D6bDLug2btw4mA7zdk1d\nvPnbC3/29ytJa9euHUz7Az7Huvhrm+9L3UaXzXDC/KE/UkrYnZ4/MHuM35zpyiuvHEyH17Jr/rHE\nRnxpavITEzbZomkMAAATQGUKAECmmW4akyoM8RH2fcwshT9n6VhKCkNafvjOD+tKzaO+xEKXYVg2\nVVNThjA8edRRRw2m/SY50nDznRKh11got8Q2csXC1374U2pu/hEb2af0QPKxnpLCwbxvuOGGwXRs\nxBc/tOv3vhVr/kPTGAAAeoDKFACATL0N8/qhVsJzaIvPzuLCDuv90Kj/xm4o1vOO3xOR/4ZtGIZd\nt27dYDq1h5rwzUs/JFf6DdMuOmxvs41YOboMG0tlBotvKmOsvP7b19JwuiDWs1EsTN2Et3kBAOgB\nKlMAADJRmQIAkKm3OVNf2FSFPBjQL7Hcl5+DDUd8ed3rXjeYDnNdfh4rlreLDVLua5Nb7KKnpDbG\n2aNS7Fz7ucRR8tOlyxjrvalEjrcNnkwBAMhEZQoAQKZiYd5p7jVomss+bTjX/bBixQqtWbNG0uND\nZn5TluXLlw/Nyw2Vhj0q+VLDc2E4eP369YPpMIzsl98fKLxtk5dxhnZLCwdzT22W5PND6rFzHSod\nHvd7vSoR1g2bwrQ5NzyZAgCQicoUAIBMVKYAAGSaiqYxXSB3h6Vs586dg9E1UpuWSPn5rrB7wtTc\nl58n9XOk0uNzdylix+GPVuLnWdtqm59NzTO2zUf6TVv8vHnsOvjXK8x/x85V6e4PY81yUgc9948l\n3J6fXyZnCgDAmFCZAgCQacmEeQnrthM7b9PQE1Xq6EOlPx/TcG5ShCN3pIY9U0N5qb3o+MuFYV0/\nJBc2z2gqbxh63b59+2B6y5Ytg+kwlOmPhpN6LlJDmeFg2CWa6/jzwvC43wTKFw627fOvwygjCqWO\nFOOvF17LpvKGTX785fzrFytvyP88pOLJFACATFSmAABkKhbm9cNasxpSbRu6m9XzMW3HNW3lHRd/\ncG1J2rRp02A6FlJt0/l6GD7ztxeG61L54cBYWNp/S3XDhg1DyzWFNsMwtH8+whDwKGHEeU3h5XD7\n4babenaKCXuO8kPW/rkPr1FT2HuU698U9h/lTee1a9cuWI7wvKWWKYYwLwAAE0BlCgBAJipTAAAy\nddI0JswtLvVcVV+aSSz164CF7dy5c+jnq6++unHZppxhLF8YywuWsHXr1sZy+HlSf99hjrSpl5v5\nXqIWWi/Mp/rzSvfCE/Kbf/i5xFF6bPKX9XPNsdxqiR6hfLHcamxZv4x++UJtyxvbZhOeTAEAyERl\nCgBAprH0gNSHMGcfyjBpS6H5EsrauHHj0M9+CLhNiHLz5s1DP6eGRv19hduIhUP9Zf19pXZeHobA\n/fXCEHC4bIpYp/KxkLLPP/7SYdiutpkitSlL2BxqUgO482QKAEAmKlMAADJRmQIAkGnJjBoDIE0s\nj+fnBf3lwvxh04Djfu5TSs9dttmXNNxVYmpO0292kXouSghzxrF9t+16sUnXecY2I+CkdjWY2nVh\nTGwUoVQ8mQIAkInKFACATIR5lxCawyBFLMzZNC8MeTaFQGOhy1Sp+2qrRBlL7Nc/12EI2O8ByQ9J\ndhGubRM27WJfqeVoE9oN1yHMCwDABFCZAgCQiTAvgCGxMGdqCDT1jdjU5Uoosa82IfC2/G2cdNJJ\nQ/OaBhkYZbDt0mHTEkqXI3aMpY+LJ1MAADJRmQIAkInKFACATORMAbQSyxGm5gzb5BbDwavH2ZSl\nRD45lX+c4agxfp60be6vy1xo6iDffcXg4AAATACVKQAAmQjzzjB6PEKKtmHTrkO0TWHktuHUSfVs\nFGo6rvDcrF+/fjB9+umnD81bvnx50r5Sm7+U7uVo2nti2rp168jr8GQKAEAmKlMAADJRmQIAkImc\nKbAE7bXXXoNBtf0BtCVpxYoVne13lH35y/rLrVu3bmi5cLBwn79sieVShWX0+dv386R+jnQxu3bt\nSlrOz63u3r07e7lUYc7Rb2rizxulCUqbMsbWiZXDb3qUiidTAAAyUZkCAJDJRmk+YWY/k3Rnd8UB\nZsYRzrmDJ12IJtzLQLKke3mkyhQAADweYV4AADJRmQIAkInKFACATFSmAABkojIFACATlSkAAJmo\nTAEAyERlCgBAJipTAAAy/X+jH+8/8mECQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f573c459110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, plots = plt.subplots(2, 2, figsize=(8, 6))\n",
    "plots[0][0].matshow(im_bw, cmap=cm.Greys_r)\n",
    "plots[0][0].set_title('Binary image')\n",
    "plots[0][1].matshow(im_noise, cmap=cm.Greys_r)\n",
    "plots[0][1].set_title('Noise image')\n",
    "plots[1][0].matshow(imgplot[-1][0], cmap=cm.Greys_r)\n",
    "# plots[2].matshow(result.reshape(im_bw.shape), cmap=cm.Greys_r)\n",
    "plots[1][0].set_title('Denoised image')\n",
    "\n",
    "marginals = inf.marginals(np.arange(nvar))\n",
    "exp_marginals = np.exp(marginals)\n",
    "marg_sum = exp_marginals.sum(axis=1)\n",
    "# print np.max(marg_sum)\n",
    "# print np.min(marg_sum)\n",
    "marginals = exp_marginals / marg_sum[:, np.newaxis]\n",
    "\n",
    "black_margs = marginals[:,0].reshape(im_bw.shape)\n",
    "white_margs = marginals[:,1].reshape(im_bw.shape)\n",
    "plots[1][1].matshow(black_margs, cmap=cm.Greys_r)\n",
    "plots[1][1].set_title('Black marginals')\n",
    "# plots[1][1].matshow(white_margs, cmap=cm.Greys)\n",
    "# plots[1][1].set_title('White marginals')\n",
    "\n",
    "for p in plots:\n",
    "    for s in p:\n",
    "        s.set_xticks(())\n",
    "        s.set_yticks(())\n",
    "plt.show()\n"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}