Image de-noising with Ising model

image_de-noising.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "from pandas import DataFrame, Series\n",
    "from numpy.random import randint, randn, rand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def spin_direction(field, image, x, y):\n",
    "    energy = float(-H - C * image[y][x])\n",
    "    for dx, dy in [(-1,0), (1,0), (0,-1), (0,1)]:\n",
    "        # Cyclic boundary condition\n",
    "        if x + dx < 0: dx += Size\n",
    "        if y + dy < 0: dy += Size\n",
    "        if x + dx >= Size: dx -= Size\n",
    "        if y + dy >= Size: dy -= Size\n",
    "        energy += -J * field[y+dy][x+dx]\n",
    "\n",
    "    p = 1/(1+np.exp(2*energy/Temp))\n",
    "    if rand() <= p:\n",
    "        spin = 1\n",
    "    else:\n",
    "        spin = -1\n",
    "    return spin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def show_image(subplot, image):\n",
    "    subplot.grid(None)\n",
    "    subplot.set_xticks([])\n",
    "    subplot.set_yticks([])\n",
    "    subplot.imshow(image, cmap=plt.cm.gray_r, interpolation='nearest')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Size = 50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "input_form = \"\"\"\n",
    "<table>\n",
    "<td style=\"border-style: none;\">\n",
    "<div style=\"border: solid 2px #666; width: 253px; height: 252px;\">\n",
    "<canvas width=\"250\" height=\"250\"></canvas>\n",
    "</div></td>\n",
    "<td style=\"border-style: none;\">\n",
    "<button onclick=\"clear_value()\">Clear</button>\n",
    "</td>\n",
    "</table>\n",
    "\"\"\"\n",
    "\n",
    "javascript = \"\"\"\n",
    "<script type=\"text/Javascript\">\n",
    "    var pixels = [];\n",
    "    for (var i = 0; i < 50*50; i++) pixels[i] = 0\n",
    "    var click = 0;\n",
    "\n",
    "    var canvas = document.querySelector(\"canvas\");\n",
    "    canvas.addEventListener(\"mousemove\", function(e){\n",
    "        if (e.buttons == 1) {\n",
    "            click = 1;\n",
    "            canvas.getContext(\"2d\").fillStyle = \"rgb(0,0,0)\";\n",
    "            canvas.getContext(\"2d\").fillRect(e.offsetX, e.offsetY, 16, 16);\n",
    "            x = Math.floor(e.offsetY * 0.2)\n",
    "            y = Math.floor(e.offsetX * 0.2) + 1\n",
    "            for (var dy = -1; dy < 3; dy++){\n",
    "                for (var dx = -1; dx < 3; dx++){\n",
    "                    if ((x + dx < 50) && (y + dy < 50) \n",
    "                     && (x + dx >= 0) && (y + dy >= 0)){\n",
    "                        pixels[(y+dy)+(x+dx)*50] = 1\n",
    "                    }\n",
    "                }\n",
    "            }\n",
    "        } else {\n",
    "            if (click == 1) set_value()\n",
    "            click = 0;\n",
    "        }\n",
    "    });\n",
    "    \n",
    "    function set_value(){\n",
    "        var result = \"\"\n",
    "        for (var i = 0; i < 50*50; i++) result += pixels[i] + \",\"\n",
    "        var kernel = IPython.notebook.kernel;\n",
    "        kernel.execute(\"image = [\" + result + \"]\");\n",
    "    }\n",
    "    \n",
    "    function clear_value(){\n",
    "        canvas.getContext(\"2d\").fillStyle = \"rgb(255,255,255)\";\n",
    "        canvas.getContext(\"2d\").fillRect(0, 0, 280, 280);\n",
    "        for (var i = 0; i < 50*50; i++) pixels[i] = 0\n",
    "    }\n",
    "</script>\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<table>\n",
       "<td style=\"border-style: none;\">\n",
       "<div style=\"border: solid 2px #666; width: 253px; height: 252px;\">\n",
       "<canvas width=\"250\" height=\"250\"></canvas>\n",
       "</div></td>\n",
       "<td style=\"border-style: none;\">\n",
       "<button onclick=\"clear_value()\">Clear</button>\n",
       "</td>\n",
       "</table>\n",
       "\n",
       "<script type=\"text/Javascript\">\n",
       "    var pixels = [];\n",
       "    for (var i = 0; i < 50*50; i++) pixels[i] = 0\n",
       "    var click = 0;\n",
       "\n",
       "    var canvas = document.querySelector(\"canvas\");\n",
       "    canvas.addEventListener(\"mousemove\", function(e){\n",
       "        if (e.buttons == 1) {\n",
       "            click = 1;\n",
       "            canvas.getContext(\"2d\").fillStyle = \"rgb(0,0,0)\";\n",
       "            canvas.getContext(\"2d\").fillRect(e.offsetX, e.offsetY, 16, 16);\n",
       "            x = Math.floor(e.offsetY * 0.2)\n",
       "            y = Math.floor(e.offsetX * 0.2) + 1\n",
       "            for (var dy = -1; dy < 3; dy++){\n",
       "                for (var dx = -1; dx < 3; dx++){\n",
       "                    if ((x + dx < 50) && (y + dy < 50) \n",
       "                     && (x + dx >= 0) && (y + dy >= 0)){\n",
       "                        pixels[(y+dy)+(x+dx)*50] = 1\n",
       "                    }\n",
       "                }\n",
       "            }\n",
       "        } else {\n",
       "            if (click == 1) set_value()\n",
       "            click = 0;\n",
       "        }\n",
       "    });\n",
       "    \n",
       "    function set_value(){\n",
       "        var result = \"\"\n",
       "        for (var i = 0; i < 50*50; i++) result += pixels[i] + \",\"\n",
       "        var kernel = IPython.notebook.kernel;\n",
       "        kernel.execute(\"image = [\" + result + \"]\");\n",
       "    }\n",
       "    \n",
       "    function clear_value(){\n",
       "        canvas.getContext(\"2d\").fillStyle = \"rgb(255,255,255)\";\n",
       "        canvas.getContext(\"2d\").fillRect(0, 0, 280, 280);\n",
       "        for (var i = 0; i < 50*50; i++) pixels[i] = 0\n",
       "    }\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML(input_form + javascript)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "original = np.array(image).reshape([Size, Size])*2 - 1\n",
    "noisy = np.copy(original)\n",
    "for _ in range(0, int(Size*Size*0.05)):\n",
    "    x, y = randint(0, Size, 2)\n",
    "    noisy[y][x] = -original[y][x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7YAAAElCAYAAADQl35cAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAADkxJREFUeJzt3DG23EQaBeDqOQTO2QCp09eLYA8DywP2wCLeS516A+Rk\nPcEbZgxI0IVU0n9L33eOA3NMu7ok/fY9Zd3b4/FoAAAAkOpfZy8AAAAAthBsAQAAiCbYAgAAEE2w\nBQAAIJpgCwAAQDTBFgAAgGiCLQAAANG+euYX3W63r1tr37bWPrfWfh25IGBaH1pr37TWfn48Hr+c\nvJbdmI/ADsxHgGVPz8engm17H0o/blwUQGutfdda++nsRezIfAT2Yj4CLPvb+fhssP3cWms//PBD\n+/jx48Y1AVf06dOn9v3337f233kykc+tmY97uN/vi//99fV102f0/P+82+Na8DzzEWBZz3x8Ntj+\n2lprHz9+bC8vL/98ZQDz/XM083GwrfvquuzHXg5nPgIs+9v5qDwKAACAaIItAAAA0Z79p8hc1O12\nO3sJ7fF4nL0EmMrSc732nPXMgK3P6h7P+pHzYub5aO7+tZ5nCK5o5vlIXU5sAQAAiCbYAgAAEE2w\nBQAAIJpgCwAAQDTlUfxPhRf9l6ytSykAlVW+b3vWUGG9FZiPfMn+wv+Zj1ThxBYAAIBogi0AAADR\nBFsAAACiCbYAAABEE2wBAACIphX5gqq218FMtC5mMh/PsbTvniGoxXykOie2AAAARBNsAQAAiCbY\nAgAAEE2wBQAAIJryKIAJKN/JtHSNrljQ4l4FYCsntgAAAEQTbAEAAIgm2AIAABBNsAUAACCaYAsA\nAEA0rcgAE+hpldWgvI+1PetpNT6yAXnt93LtAZiBE1sAAACiCbYAAABEE2wBAACIJtgCAAAQTXkU\nQFGjSp6UBe3jyOKnPcxy3Xv2XakawHU4sQUAACCaYAsAAEA0wRYAAIBogi0AAADRBFsAAACiaUUG\nKEojK0eYoQ1YWzgATmwBAACIJtgCAAAQTbAFAAAgmmALAABANOVRAExtqRypNWVBv0nbh6X1usYA\nOLEFAAAgmmALAABANMEWAACAaIItAAAA0QRbAAAAomlFXrHWsHgkbY5ARWnz0Sydy9L9t3aNe34t\n7CFtPsJMnNgCAAAQTbAFAAAgmmALAABANMEWAACAaFHlURVeyD/SHt935gKBqvfDzHvO+are91X0\n7M/Mz2rV+2SPPVccxpqq9/3R7MNfq7o/5tV2TmwBAACIJtgCAAAQTbAFAAAgmmALAABANMEWAACA\naCVbkau2lSWyl8fTygr/zNrzsPRM7THbjpyPPd9tZmnzcW29FdYGzCVtPlbkxBYAAIBogi0AAADR\nBFsAAACiCbYAAABEO7086mrFGfAlxST8lRnm49GFUFXN/N32sHSfVJiPPffvHr+W53mmuLIK87Ei\nJ7YAAABEE2wBAACIJtgCAAAQTbAFAAAgmmALAABAtNNbkYE/06LJmp7m1QqqrquytGu8h57vVmE+\n9vx+ZjdwlArz8UxObAEAAIgm2AIAABBNsAUAACCaYAsAAEC0w8qjKpdeVH2puvKeVdVzLe0vVWwt\nzhnJfDze2ndbuhZ77MOoz4U9uBeBZzmxBQAAIJpgCwAAQDTBFgAAgGiCLQAAANEEWwAAAKId1opc\nQdV2zzVr69UQ+G7r9UxrUO5pSoVeo+6jUfftFedjz3czH99VnY8V9gyu6si5UOFZT5uPWzixBQAA\nIJpgCwAAQDTBFgAAgGiCLQAAANGmLY+a8YXo36SVemxV4VpWLqpZWkOFPaOuI++Po+/FqvOxZ4ak\nPb8zz8ejys/e3t7a/X7f9JnsI+35G6XC89uj6nWbeT5W5MQWAACAaIItAAAA0QRbAAAAogm2AAAA\nRBNsAQAAiDZtKzLvltrNKjSxrUlrY6u6v6OaPGEmRz6/PZ876vk9ei7MMB/32JsZm0eZX9XndxZV\n9zf9749ObAEAAIgm2AIAABBNsAUAACCaYAsAAEC0+PKolJeZr8C1eLe2DxVKARjnfr//7ucVnoej\n1zBDSc6o57fCPlRewwzzMb105Wpcl34zz8cKZp6PR3FiCwAAQDTBFgAAgGiCLQAAANEEWwAAAKIJ\ntgAAAESLb0VOM6rZrKdRTvvcOZb2XdPdPF5fX9vLy8vZyzjV1tlSYT6O/AzW9exv1bnpHuGq3Ptj\n+fvj85zYAgAAEE2wBQAAIJpgCwAAQDTBFgAAgGjKowbyYjcAM1n6c21UcczMf4bO/N0AzuLEFgAA\ngGiCLQAAANEEWwAAAKIJtgAAAEQTbAEAAIimFRlOtNYmqjGTXlubadfuuZkbbys0/I76/UY5cr0z\nz8c/fre3t7d2v99PWg0JRt33aTOIdzPPxy2c2AIAABBNsAUAACCaYAsAAEA0wRYAAIBoyqMAJrC1\nCKmniKLnc3tKk44uvahQhAQA7MOJLQAAANEEWwAAAKIJtgAAAEQTbAEAAIgm2AIAABBNK/IOjm7y\nBHKNmhejWndHfe7Mc3NrkzRwDTPPQTiDE1sAAACiCbYAAABEE2wBAACIJtgCAAAQTXlUICUkkGvp\n+d2jQERh0bsK37nCGgDgapzYAgAAEE2wBQAAIJpgCwAAQDTBFgAAgGiCLQAAANG0IsOJ9mjDBZjR\nzPNx5u8GjGeGLHNiCwAAQDTBFgAAgGiCLQAAANEEWwAAAKIJtgAHut1uf/oBPO/xePzpxyhLz+se\nz+wf1//6+rrDavOZj1DT2iys9swKtgAAAEQTbAEAAIgm2AIAABBNsAUAACCaYAsAAEC0r85eAP2W\nGsdGtkICpDAf53dk66Z7ByCHE1sAAACiCbYAAABEE2wBAACIJtgCAAAQLb48aq1EQuEDUNHSbBpV\nhmM+XpcSrX49z6G9HMN8BLZwYgsAAEA0wRYAAIBogi0AAADRBFsAAACiCbYAAABEi29FhhSjmh2X\naHVkL6PadY9sP10zc3Nw1e9x9DXuUXXPgGurMDdT5qMTWwAAAKIJtgAAAEQTbAEAAIgm2AIAABBN\nedQO1l6orvCyN8AWWwsjRs3HPT43pQxjJhX+vHTdqaLC8wBfSp+PTmwBAACIJtgCAAAQTbAFAAAg\nmmALAABANMEWAACAaNO2Ii81yh3d9JXeLAbMaeb5aO4CW8w8H2F2TmwBAACIJtgCAAAQTbAFAAAg\nmmALAABAtMPKo9ZehF96SZ+xjt7zq5Ug2F96mY91VH5+R5XaHFmWU3l/qcl8rMPzO5b93c6JLQAA\nANEEWwAAAKIJtgAAAEQTbAEAAIgm2AIAABDtsFbkCtbaxmZsBQPoccX5mNaqunQt9rhuM1/jCo5s\nnWYM8xEyOLEFAAAgmmALAABANMEWAACAaIItAAAA0U4vj+opwxhlhmKHyi/521+uaFSpj/nYr/Lz\nu3V/K1yLPfa3wr0+SoVrNKMK94z5OJb9pZcTWwAAAKIJtgAAAEQTbAEAAIgm2AIAABBNsAUAACDa\n6a3IVe3RaDrKDA1r9rdfhb3heTNfL8/vWBX2N20fR+1Nz7WocN04X+X7IO25XlJhf9P2scK9dxQn\ntgAAAEQTbAEAAIgm2AIAABBNsAUAACCa8qhOSy+Mr72UnfZyeQU9e9bzMrxrwSwqz5sKa5jZzPOx\nwhqW9OzjlQpaqjIfr2vm+cjznNgCAAAQTbAFAAAgmmALAABANMEWAACAaIItAAAA0Uq2IldutVtS\ndV2zm3nftWuSaOm+nfk5rWzmfR/VaGruwjWYj/NyYgsAAEA0wRYAAIBogi0AAADRBFsAAACilSyP\nWqOYhNlc/SV/9lNhPprH7GmP+WjG0lqN+Qh7MtuWObEFAAAgmmALAABANMEWAACAaIItAAAA0QRb\nAAAAokW1Ii9ZawXTdrdujyY1+7tOU9013e/33/28wn1QYT5WWMMeeq5n2nc7UoXnAlqbZzZVYD7u\nw3zczoktAAAA0QRbAAAAogm2AAAARBNsAQAAiBZfHrVm5hfZK7xcPvP+Lqmw59T2+vraXl5ezl7G\nU458fo9+/is8q1v3N63UpsKer+nZX2ht7r/fVLj3Z97fJRX2/Eqc2AIAABBNsAUAACCaYAsAAEA0\nwRYAAIBogi0AAADRpm1F7qGxbCz7C7k8v2P17O/Sr11rDT3yulVYw5oKa9DMPC/XcSz7Sy8ntgAA\nAEQTbAEAAIgm2AIAABBNsAUAACCa8igAulQuC7qaCnves4aj750KxU0VrhHAFTixBQAAIJpgCwAA\nQDTBFgAAgGiCLQAAANEEWwAAAKJpRQagi5bX/VRo7T2SRmIARnFiCwAAQDTBFgAAgGiCLQAAANEE\nWwAAAKIpjwJgGlXLmJbW1dq4tVXdh8rsGUA2J7YAAABEE2wBAACIJtgCAAAQTbAFAAAgmmALAABA\nNK3IAEXN3NI66rtt/YxR7cVHXzdty/2WvsfRbdYA/HNObAEAAIgm2AIAABBNsAUAACCaYAsAAEA0\n5VEARc1cUFP1u1VdVxVXK6XaYw1VvxvAbJzYAgAAEE2wBQAAIJpgCwAAQDTBFgAAgGiCLQAAANG0\nIgMwtaVW2tY001Yy87WY+bsBVOLEFgAAgGiCLQAAANEEWwAAAKIJtgAAAERTHgXA4ZYKnUaV7Cjv\nAYD5ObEFAAAgmmALAABANMEWAACAaIItAAAA0QRbAAAAomlFBhhgqfW3NQ29v6mwDz3NzEe2OAMA\n/ZzYAgAAEE2wBQAAIJpgCwAAQLRn37H90Fprnz59GrgUYGZfzI8PZ65jgK75+Pb2NnQxbNNzfVxL\n9mI+AizrmY+3Z8ovbrfbv1trP25bFkBrrbXvHo/HT2cvYi/mI7Aj8xFg2d/Ox2eD7dettW9ba59b\na7/usjTgaj601r5prf38eDx+OXktuzEfgR2YjwDLnp6PTwVbAAAAqEp5FAAAANEEWwAAAKIJtgAA\nAEQTbAEAAIgm2AIAABBNsAUAACCaYAsAAEC0/wAw78Fss4jP+QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x556b490>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "J = 8.0\n",
    "C = 14.0\n",
    "H = 0.0\n",
    "Temp = 0.01\n",
    "\n",
    "field = randint(0, 2, [Size, Size])*2 - 1\n",
    "iternum = 50\n",
    "for _ in range(iternum):\n",
    "    lattice = DataFrame([(y,x) for x in xrange(Size) for y in xrange(Size)])\n",
    "    lattice.reindex(np.random.permutation(lattice.index))\n",
    "    for x, y in lattice.values:\n",
    "        field[y][x] = spin_direction(field, noisy, x, y)\n",
    "\n",
    "fig = plt.figure(figsize=(12,4))\n",
    "subplot = fig.add_subplot(1,3,1)\n",
    "show_image(subplot, original)\n",
    "subplot = fig.add_subplot(1,3,2)\n",
    "show_image(subplot, noisy)\n",
    "subplot = fig.add_subplot(1,3,3)\n",
    "show_image(subplot, field)"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}