{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Image Processing with Wavelets", "==============================", "", "*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 overviews the use of wavelets for image approximation and", "denoising."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import introduction_5_wavelets_2d as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Wavelet Approximation", "---------------------", "Image approximation is obtained by thresholding wavelets coefficients.", "", "", "First we load an image $f \\in \\mathbb{R}^N$ of $N = n_0 \\times n_0$ pixels."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["name = 'cortex'", "n0 = 512", "f = load_image(name, n0)", "f = rescale(sum(f, 3))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display it."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["An orthogonal wavelet basis $ \\mathcal{B} = \\{ \\psi_{j,n}^k \\}_{j,n} $ of", "$\\mathbb{R}^N$ is composed of multiscale atoms parameterized by", "their scale $2^j$, position $2^j n \\in [0,1]^2$ and", "orentation $ k \\in \\{H,V,D\\}$.", "", "", "A forward wavelet transform computes the set of inner products", "$$ \\Psi f = \\{ \\langle f,\\psi_{j,n}^k\\rangle \\} \\in \\mathbb{R}^N $$", "where the wavelet atoms are defined as", "$$  \\psi_{j,n}^k(x) = \\psi^k\\left( \\frac{x - 2^j n}{2^j} \\right). $$", "", "", "Set the minimum scale for the transform."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Jmin = 0"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A short-cut for the wavelet transform $\\Psi$:"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Psi = lambda f: perform_wavelet_transf(f, Jmin, + 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A short-cut for the inverse wavelet transform $\\Psi^{-1} = \\Psi^*$:"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["PsiS = lambda fw: perform_wavelet_transf(fw, Jmin, -1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform the wavelet transform to compute $\\Psi f$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fW = Psi(f)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the transformed coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot_wavelet(fW)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["To perform non-linear image approximation, one remove the small amplitude", "coefficients. This is performed using a hard thresholding", "$$ H_T(f,\\mathcal{B}) = \\Psi^{-1} \\circ H_T \\circ \\Psi (f) =", " \\sum_{|\\langle f,\\psi_{j,n}^k\\rangle| > T} \\langle f,\\psi_{j,n}^k\\rangle \\psi_{j,n}^k. $$", "", "", "$T$ should be adapted to ensure a given number $M$ of non-zero", "coefficients, and then $f_M = H_T(f,\\mathcal{B})$ is the best $M$", "terms approximation of $f$ in $\\mathcal{B}$.", "", "", "Select a threshold."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["T = .5"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Shortcut for the thresholding operator $H_T$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Thresh = lambda fW, T: fW .* (abs(fW) >T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform hard thresholding of the coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fWT = Thresh(fW, T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Compute the ratio $M/N$ of non-zero coefficients."]}, {"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": ["Display the thresholded coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "plot_wavelet(fW)", "title('Original coefficients')", "subplot(1, 2, 2)", "plot_wavelet(fWT)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform reconstruction using the inverse wavelet transform $\\Psi^*$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = PsiS(fWT)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display approximation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f, 'Image', 1, 2, 1)", "imageplot(clamp(f1), strcat(['Approximation, SNR = ' num2str(snr(f, f1), 3) 'dB']), 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Number of coefficients for the approximation."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["M = n0^2/ 16"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Compute a threshold $T$ to keep only $M$ coefficients."]}, {"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": ["Perform hard thresholding."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fWT = Thresh(fW, T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Check the number of non-zero coefficients in |fWT|."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["disp(strcat(['      M = ' num2str(M)]))", "disp(strcat(['|fWT|_0 = ' num2str(sum(fWT(: )~ = 0))]))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Compute an approximation with an decreasing number of coefficients."]}, {"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": ["Orthognal Wavelet Denoising", "---------------------------", "Image denoising is obtained by thresholding noisy wavelets coefficielts.", "", "", "Here we consider a simple setting where we intentionnaly add some noise", "$w$ to a clean image $f$ to obtain $ y = f + w $.", "", "", "", "A denoiser computes an estimate $\\tilde f$ of $f$ from the", "observations $y$ alone. In the mathematical model, since $y$ is a random variable depending on", "$w$, so is $\\tilde f$. A mathematical evaluation of the efficiency of", "the denoiser is the average risk $E_w( \\|f-\\tilde f\\|^2 )$.", "", "", "", "Here we consider a single realization of the noise, so we replace the risk", "by the oracle error $ \\|f-\\tilde f\\|^2$.", "This allows us to bench the efficiency of the denoising methods by", "comparing the result to $f$. But you have to keep in mind that for real", "application, one does not have access to $f$.", "", "", "", "We consider a Gaussian white noise $w$ of variance $\\sigma^2$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["sigma = .1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate a noisy image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["y = f + randn(n0, n0)*sigma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f, 'Clean image', 1, 2, 1)", "imageplot(clamp(y), ['Noisy image, SNR = ' num2str(snr(f, y), 3) 'dB'], 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A denoising is obtained by thresholding the wavelet coefficients", "$$ \\tilde f = H_T(f,\\mathcal{B}). $$", "", "", "The asymptotically optimal threshold of Donoho and Johnstone is", "$T = \\sqrt{2 \\log(N)} \\sigma$. In practice, one observes that", "much better result are obtained using $T \\approx 3 \\sigma$.", "", "", "Compute the noisy wavelet coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fW = Psi(y)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the threshold value using the $3\\sigma$ heuristic."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["T = 3*sigma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform hard thresholding."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fWT = Thresh(fW, T)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the thresholded coefficients."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["subplot(1, 2, 1)", "plot_wavelet(fW)", "title('Original coefficients')", "subplot(1, 2, 2)", "plot_wavelet(fWT)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Perform reconstruction."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = PsiS(fWT)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display denoising."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(clamp(y), 'Noisy image', 1, 2, 1)", "imageplot(clamp(f1), strcat(['Denoising, SNR = ' num2str(snr(f, f1), 3) 'dB']), 1, 2, 2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 4__", "", "Try to optimize the value of the threshold $T$ to get the best possible", "denoising 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": ["Translation Invariant Denoising", "-------------------------------", "The quality of orthogonal denoising is improved by adding translation", "invariance. This corresponds to denoising translated copies of the image.", "", "", "The translation of an image is $(\\theta_\\tau f)(x) = f(x-\\tau)$, where", "we use periodic boundary conditions.", "", "", "Given a set $ \\Omega \\subset \\mathbb{R}^2 $, the $\\Omega$-translation invariant", "denoising is defined as:", "$$ \\tilde f = \\frac{1}{\\Omega}\\sum_{\\tau \\in \\Omega} \\theta_{-\\tau} \\left( H_T( \\theta_\\tau y, \\mathcal{B} ) \\right). $$", "", "", "Here we consider translation of integer pixels in", "$\\{0,\\ldots,\\tau_{\\max}-1\\}$. The number of translations is thus", "$ \\tau_{\\max}^2$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["tau_max = 8"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Generate a set of translation vectors $\\Omega = \\{ \\tau_i = (X_i,Y_i) \\}_i$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[Y, X] = meshgrid(0: tau_max-1, 0: tau_max-1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["A \"trick\" to compute the full denoising image after all translations is to", "initialize $\\tilde f = 0$, and then accumulate each denoising with translate $\\tau_i$ in the", "following way:", "$$ \\tilde f \\longleftarrow \\frac{i-1}{i} \\tilde f + \\frac{1}{i} \\theta_{-\\tau_i} \\left( H_T( \\theta_{\\tau_i} y, \\mathcal{B} ) \\right) $$", "", "", "", "Initialize the denoised image $\\tilde f$ as 0."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = zeros(n0, n0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Initialize the translation index."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["i = 1"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Translate the image to obtain $\\theta_{\\tau_i}(f)$", "for $\\tau_i = (X_i,Y_i)$, with periodic boundary conditions."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fTrans = circshift(y, [X(i) Y(i)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Denoise this translated image, to obtain $H_T(\\theta_{\\tau_i} f,\\mathcal{B})$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fTrans = PsiS(Thresh(Psi(fTrans) , T))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Translate back."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["fTrans = circshift(fTrans, -[X(i) Y(i)])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Accumulate the result."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = (i-1)/ i*f1 + fTrans/ i"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 5__", "", "Compute the full denoising by cycling through the $i$ indices."]}, {"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": ["__Exercise 6__", "", "Determine the optimal threshold $T$ for this translation invariant", "denoising."]}, {"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__", "", "Test on other images."]}, {"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."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}