{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Data Dependent 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": ["The simplest noise model is Gaussian additive noise, where the variance", "of the pixel value is independent of the mean (which is the value we look", "for).", "", "", "Unfortunately, most real-life data corresponds to noise model that are", "much more complicated, and in particular the variance of the noise often", "depends on the parameter of interest (for instance the mean)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import denoisingwav_5_data_dependent as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Poisson Noise", "-------------", "A Poisson model assume that each pixel $x$", "of an image $f(x)$ is drawn from", "a Poisson distribution of parameter $\\lambda=f_0(x)$, where", "$f_0$ is the clean intensity image to recover.", "", "", "$$ \\PP(f(x)=k)=\\lambda^k e^{-\\lambda}/k! $$", "", "", "Display the Poisson distribution for several value of $\\lambda$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["lambda = [4 10 20]", "[k, Lambda] = meshgrid(1: 50, lambda)", "P = Lambda.^k .* exp(-Lambda)./ factorial(k)", "h = plot(P'); axis('tight')", "if using_matlab()", "    set(h, 'LineWidth', 2)", "", "legend('\\lambda = 2', '\\lambda = 10', '\\lambda = 20')", "set_label('k', 'P(k)')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["This model corresponds to a photon count, where $\\lambda$", "is proportional to the number of photons that hits the receptor during", "the exposition time. This is useful to model medical imaging (confocal", "microscopy), TEP and SPECT tomography, and digital camera noises.", "", "", "The goal of denoising is to retrieve the value of", "$f_0(x)=\\lambda(x)$, the true", "image value, which corresponds to a clean image.", "", "", "Note that $\\lambda$ is the mean of the distribution, but is also its variance,", "so that the noise intensity perturbating the image pixel $f(x)$", "is proportional to $f_0(x)$.", "", "", "We load a clean, unquantized image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256", "name = 'lena'", "f0u = rescale(load_image(name, n))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Quantize the values to given $\\lambda$ range."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["lmin = 1", "lmax = 40", "f0 = floor(rescale(f0u, lmin, lmax))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate the noisy image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f = poissrnd(f0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Intensity map f0', 1, 2, 1)", "imageplot(f,  'Observed noisy image f', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the difference, which shows that the noise level depends on the", "intensity. The noise is larger in bright areas."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Intensity map f0', 1, 2, 1)", "imageplot(f-f0,  'f-f0', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Display noisy image contaminated by Poisson noise of varying range."]}, {"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": ["Parameters for the wavelet transform."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Jmin = 4", "options.ti = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A translation invariance wavelet transform denoising computes", "an estimate $\\tilde f$ of the clean image", "$f_0$ as", "$$ \\tilde f = \\Psi^+ \\circ S_T \\circ \\Psi f $$", "where $\\Psi$ is the wavelet transform and $S_T$ the hard", "thresholding operator using a well chosen threshold $T$."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Perform translation invariant wavelet hard thresholding directly on the", "Poisson noisy image $f$. Check for an optimal threshold that maximize", "the SNR. Record the optimal result in |fPoisson|.", "isplay", "est result"]}, {"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": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Original image', 1, 2, 1)", "imageplot(clamp(fPoisson, min(f0(: )), max(f0(: ))), strcat(['Denoised, SNR = ' num2str(snr(f0, fPoisson), 4)]), 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Variance Stabilization Poisson Denoising", "----------------------------------------", "A variance stabilization transform (VST)", "is a mapping $\\phi$ applied to a noisy image $f$", "so that the distribution of each pixel $\\phi(f(x))$ is", "approximately Gaussian.", "", "", "The Anscombe variance stabilization transform (VST) is given by", "the non-linear mapping.", "", "", "$$\\phi(a) = 2 \\sqrt{a+3/8}$$", "", "", "It transforms a Poisson distribution $P(\\lambda)$ to a", "distribution with approximate", "variance 1, whatever $\\lambda$ is, for $\\lambda$ large enough", "(say $\\lambda>20$).", "", "", "Other VST have been proposed, for instance the Freeman & Tukey VST, given", "by", "", "", "$$\\phi(a) = \\sqrt{a+1}+\\sqrt{a}$$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Display the estimated variance of a Poisson distribution for various", "$\\lambda$ (e.g. 10000 realization for each $\\lambda$) and display", "the variance of a stabilized distribution (here the green", "curve corresponds to 'Freeman & Tukey' and the blue curve to 'Anscombe'.", "lot", "bplot(2,1,1);", "= plot(v); axis('tight');", "using_matlab()", " set(h, 'LineWidth', 2);", "d", "tle('Poisson variance');", "t_label('\\lambda', 'Variance');", "bplot(2,1,2);"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo3()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["To perform denosing, one applies the VST, performs", "wavelet thresholding as if the data was corrupted by an additive Gaussian", "noise of variance $\\sigma=1$, and then applies the inverse VST.", "", "", "This corresponds to computing an estimate $\\tilde f$ of the clean image", "$f_0$ as", "$$ \\tilde f = \\phi^{-1} \\pa{ \\Psi^+ S_T \\circ \\Psi \\circ \\phi(f) } $$", "where $\\Psi$ is the wavelet transform and $S_T$ the hard thresholding."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Perform translation invariance wavelet hard thresholding on the", "variance stabilized image. Use for instance the Anscombe VST.", "Check for an optimal threshold that maximize", "the SNR. Record the optimal result in |fVST|.", "isplay", "est result"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo4()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(clamp(fPoisson, min(f0(: )), max(f0(: ))), ...", "        strcat(['Un-stabilized, SNR = ' num2str(snr(f0, fPoisson), 4)]), 1, 2, 1)", "imageplot(clamp(fVST, min(f0(: )), max(f0(: ))), ...", "        strcat(['Stabilized, SNR = ' num2str(snr(f0, fVST), 4)]), 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["There is no close form solution for the Freeman VST.", "For each denoised stabilized pixel value $b \\in \\RR$,", "one should use a Newton algorithm to solve the equation", "$\\phi(a)=b$ and recovered the denoised value $a \\in \\RR$.", "This reads", "$$ a_{k+1} = a_k - \\frac{\\phi(a_k)-y}{\\phi'(a_k)} $$"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 5__", "", "Perform VST denoising using the Freeman VST."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo5()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Multiplicative Noise", "--------------------", "A multiplicative noise corresponds to a noise model where the clean image", "$f_0$ is multiplied by the noise $W$ to obtained the noisy image $f(x)=W(x) f_0(x)$.", "", "", "This model is useful for data acquired with an active acquisition device,", "for instance SAR imaging and ultra-sound imaging.", "", "", "The distribution of $f(x)$ thus has mean $f_0(x)$ and variance", "proportional to $f_0(x) \\sigma$ where $\\sigma$ is the variance of $W$.", "", "", "Load a clean image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256", "name = 'boat'", "f0 = rescale(load_image(name, n), 1e-2, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A classical model for the multiplier noise $W$, that is used for instance in", "SAR imaging, is to assume that the noise has a Gamma law", "of mean 1 and variance parameterized by the noise level $L$.", "$$ P(W=x) \\sim x^{K-1} e^{ -x/\\theta } $$", "", "", "where the mean of $P$ is $s=K \\theta$ and the variance is", "$\\sigma^2=K \\theta^2$.", "", "", "This corresponds to a acquisition device which is", "averaging $K$ measures with 1-sided exponential distribution."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["K = 4", "sigma = 1/ sqrt(K)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate the random multiplier."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = gamrnd(1/ sigma^2, sigma^2, n, n)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate the noisy signal."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f = f0.*W"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Intensity map f0', 1, 2, 1)", "imageplot(f,  'Observed noisy image f', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the difference, which shows that the noise level depends on the", "intensity. The noise is larger in bright areas."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Intensity map f0', 1, 2, 1)", "imageplot(f-f0,  'f-f0', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 6__", "", "Generate several noisy images for several noise levels."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo6()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 7__", "", "Perform translation invariance wavelet hard thresholding directly on the", "noisy image $f=f_0 W$. Check for an optimal threshold that maximize", "the SNR. Record the optimal result in |fMult|.", "isplay", "est result"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo7()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0, 'Original image', 1, 2, 1)", "imageplot(clamp(fMult, min(f0(: )), max(f0(: ))), strcat(['Denoised, SNR = ' num2str(snr(f0, fMult), 4)]), 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Variance Stabilization for Multiplicative Noise", "-----------------------------------------------", "An approximate variance stabilization transform consist in taking the", "log. The $\\log(f)-a$ image is then equal to $\\log(f_0)$ contaminated by $log(W)-a$ which is not", "too far from a centered Gaussian distribution if the variance of $W$ is", "not too large.", "", "", "The value of $a$ should be chosen as the mean value of the", "random variable $\\log(f)$ so that $\\log(W)-a$ has zero mean.", "There exists close form for the value of $a$ as a function of $\\sigma=1/\\sqrt{K}$,", "which we use here, where $\\psi$ is the polygamma function.", "", "", "_Important:_ Scilab user should replace |psi| by |dlgamma|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["a = psi(K) - log(K)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display stabililized image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f, 'f', 1, 2, 1)", "imageplot(clamp(log(f)-a, -2, 2), 'log(f)', 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Distribution of the noise multiplier and of the log."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(2, 1, 1)", "hist(W(: ), 100)", "axis('tight')", "subplot(2, 1, 2)", "hist(log(W(: ))-a, 100)", "axis('tight')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 8__", "", "Perform translation invariance wavelet hard thresholding on the", "variance stabilized image using the log.", "Check for an optimal threshold that maximize", "the SNR. Record the optimal result in |fVST|.", "isplay", "est result"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["solutions.exo8()"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["## Insert your code here."]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(clamp(fMult, min(f0(: )), max(f0(: ))), strcat(['Un-stabilized, SNR = ' num2str(snr(f0, fMult), 4)]), 1, 2, 1)", "imageplot(clamp(fVST, min(f0(: )), max(f0(: ))), strcat(['Stabilized, SNR = ' num2str(snr(f0, fVST), 4)]), 1, 2, 2)"]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}