{"nbformat_minor": 0, "worksheets": [{"cells": [{"source": ["Color Image Denoising with Median Filtering\n", "===========================================\n", "\n*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_matlab/) for details about how to install the toolboxes.\n", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$\n", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$\n", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$\n", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$\n", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$\n", "$\\newcommand{\\norm}[1]{\\|#1\\|}$\n", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$\n", "$\\newcommand{\\choice}[1]{ \\left\\{ \\begin{array}{l} #1 \\end{array} \\right. }$\n", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$\n", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$\n", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$\n", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$\n", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$\n", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$\n", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$\n", "$\\newcommand{\\CC}{\\mathbb{C}}$\n", "$\\newcommand{\\RR}{\\mathbb{R}}$\n", "$\\newcommand{\\EE}{\\mathbb{E}}$\n", "$\\newcommand{\\Zz}{\\mathcal{Z}}$\n", "$\\newcommand{\\Ww}{\\mathcal{W}}$\n", "$\\newcommand{\\Vv}{\\mathcal{V}}$\n", "$\\newcommand{\\Nn}{\\mathcal{N}}$\n", "$\\newcommand{\\NN}{\\mathcal{N}}$\n", "$\\newcommand{\\Hh}{\\mathcal{H}}$\n", "$\\newcommand{\\Bb}{\\mathcal{B}}$\n", "$\\newcommand{\\Ee}{\\mathcal{E}}$\n", "$\\newcommand{\\Cc}{\\mathcal{C}}$\n", "$\\newcommand{\\Gg}{\\mathcal{G}}$\n", "$\\newcommand{\\Ss}{\\mathcal{S}}$\n", "$\\newcommand{\\Pp}{\\mathcal{P}}$\n", "$\\newcommand{\\Ff}{\\mathcal{F}}$\n", "$\\newcommand{\\Xx}{\\mathcal{X}}$\n", "$\\newcommand{\\Mm}{\\mathcal{M}}$\n", "$\\newcommand{\\Ii}{\\mathcal{I}}$\n", "$\\newcommand{\\Dd}{\\mathcal{D}}$\n", "$\\newcommand{\\Ll}{\\mathcal{L}}$\n", "$\\newcommand{\\Tt}{\\mathcal{T}}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\al}{\\alpha}$\n", "$\\newcommand{\\la}{\\lambda}$\n", "$\\newcommand{\\ga}{\\gamma}$\n", "$\\newcommand{\\Ga}{\\Gamma}$\n", "$\\newcommand{\\La}{\\Lambda}$\n", "$\\newcommand{\\si}{\\sigma}$\n", "$\\newcommand{\\Si}{\\Sigma}$\n", "$\\newcommand{\\be}{\\beta}$\n", "$\\newcommand{\\de}{\\delta}$\n", "$\\newcommand{\\De}{\\Delta}$\n", "$\\newcommand{\\phi}{\\varphi}$\n", "$\\newcommand{\\th}{\\theta}$\n", "$\\newcommand{\\om}{\\omega}$\n", "$\\newcommand{\\Om}{\\Omega}$\n"], "metadata": {}, "cell_type": "markdown"}, {"source": ["This numerical tour explores denoising of color images using a local\n", "multi-dimensional median. This is the sequel to the numerical tour\n", "<../tv_median/ _Outliers and Median Denoiser_>."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 2, "cell_type": "code", "language": "python", "metadata": {}, "input": ["addpath('toolbox_signal')\n", "addpath('toolbox_general')\n", "addpath('solutions/multidim_7_median')"]}, {"source": ["Multidimensional Median\n", "-----------------------\n", "The median of |n| real values |x| is obtained by taking |v(n/2)| with\n", "|v=sort(x)| (with a special care for an even number |n|).\n", "It can alternatively obtained by minizing over |y| the sum of _absolute\n", "values_.\n", "\n", "\n", " |\\sum_i abs(x(i)-y)|\n", "\n", "\n", "This should be contrasted with the mean that minimizes the sum of\n", "_squares_.\n", "\n", "\n", " |\\sum_i (x(i)-y)^2|\n", "\n", "\n", "This allows one to define a mutidimensional median for set of points |x(i)| in\n", "dimension |d| by replacing |abs| by the d-dimensional norm.\n", "\n", "\n", "We define a Gaussian point cloud in 2D."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 3, "cell_type": "code", "language": "python", "metadata": {}, "input": ["d = 2; % dimension\n", "n = 1000; % number of points\n", "X = randn(d,n);"]}, {"source": ["We modify some points as positive outliers (to shift the mean)."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 4, "cell_type": "code", "language": "python", "metadata": {}, "input": ["p = 100; % number of outliers\n", "sel = randperm(n); sel = sel(1:p); % index of outliers\n", "X(:,sel) = rand(d,p)*50;"]}, {"source": ["We can compute the mean point."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 5, "cell_type": "code", "language": "python", "metadata": {}, "input": ["m = mean(X,2);"]}, {"source": ["To compute the median in 2D, one needs to minimize the sum of norms.\n", "This is not as straightforward as the sum of squares, since there\n", "is no close form solution. One needs to use an iterative algorithm, for\n", "instance the re-weighted least squares, that computes weighted means.\n", "\n", "number of iterations of the method"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 6, "cell_type": "code", "language": "python", "metadata": {}, "input": ["niter = 30;"]}, {"source": ["initialize the median using the mean"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 7, "cell_type": "code", "language": "python", "metadata": {}, "input": ["med = m;\n", "energy = [];\n", "for i=1:niter\n", " % comute the distance from med to the points\n", " dist = sqrt( sum( (X-repmat(med,[1 n])).^2 ) );\n", " % compute the weight, take care of not dividing by 0\n", " weight = 1./max( dist, 1e-10 ); weight = weight/sum(weight);\n", " % compute the weighted mean\n", " med = sum( repmat(weight,[d 1]).*X, 2 );\n", " energy(end+1) = sum( dist );\n", "end"]}, {"source": ["We can display the decay of the L1 energy through the iterations."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAARm0lEQVR4nO3d0Zab\nuBIFUHRX/v+XdR+8whAESGALU/LeD7M6dBmbMOnTJYRIOecJAKL537c/AABcIcAACEmAARCSAAMg\nJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAA\nCEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBg+wlKaU\nvv0hAOjgz7c/wB2SEAM4Kef87Y9QMXiA5TylFOA0lFJKET/2WQ5zJA5zJCF+7x98CBGAUQkwAEIS\nYACENP5gbkrT6IcI8GEhLvXpwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYAB\nEJIAAyCkXwiwEM+1AeCcXwgwAAYkwAAISYABEJIAAyAkAQZASAIMgJD+9H6D5Rz2zQdU31AAwHj6\ndmCrO7DKG7KOC1JKbxYAMKqOHdgrS+aWaC+9Dgpe3i8AYDy9OrBVOM1fzxnTWLC355YCAAZ26ySO\nnPPxNaqyoEy48iXHBQAMKcAkDgAo9Q2wcobFKqKqBR//JDISYFO46y/dO7DVDIsyoqoFn/0YAGxa\n/pwMEWZ9r4G5ggVAJ1biACCkSAFWbWlD9LwAfESvAFvd1DUV6XKhYLm9pQCAgd06jX7aSpezBdW3\nAOAXdBxCLKPo4CblTgUAjKpvB1aNkxsKABhSpEkcl+U8GWUEGMxPBBgA4xFgAIQkwAAISYABEJIA\nAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQk\nwAAISYABEJIAAyCkXwmwnKeUvv0hAPicP73fIC1yI+d8XLZZUN1D41sAMJK+HVj6t+tJJ5uglNLx\nHqoFAIyqYwe2aqoOouU4dap7aHkLAAbTqwMrhwRfX5cZk1KqDgxubq8WADCwWydx5JzPXqMqI/Bs\nAQBD+vIkjr32CwCO9Q2wcobFMq5uG+v7+0ZZXgLsCXf9pXsHtpphsYqQe+Lk7+U3Y4wAu77SYLyj\n7zWwg6zSDAHwju4d2IHNGYmTuRgANIi0lFS1pQ3R8wLwEb0CrLzra5Uu+V/LjXt7WG5vKQBgYLdO\no58upYvGC4BSxyHEMqveuYt5c8v7bwFAUH07sPY42aus7kFiAfymSJM4AGAmwAAISYABEJIAAyAk\nAQZASAIMgJB+KMByntzxDDCMHwowAEYiwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQk\nwAAISYABEJIAAyAkAQZASAIMgJAEGAAhCTAAQvrT+w3S4iGSOeeD724WVPfQUgDAePp2YKt8Ov7j\nZkF1D9V9AjCkjh3YK0vmlmgvaVYFKaVVF7W3h/YCAMbTqwNbpdf89cGY4Sq39tJo3l4tAGBgt07i\nyDmfvUZ1kHCNBQAM6WuTOL4SNjlPKU1iDmAAfQOsvO51kFvlqGOfjyG+ADaEu/7SvQOrztFYfrdT\nW7bcbbQTBHCTf39UBvhZ2TfAVmOGm38jXaMLgFF178AOlDPpAaDRNwNsOhld1ZY2RM8LwEf0mkZf\n3vV1dsmMzfvGpkXmVQsAGNit0+intruVG+9obi8AYDwdb2QuO6F37mLe3PL+WwAQ1H2zENu/dapM\nYgH8Js8DAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAkAQZASAIMgJAEGAAh/VyA\n5Tx5+grAAH4uwAAYgwADICQBBkBIAgyAkAQYACH9Ofhe+jtdL+d8y4cBgFZHHdicWymlZO45AE9y\n1IFNfzPslV4aMgCeoxJgL8tWbP6vGAPgi5oCbJOGDIAvag2w1TWw1dDiQYYtX7hZdkMBAOOpBNhm\nbi3/eDy5Y/XdVdSVr/14AQCj6jiNfnWpbC/qbigAYDyVafQv1ZpyeznRYznqOO2HzQcLABjYrStx\nlGlXJlz5kjcLABhS0xDigWpmmGEBQA/Xp9G/tE9BrBb3I0QBqsJdfzkKsHmSYTnxr30afXkT9P0R\nIrQAqo7neD/Q0TWwzRU3ltFVDQYXqADoxONUAAjpQQFW7VjfL3jJeYrQHANwpFeAre76mop0KQuW\n2z9SAMDA6rMQD1ZvarkSdrwYVVlQ3cOFAgDGc+KBlsvHWrZ0OWXNwV3MnQoAGNWJB1out5RfH7z8\nuwUADOncAy0B4CEq94G5vATAMz1oGj0AtGuaxAEAT1Nfjb76GEkAuJ8hRABCqqxGf9vnAIBTdGAA\nhNQUYMtlOEysB+AJKgG2eSuYDAPg6+oPtJy2lo+SYQB8V30I0VQOAB7IJA4AQhJgAIRUD7DV5S5X\nvwB4gtYHWq6+cGEMgO+qdGA5Z089BuCBfvSBljlPKU3DHRbADzGJA4CQKh3YwZSN8doyAAJpWokD\nAJ7mR6+BARCdpaQACKmpA7tmcwRyFYerms2wXNZcKwBgPN+chVh9UEv5MJezBQCM6qgDyzm/EuJy\nW3PwwnJRj/mBmatXlQuC7L2L9AL4HUcBVq4gtfSRwbrVk8aWb7SXRnPCVQsAGFj3IcS0cOHl5bM0\nzxYAMKTKEOLl/c7jgauN8z7L8clXsRACoEXHWYgvB5e45gwriz/LNEWAqnDTCJoCbDnh4tQVpoNL\nXNPOLMQeASO0AKqWPypDhNmVtRBbYqZaUA4Y7s1CBIBS01qI5USJHuEstwBoF2kpqWpqhuh5AfiI\nXtPo99bIaLm1+WWv1VvOATkuAGBgvWYhljMMGwvK+NF4AVCqd2CXFxssV4TaWyPq4CVvFgAwqsqU\nv+ryg8+3N60xpSnOQQDcKsSE8EoHVrZNU6j0AmBUv/tE5pw1YQCBffN5YABwmQADICQBBkBIAgyA\nkAQYACEJMABCEmAAhCTAAAhJgAEQkgADIKT1UlKXF5sHgDutOzCxBEAIG4v5vh41OQkzAB7MNTAA\nQtoOML0XAA+324HJMACezBAiACFdCbCUUvtsewDoQQcGQEgCDICQBBgAIW3cyPwpm9fJysmNy7LN\nqY/vF+zJeUppMt0SIKIvd2CrkCv/+GYBAKPq2IFNtX5otWDVXva8XwDAeP4JsDsDoFxu8bUGY0rp\ntXHvw7QXADCwvh3YdPICVVlQJtzZAgCG9E+AfbBxeQVJeYHqwiQOACh178BWF6hWGdaScO+TkQBV\n4Uaw7pvEsTe+d5xwH/8YAGw66C6eqdc0+pxzS2ysEq7ThwFgPFbiACCkSAFWbWlD9LwAfESvANtb\nI2MeJyzv5VrV793sdbCHVQEAA+s1iWO+K7lauZlbBwXVPQDwCzoOIa6iqJzWcXzb8kcKABjVN9dC\nvKcAgCFFmsQBADMBBkBIAgyAkAQYACEJMABCEmAAhCTAAAjp1wMs58k6HgAR/XqAARCUAAMgJAEG\nQEgCDICQBBgAIQkwAEISYACEJMAACEmAARCSAAMgJAEGQEgCDICQBBgAIQkwAEISYACE9KffrtPW\ng7ZyzgfFm99d7udaAQDjeXQHllJapWD5x+MCAEbVsQObmvuh49SZd7JXVi0AYDzf78BSStWBwc3t\n1QIABtY9wNLChZcvs20z56oFAAypV4C94qp6gWqv/QKAY32vgU3FBaplYt021nc8TTHnKaVJjAI/\nLtz1l/smceScy7+de9ovTR5A1fJHZYgw6xVg1cwweAjAO7oPIR7YvCQ2aZgAaPD9afTtqi1tiJ4X\ngI/oOAtxcwri3F3lfy03zl9PW5nUXgDAwDpeA7t879eKxguAUschxFUntOyuru1hc5/VlwAwpEes\nhXhQWd2DxAL4TZEmcQDATIABEJIAAyAkAQZASAIMgJAEGAAhCTAAQhJgAIQkwAAISYABEJIAAyAk\nAQZASAIMgJAE2DRNU86TZ4oBxCLAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAjp\nT79dp63FLXLOBzWr75Y11woAGM83O7Ay4VZbUkrlllMFAIyqYwc2HfZDc9LMNa8tKaXVq1YFB+8i\nvQB+x5evgS2z6nh0sdxeLQBgYN0DLC1cePlBwjUWADCkXgE2jweWG19yzvIGgMv6XgObGi5xzV4F\nPVLNNEWAqnDXX+6bxJFzPvjb6Zde/XYLMJLlj8oQYdYrwNozo2t0ATCq7kOIB8qZ9ADQ6JsBNp2M\nrmpLG6LnBeAjOs5C3JyC2H7T8auyLJv3UC0AYGAdr4G13Pt1Yb3Exp2clfOU0iT7AKLoeCNzuSLU\n2d6orN9bZepgCwBD+tpaiI1JUy2TWAC/yfPAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAA\nhCTAAAhJgAEQkgADICQBBkBIAgyAkAQYACEJMABCEmAAhCTAAAhJgAEQkgADICQB9p+cp5S+/SEA\naCPAAAhJgAEQkgADIKQ//Xadti4o5ZwPysrvfqQAgPF8swNLKa1CrvzjmwUAjKpjBza19UNzzV72\nvF8AwHi+1oHthc28/f0CAAbWPcDSQvndZYu22a69XwDAkHoF2CuuXKACoJO+18Cm4gJVSun+Psk0\nRYCqcD3GfZM4cs7f+tsRWgBVyx+VIcKsV4BFzIzX+UppCvjZAX6OlTj+I7cAAnlQgFU71vcLqnK2\nJj1ADB1nIW5OQZyHFl9flJHzwYLLXhkmxgCerOM1sL17v1ae0HiVXiHoehjAY3UcQlx1Qjnnckv1\nJW8WvMlwIsBjfXktxBsK3vTKMH0YwNM8aBLHY+nDAB5IgDUxrQPgaQRYKzPsAR5FgJ0jwwAeQoCd\nJsMAnkCAXTFfEuuXZCFW0nyfwxyJw+RmAuyieWK9yR0AX9H9eWDDm9fs2NwOQCcC7LplRJVxtdmW\nnUm1jgHY9dbskzuPepgnOcx3nTrMzn8nTznMzp7yOQ584fnIN3vMgPXgf8/AeB6eD+N3YEET+hW7\njZ/9QnHXnZ/6TbnrJ+m68+cc5hM+SdedP+2TdN35ow7z4cbvwAAYklmIAIQkwAAISYABEJIAAyAk\nAQZASONPo49l8661wWaKprQ793V5+NGPevMwhzm/qwPZPIToZ7N6jGOczZZTOT31bOrAuNXefeUp\npdW3HnMH+hWhP3xVeXTluYt+NqvHOIaWw3zy2dSBPc5zfrv5rJb/6edjf86/kLOqnzz6+Z0PcHWy\nyo4z7tm8cIwRvQ5qeQh7hzk99WzqwLjD8f/0B21Zn4/TS7gPfNnyB9zqh90wZ/PgGEfVOEx6sP1m\nOrAneuZw8ztafn1b/bx4yL+QUxp/Sx3v/JYGOJuN4p7N9k/72LOpA3uQuX8vNzKAMc5vzjnWj+kL\nWo5xjLO5Ug4qPpwO7ImqI++ENtj5DfdT74KDYxzgbMZtInVgj/ODI+8/ZbDzK702v44l/zVtzTl8\nMh3Yg8T9B0CLwc6v6Lr943T3qOtbLQQYcE45y3w8v3CMAxBgPFGsXwN/0Kkf60HP5vDRde1y3aPO\npmtgD7J3x/vw/5DmwffN7cMY4/w23qkd+my2/Iwe8mzunbXHnk0d2IO8BqAf9QvOzcY+9pHOb8sy\ngNGP9PgYBzibe4dQXdnyOXRgz1IuVPOQ33R6O1i6ZiQ/cn6dzShaztSTz2a8WxYAYNKBARCUAAMg\nJAEGQEgCDICQBBhseMgM6Yd8DHgmAQZASG5khib3rLOwehd3ucABHRgAIQkweC7XwOCAIUSoWEbI\ncgHvVbQsh/vmkcBy4LFlvdSDZcLb33SzBkaiA4Mryhy6tmVvY+83hQHowKBi2SGV67e+vmhZnX1z\n4+rBiQdTRcpHLL62lO3a8aeCYejA4LQyZvYem7Ty/mhe+aZnC2AYOjC46HJ/ozGCjxBgcJ+WOR1A\nIwEGF50doIv4yHl4MtfA4KJV83Thnq0L7ddqtv3Zl8NIdGBwwmvK33yv1YU7rvZm0q/GFTef4375\nTWFIOjBosjdV/WBLteDsTi68KQxs94Z/AHgyHRgAIQkwAEISYACEJMAACEmAARDS/wFR1OLEh5jZ\nsAAAAABJRU5ErkJggg==\n", "output_type": "display_data"}], "prompt_number": 8, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "plot(energy, '.-'); axis('tight')\n", "set_label('Iteration', 'L1 energy');"]}, {"source": ["We can display the points, the mean and the median."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [{"metadata": {}, "png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAGwCAIAAADOgk3lAAAACXBIWXMAAAsSAAALEgHS3X78AAAA\nIXRFWHRTb2Z0d2FyZQBBcnRpZmV4IEdob3N0c2NyaXB0IDguNTRTRzzSAAAQL0lEQVR4nO3d63Li\nyLaFUeWJ/f6vnOeH2ipZAixAt5k5RnR0uNx2Gyisj5W6UGqtAwCk+b+rbwAAfELAAIgkYABEEjAA\nIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRg\nAEQSMAAiCRgAkQQMgEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAsbFShlKufpGAIEa\nD1ixaQRoVOMB4/5qvfoWwDZWC+5GwLiehhFEw+7jf1ffAIAMXmndjQkMgEgCBkAkAQMgkoCRxGFg\nwETAAIjkKESSOAwMmJjAAIgkYABEsoTYu+mYCKtzQBYTGNAah6p2wgTWO4MXjRnrVYrndvtMYECD\n1KsHJjCgKdLVDxMYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwYBi8WSiBBCyYLQ7QM1fiAIbB\nBSwIFB+wUkqtdfxg+mTt43exj3sJ8Fh2wMrvFbROugUtcdl4Pha8D2yaveafKXYKQY7prU/gA9kT\n2MK0ljgPmykNbqtW9eJzqQEbszT9u9b6rEyKBXfmF5SPpS4h1h/jx1YOAXqTOoEtzBtm5KIl4/Pa\nkxrW4gM25Uq3OIiEwD3FBwwOdfnitHDCMwIGr+gH3FbqQRwAdE7AANrU/PW+BQwa1PyWiz/18ASw\nDwxa08OWiz/1sPtWwKA1PWy5YLCEyEMWoCBUV7+5AsZSV78A0JLeru5vCZElC1AQqrer+wsY3IjL\nVvGlrp48lhDhDHYrwu4EDM6jYbAjAduH19e8Ni7s/Lm8U2tfS0DwDQGDkygT7MtBHPuwbQI4mQkM\ngEgCBkAkAQMgkoABh+vhMN3m7+ANCRjAPjTsZI5CBA7nMF2OIGC0zKUFOY2n2fksIfKGHvZkAClM\nYLTMi2JomIDxBj0A7iN+CbH8LGmVH9feHgDOkT2BzetVf6aD+ccAtCo4YGOo/hy5Fl+gbQBtCA7Y\nRooFfMyZGHeWug9snKvm/4ZOOJkBRqkT2GKPl4YBRzB73VlqwBbmDbNmCNCD+IBNudItOuGZDqPU\nfWAAdE7AgDc4hIT7EDAAIgkY8B474bgJAQPeoF7ch4ABEEnAAIgkYBzIEWvAcQSMw2kYcIT4K3HE\n6eri1p3cTeASJrBrGEoAvmQCO5uhBGAXJjAAIgkYQLCej/UVMHis5+0CRBAwACIJGDzmcBsi1Nrv\nc1XA4Kl3tws7rjpawIQ/CRjcjnTBFgIGu9lrMafbFSF4ixOZ4Y40DP5kAgMgkoABEEnA3uPYMICb\nELA37JIuCQTYRfBBHOWnA7XW+R+nz+zOfnWA+wgO2DBL1/jBQd3aXcjNBLi14ICtczUOYTfP2L1v\nHXCZaRXJVmKj4IANv5cNh9VA9uJrAEiXHbD53q9nZVIsIIJt1btSj0JczFXFgX0AnUmdwGqtz45C\nPGHkGn+Ul0sAF0oN2LAKlaVCgK4EB+xCWglwudR9YAB0TsAAiCRgAEQSMAAiCRgQzNs79EzAAIjk\nMHogmHNaemYCAyCSgAEQScAAiCRgAEQSMAAiCRjQMieKNUzA4Aw2o7A754EBLXOiWMMEDM5gMwq7\ns4QIQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQMgEjBJzKXnyvz1FrXfwSgbcEBG36na+pWKUXD\nAJoXvISoUgA9y57Ayobrey++5szsjT/5oB9YisvrAV3LDthi79eLr7nQEaUZ77GGAdtNW8pmthup\nAYvY0XX7Gwj0LrpqqQGrtfZ8FGIHdxHYWXvbjdSADatQ9dAtgH1FbziDj0IEoGcCBkAkAQMgkoDR\noFKGDacItumS+97zA86FBIzW2JJCJ4KPQoSHoo+q+t4ld7/zx5yrmMAAiCRgAESyhMiBoq9SA9yc\nCYzDqRdwBBMYB5Iu4DgmMH5xQg+QQsD4R7qAIJYQ+ceKHxDEBAZAJAEDIJKAARBJwACIJGAARBIw\nACI5jB6gX9EXLDWBQTtcSIWumMCATUqJfJHOa9F/pyYwaEetx26PjHfcioABEMkSIrBJ9FoTTQoO\nWPlZzqi1zv84fQZoQ/SRchwnOGDDLF3jB7oFTarV7rdTpRywE7wPbJ2rUkrxNIcWHX18CpNxIxqx\nKc2ewIbZ+DWsBrLpC+Zfb0oD3tLnAmbEnQ0O2FimRb3WFAtgu6BNZnDAhlmcFlMX8No4Vfil2cKj\ndFupARvHr/mBiIuDEgFoW2rA1pXSLYCupAYM+IbXe/dnmfdPwYfRwzlc4h3uyQQGcEdmrz8J2FfM\n+D3w9wv3ZAkRgEgmsK94bQ5wFRMYAJEEDIBIAgZAJAEDIJKAHeXCs1+deAv0QMAOoR8ARxOwQ4yH\n1191kL33rgV6IGBHkRCAQwkYAJEEDIBIAgZAJNdChKXpIFI7MuHOTGAARDKBwZLBCyKYwACIJGAA\nRBIwACIJGACRBAyASMFHIZafs3Vqres/AtC24IANv9M1dauUomFtK8WR7kDyEqJK9WmctL3jGpA9\ngQ0/81Z5vj1b/CfZa4O/RiA4YIuVw2cUqzH+PoFR8BLiIE4AHUudwMbxa37koaMQAbqSGrB1pXQL\noCvZS4gAdEvA4CSlOPof9iRgAERK3Qd2Ie83z2c8YWBfJjAAIpnA3uZ1NMAdmMAAiCRgAEQSMAAi\nCViwZycVOd8I6IGAZVuHSrqATjgKsTUOkoTPjC/+/AYFEbBgftO4CZt+LiFgAMMgwIHsAwN2YOvP\n+QQM+JZ6veCo4OMIGACR7AM7hH3awMh24DgmsP1ZLkhknQfimMD25wUXwAkEDIbByw4IZAkRgEgC\nBkAkAQMgkn1gcJTpsEY72OAIJjAAIsVPYKWUWuv4wfTJ6hUvl3ImO5wgOGBldd6pbsE5SpFnrhcc\nsPXgNX68yNiicyL3GSMFk/HJoGFcLjhga1PS5pVSLM7Xw5Ouh/vIzbUTMKGCc/hV4yYaOQpxvT+M\nfdVqswXcSyMTWK11aphRDKAH8QGbcqVbAF1pZAkRgN4IGPflTSbnPBSwIGAQQL1gTcAAiBR/EAcN\nc1zOxEMBayYwACIJGACRBAyASAIGQCQBAyCSgHFTzmL+nseQtgkYQMsafh3jPDBolrPHaJsJjJsa\nN76tvnKE0zT8Zn4CBkAkS4jcV6svG4FdmMDoTsP7tKErAgZAJEuIdMfKJLTBBAZAJAEDIJKAARBJ\nwHiD4/eA+xAwgLvwGvEtjkLkDY7fA+4jPmCllFrr+MH4mWory0XG56AnIB/z5HlLcMDKbNKeMrb4\nGIBWBQdsMXjB5bxwgjMFB2yjReEMZ8C1LDXvpf2AKRZAkxxGD3Cq7e8waQ/Ja41MYLVWRyECLRk3\naaVYbHwqPmBTrnQLONr5u69s2F6IDxhAk6TrTwIGsJWo3IqDOACIZAKjKdNRW14pQ/NMYABEMoHR\nFIMX9MMEBkAkExgAv6TsSzaBARDJBAbALzcfvCYmMKBxpbgqbpsEDIBIlhCBxqUsiPEuExh9sZoE\nzRAwACJZQqQvVpOgGSYwACIJGACRBAyASAIGQCQBAyCSgAFdu+rUQOcjfk/AAM421kvDvuQ8MKBr\nF54a6KzELwkYLRhfydockMJzdRdNBazMBvLqCQLQtKYCNuhWr/y1Q4daO4ijlFLsGAXoQJsTWCll\nGsUWPTOiAbShqYA9jJNi3VMp1v2Ar7SzhGjlMIiTYIDvtTOB1Vqnhpm6bq5WExjwrXYCNuhWFH9X\nwJfaWUIEoCsCBkAkAQMgkoABEEnAAIgkYLdw1VvqAeQSMPjFKwlIIWC34KSom3CJEAgiYP8ZF/HW\nW675Jw/drmnYHYx/C/4uIIKAvTIv1uK1+bt7rUr542qNXvXfhHpBCgH7T63//bP45LB6Vf5Baebp\nepixtuv1bLq9v9xbDj0QsF/WW6t50tavzTdv3eowlOH3G5UB8I3Gt6fvBmOq0ZZvmqfrxdc/HLna\nftgBTtDU1ejPtH1ZaXqT6PkfAfiSgP0yvk/V8GQUezafrXePrSM1vl2ZegHsRcD+q85Ulqlh8//6\n8FsAuJCDOJYWPdv4LVu+y/gFsCMT2ONFwl3+bwAcp/0JbPt5PH+uHG7/9rd+LgAfMIH9s9j79YGH\nV6IajGXs6q2TPaBh7U9gDz277OEL32wsnv2fTWl8TL2g/Qns5SnG/77gz5C8c+LXqzra7vAlTyEY\n9TiBbdxZ9c1m4s9qTu20MQL4THcBW1/q8FlCPlvcG8ev19/rTacaYPkXLtd8wH7Vab2Ot/tmaHGB\nRANWk6QL7qDxixu92Mv1wTGHdZi9K8rw9+P27Ayzph9ygJM0NYGVH0f8z8d6laGO6ZrH7PntWX7w\n8I8M3nkLeF87ARsvlTuaNazML3I4X9P75FTln6lry/i1/ilWFAF21P5h9Acp/61B/l0k5zJv4fEB\n3tV+wBYrinvt89tYr9nN+G8bbUsNsIv2A/asWN9fOGrzDfh33Pz0sYwBfKn9gD3zbr3G+sz/+PDL\nTusiQOeaOox+Wi2c7tSzN0G+qjENPdgAF2tqAvssxkfPTGYygCM0FbB3nTAPOXAD4CBdB+xQogVw\nqHZOZP7A+St7rjcBsJeuA7ajjfOWsQxgL50GbHFZqemTOzJmARyqx31gR5xKvCVXkgawo+4msH0r\n8qKCz/6TVUSAXXQ3ga378exNTxbf9ecbev15vtf8uvgAfKm7CWyjLW99sjiecLpc7yBRAMfrNGDr\ncWqenHWH1u9I+WzYWqTLQfMAB+kxYNO14YdZYKbMPDs08fVRi/P/avwCOEF3+8Am6xptGaoefpdi\nAZyvxwns2bT0536v8ihxr79ly7607R7egNNc+9MvvwE9//TLb0DPP/0ON+C2+p3AJuYngEQ9TmBf\nqF4JAdyEgC05bhAgQlPvyLz20drx+ICU1Wf+/V8/v0EAIe5fh8YDtosjrp0IwJccxPE36QK4IfvA\nAIgkYABEEjAAItkHttV0QONVh73Mj6g8+TaU8u9gn0seh+kGnPwgLO7syfd9/eN6vvuXPP+7fean\nELBNFs/jq55A5//cxXkI5z8O6xMhTn4Q5puP858D00+fflzPd//8123zj8+/+4snv26tWUJMUko5\n+apotdZrf23WN+DMB+Hy+77+ZM93/8z7fuHr1Gc34Pxf//szgSVZvx7v0PkPwvizrtp2LO7pyXf/\nVlfR7fz53/ndf8gEFsOzdrhiEenC7cX6p59/S64dwec//eQ9T4t/n2x9A/z6PyRgGSwdDBc9CPdZ\nRjv57t9q9rpk5fySHW8Pb4Bf/2dMo1vd5yjESw7luNuxWOcfPzLfjlzy04fTH/x7HoV4yRPv8htw\n+fbnngQMgEiWEAGIJGAARBIwACIJGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAEQSMAAiCRgAkQQM\ngEgCBkAkAQMgkoABEEnAAIgkYABEEjAAIgkYAJEEDIBIAgZAJAEDIJKAARBJwACIJGAARBIwACIJ\nGACRBAyASAIGQCQBAyCSgAEQScAAiCRgAET6fygCX480EVRFAAAAAElFTkSuQmCC\n", "output_type": "display_data"}], "prompt_number": 9, "cell_type": "code", "language": "python", "metadata": {}, "input": ["clf;\n", "hold('on');\n", "plot(X(1,:), X(2,:), '.');\n", "plot(m(1,:), m(2,:), 'k*');\n", "plot(med(1,:), med(2,:), 'ro');\n", "axis('tight');"]}, {"source": ["Color Image Denoising using 1D Median\n", "-------------------------------------\n", "A median filter can be used to denoise a color image, by applying it to each channel of the image.\n", "\n", "\n", "We load a color image, which is an array of size |[n,n,3]|."], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 10, "cell_type": "code", "language": "python", "metadata": {}, "input": ["name = 'flowers';\n", "options.nbdims = 3;\n", "n = 256;\n", "M0 = load_image(name, n, options);\n", "M0 = rescale(M0);"]}, {"source": ["We create a colored impulse noise by taking two Gaussians of different\n", "standard deviations.\n", "\n", "percent of strong Gaussian"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 11, "cell_type": "code", "language": "python", "metadata": {}, "input": ["rho = .4;"]}, {"source": ["mask of pixel corrupted by strong gaussian"], "metadata": {}, "cell_type": "markdown"}, {"collapsed": false, "outputs": [], "prompt_number": 12, "cell_type": "code", "language": "python", "metadata": {}, "input": ["mask = repmat(rand(n,n)