{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Signal and Image Noise Models", "=============================", "", "*Important:* Please read the [installation page](http://gpeyre.github.io/numerical-tours/installation_python/) for details about how to install the toolboxes.", "$\\newcommand{\\dotp}[2]{\\langle #1, #2 \\rangle}$", "$\\newcommand{\\enscond}[2]{\\lbrace #1, #2 \\rbrace}$", "$\\newcommand{\\pd}[2]{ \\frac{ \\partial #1}{\\partial #2} }$", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$", "$\\newcommand{\\umax}[1]{\\underset{#1}{\\max}\\;}$", "$\\newcommand{\\umin}[1]{\\underset{#1}{\\min}\\;}$", "$\\newcommand{\\uargmin}[1]{\\underset{#1}{argmin}\\;}$", "$\\newcommand{\\norm}[1]{\\|#1\\|}$", "$\\newcommand{\\abs}[1]{\\left|#1\\right|}$", "$\\newcommand{\\choice}[1]{ \\left\\{  \\begin{array}{l} #1 \\end{array} \\right. }$", "$\\newcommand{\\pa}[1]{\\left(#1\\right)}$", "$\\newcommand{\\diag}[1]{{diag}\\left( #1 \\right)}$", "$\\newcommand{\\qandq}{\\quad\\text{and}\\quad}$", "$\\newcommand{\\qwhereq}{\\quad\\text{where}\\quad}$", "$\\newcommand{\\qifq}{ \\quad \\text{if} \\quad }$", "$\\newcommand{\\qarrq}{ \\quad \\Longrightarrow \\quad }$", "$\\newcommand{\\ZZ}{\\mathbb{Z}}$", "$\\newcommand{\\CC}{\\mathbb{C}}$", "$\\newcommand{\\RR}{\\mathbb{R}}$", "$\\newcommand{\\EE}{\\mathbb{E}}$", "$\\newcommand{\\Zz}{\\mathcal{Z}}$", "$\\newcommand{\\Ww}{\\mathcal{W}}$", "$\\newcommand{\\Vv}{\\mathcal{V}}$", "$\\newcommand{\\Nn}{\\mathcal{N}}$", "$\\newcommand{\\NN}{\\mathcal{N}}$", "$\\newcommand{\\Hh}{\\mathcal{H}}$", "$\\newcommand{\\Bb}{\\mathcal{B}}$", "$\\newcommand{\\Ee}{\\mathcal{E}}$", "$\\newcommand{\\Cc}{\\mathcal{C}}$", "$\\newcommand{\\Gg}{\\mathcal{G}}$", "$\\newcommand{\\Ss}{\\mathcal{S}}$", "$\\newcommand{\\Pp}{\\mathcal{P}}$", "$\\newcommand{\\Ff}{\\mathcal{F}}$", "$\\newcommand{\\Xx}{\\mathcal{X}}$", "$\\newcommand{\\Mm}{\\mathcal{M}}$", "$\\newcommand{\\Ii}{\\mathcal{I}}$", "$\\newcommand{\\Dd}{\\mathcal{D}}$", "$\\newcommand{\\Ll}{\\mathcal{L}}$", "$\\newcommand{\\Tt}{\\mathcal{T}}$", "$\\newcommand{\\si}{\\sigma}$", "$\\newcommand{\\al}{\\alpha}$", "$\\newcommand{\\la}{\\lambda}$", "$\\newcommand{\\ga}{\\gamma}$", "$\\newcommand{\\Ga}{\\Gamma}$", "$\\newcommand{\\La}{\\Lambda}$", "$\\newcommand{\\si}{\\sigma}$", "$\\newcommand{\\Si}{\\Sigma}$", "$\\newcommand{\\be}{\\beta}$", "$\\newcommand{\\de}{\\delta}$", "$\\newcommand{\\De}{\\Delta}$", "$\\newcommand{\\phi}{\\varphi}$", "$\\newcommand{\\th}{\\theta}$", "$\\newcommand{\\om}{\\omega}$", "$\\newcommand{\\Om}{\\Omega}$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This numerical tour show several models for signal and image noise.", "It shows how to estimate the noise level for a Gaussian additive noise on a natural image.", "It also shows the relevance of thresholding to remove Gaussian noise", "contaminating sparse data."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import denoisingsimp_1_noise_models as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Additive Gaussian Noise Model ", "------------------------------", "The simplest noise model consist in adding a realization of a zero mean", "random vector to a clean signal or image.", "", "", "Load a clean image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["N = 128", "name = 'boat'", "M0 = load_image(name, 256)", "M0 = rescale(crop(M0, N))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Load a clean signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 1024", "name = 'piece-regular'", "f0 = rescale(load_signal(name, n))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The simplest noise model is Gaussian white noise.", "Here we generate a noisy signal or image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma = .1", "M = M0 + randn(N, N)*sigma", "f = f0 + randn(n, 1)*sigma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the signals."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(3, 1, 1)", "plot(f0); axis([1 n 0 1])", "title('Clean signal')", "subplot(3, 1, 2)", "plot(f-f0); axis([1 n -3*sigma 3*sigma])", "title('Noise')", "subplot(3, 1, 3)", "plot(f); axis([1 n 0 1])", "title('Noisy signal')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the images."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(M0, 'Clean image', 1, 3, 1)", "imageplot(M-M0, 'Noise', 1, 3, 2)", "imageplot(clamp(M), 'Noisy image', 1, 3, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the statistics of the noise"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["nbins = 51", "[h, t] = hist(M(: )-M0(: ), nbins); h = h/ sum(h)", "subplot(3, 1, 2)", "bar(t, h)", "axis([-sigma*5 sigma*5 0 max(h)*1.01])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Additive Uniform Noise", "----------------------", "A slightly different kind of noise is uniform in a given interval.", "", "", "Generate noisy data with uniform", "noise distribution in |[-a,a]|, with |a| chosen", "so that the variance is sigma."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = sqrt(3)*sigma", "M = M0 + 2*(rand(N, N)-.5)*a", "f = f0 + 2*(rand(n, 1)-.5)*a"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the signals."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(3, 1, 1)", "plot(f0); axis([1 n 0 1])", "title('Clean signal')", "subplot(3, 1, 2)", "plot(f-f0); axis([1 n -3*sigma 3*sigma])", "title('Noise')", "subplot(3, 1, 3)", "plot(f); axis([1 n 0 1])", "title('Noisy signal')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the images."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(M0, 'Clean image', 1, 3, 1)", "imageplot(M-M0, 'Noise', 1, 3, 2)", "imageplot(clamp(M), 'Noisy image', 1, 3, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the statistics of the noise"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["nbins = 51", "[h, t] = hist(M(: )-M0(: ), nbins); h = h/ sum(h)", "subplot(3, 1, 2)", "bar(t, h)", "axis([-sigma*5 sigma*5 0 max(h)*1.01])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Impulse Noise", "-------------", "A very different noise model consist in sparse impulsions, generate by a", "random distribution with slowly decaying probability.", "", "", "Generate noisy image with exponential distribution,", "with variance |sigma|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = log(rand(N, N)).*sign(randn(N, N))", "W = W/ std(W(: ))*sigma", "M = M0 + W"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate noisy signal with exponential distribution,", "with variance |sigma|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = log(rand(n, 1)).*sign(randn(n, 1))", "W = W/ std(W(: ))*sigma", "f = f0 + W"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the signals."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(3, 1, 1)", "plot(f0); axis([1 n 0 1])", "title('Clean signal')", "subplot(3, 1, 2)", "plot(f-f0); axis([1 n -3*sigma 3*sigma])", "title('Noise')", "subplot(3, 1, 3)", "plot(f); axis([1 n 0 1])", "title('Noisy signal')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the images."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(M0, 'Clean image', 1, 3, 1)", "imageplot(M-M0, 'Noise', 1, 3, 2)", "imageplot(clamp(M), 'Noisy image', 1, 3, 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the statistics of the noise"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["nbins = 51", "[h, t] = hist(M(: )-M0(: ), nbins); h = h/ sum(h)", "subplot(3, 1, 2)", "bar(t, h)", "axis([-sigma*5 sigma*5 0 max(h)*1.01])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Thresholding Estimator and Sparsity", "-----------------------------------", "The idea of non-linear denoising is to use an orthogonal basis in which", "the coefficients |x| of the signal or image |M0| is sparse (a few large", "coefficients). In this case, the noisy coefficients |x| of the noisy", "data |M| (perturbated with Gaussian noise) are |x0+noise| where |noise|", "is Gaussian. A thresholding set to 0 the noise coefficients that are", "below |T|. The threshold level |T| should be chosen judiciously to be", "just above the noise level.", "", "", "First we generate a spiky signal.", "", "dimension"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 4096"]}, {"cell_type": "markdown", "metadata": {}, "source": ["probability of spiking"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["rho = .05"]}, {"cell_type": "markdown", "metadata": {}, "source": ["location of the spike"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["x0 = rand(n, 1) <rho"]}, {"cell_type": "markdown", "metadata": {}, "source": ["random amplitude in [-1 1]"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["x0 = 2 * x0 .* (rand(n, 1)-.5)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We add some gaussian noise"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma = .1", "x = x0 + randn(size(x0))*sigma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(2, 1, 1)", "plot(x0); axis([1 n -1 1])", "set_graphic_sizes([], 20)", "title('Original signal')", "subplot(2, 1, 2)", "plot(x); axis([1 n -1 1])", "set_graphic_sizes([], 20)", "title('Noisy signal')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "What is the optimal threshold |T| to remove as much as possible of noise", "? Try several values of |T|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo1()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["In order to be optimal without knowing in advance the amplitude of the", "coefficients of |x0|, one needs to set |T| just above the noise level.", "This means that |T| should be roughly equal to the maximum value of a", "Gaussian white noise of size |n|."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "The theory predicts that the maximum of |n| Gaussian variable of variance |sigma^2|", "is smaller than |sqrt(2*log(n))| with large probability (that tends to 1", "when |n| increases). This is also a sharp result. Check this numerically", "by computing with Monte Carlo sampling the maximum with |n| increasing", "(in power of 2). Check also the deviation of the maximum when you", "perform several trial with |n| fixed."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo2()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Estimating the noise level", "--------------------------", "In practice, the noise level |sigma| is unknown.", "For additive Gaussian noise,", "a good estimator is given by the median of the wavelet coefficients at", "the finer scale. An even simple estimator is given by the normalized", "derivate along X or Y direction", "", "", "Load a clean image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256", "M0 = rescale(load_image('boat', n))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate a noisy image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma = 0.06", "M = M0 + randn(n, n)*sigma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["First we extract the high frequency residual."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["H = M", "H = (H(1: n-1, : ) - H(2: n, : ))'/ sqrt(2)", "H = (H(1: n-1, : ) - H(2: n, : ))'/ sqrt(2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(clamp(M), 'Noisy image', 1, 2, 1)", "imageplot(H, 'Derivative image', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Histograms."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[h, t] = hist(H(: ), 100)", "h = h/ sum(h)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display histogram."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["bar(t, h)", "axis([-.5 .5 0 max(h)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The mad estimator (median of median) must be rescaled so that", "it gives the correct variance for gaussian noise."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma_est = mad(H(: ), 1)/ 0.6745", "disp(strcat(['Estimated noise level = ' num2str(sigma_est), ', true = ' num2str(sigma)]))"]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}