{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Inpainting using Variational Regularization", "===========================================", "", "*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 explores the use of variational energies (Sobolev,", "total variation) to regularize the image inpaiting problem.", "", "", "Here we consider inpainting of damaged observation without noise."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import inverse_4_inpainting_variational as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Missing Pixels and Inpainting", "-----------------------------", "Inpainting corresponds to filling holes in images.", "", "", "First we load the image $f_0 \\in \\RR^N$ of $N=n\\times n$ to be inpainted."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["name = 'cameraman'", "n = 256", "f0 = rescale(load_image(name, n))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the original image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Ratio of removed pixels."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["rho = .7"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Construct a random mask $\\Ga = \\chi_{\\Om}$", "so that $\\Ga_i=0$ for removed pixels $i \\notin \\Om$,", "and $\\Ga_i=1$ for kept pixels $i \\in \\Om$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Gamma = rand(n) >rho"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We create the masking operator $\\Phi$ which is a diagonal operator:", "$$ (\\Phi f)_i = \\Ga_i f_i $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Phi = lambda f: f.*Gamma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the damaged observation $y=\\Phi(f_0)$ (no noise is added)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["y = Phi(f0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the observations."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(y)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Sobolev Impainting", "------------------", "We solve the inpainting problem by minimzing the Sobolev norm of the", "image under the constraint of matching the observation", "$$ f^\\star = \\uargmin{ \\Phi(f) = y } E(f) = \\norm{\\nabla f}^2 $$", "where $\\nabla$ is a finite difference approximation of the gradient.", "", "", "It can be shown that the solution to this problem is an harmonic function", "with prescribed boundary condition", "$$ \\forall i \\notin \\Om, \\quad (\\Delta f^\\star)_i=0", "      \\qandq \\forall i \\in \\Om, \\quad f^\\star_i = y_i. $$", "", "", "This problem requires the constrained minimization of a smooth function,", "it can thus be solved using a projected gradient descent", "$$ f^{(\\ell+1)} = \\Pi \\pa{ f^{(\\ell)} + \\tau \\Delta(f^{(\\ell)}) }$$", "where $ \\Pi $ is the orthogonal projector on the constraint $y=\\Phi f$", "$$ (\\Pi f)_i =", "  \\choice{", "      y_i \\qifq i \\in \\Om, \\\\", "      f_i \\qifq i \\notin \\Om, \\\\", "  }", " $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Pi = lambda f: f.*(1-Gamma) + y.*Gamma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Here $ \\Delta = -\\nabla^* \\circ \\nabla = \\text{div} \\circ \\nabla $ is the", "gradient of the Sobolev energy $ E $."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Delta = lambda f: div(grad(f))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["For convergence, the gradient descent step size should satisfy:", "$$ \\tau<\\frac{2}{\\norm{\\Delta}}=\\frac{1}{4} $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["tau = .8/ 4"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Perform the projected gradient descent.", "Record in a variable |E| the evolution of the Sobolev energy $E$."]}, {"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 decay of the energy $E(f^{(\\ell)})$", "with the iterations."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot(E); axis('tight')", "set_label('Iteration #', 'E')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the result."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f, strcat(['Inpainted, SNR = ' num2str(snr(f0, f), 3) 'dB']))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Inpainting with TV Regularization", "---------------------------------", "A non-linear prior replaces the Sobolev energy by the TV norm, that", "tends to better reconstruct edges. Here we use a smoothed TV norm to", "avoid convergence issue with gradient descent algorithms.", "", "", "The smoothed TV norm reads:", "$$ J_\\epsilon(f) = \\sum_x \\sqrt{\\norm{ \\nabla f(x) }^2+\\epsilon^2} $$", "", "", "We use a projected gradient descent to solve this problem", "$$ f^{(\\ell+1)} = \\Pi \\pa{ f^{(\\ell)} + \\tau G_\\epsilon(f^{(\\ell)}) }$$", "where $ G_\\epsilon $ is the gradient of $J_\\epsilon$, that is defined", "as", "$$ G_\\epsilon(f) = -\\text{div} N_\\epsilon( \\nabla f ) $$", "where $ N_\\epsilon $ is the following normalization operator", "$$ N_\\epsilon(u)_i = \\frac{u_i}{ \\sqrt{\\norm{u_i}^2 + \\epsilon^2} } $$", "that is applied to any vector field $u=(u_i)_i \\in \\RR^{N \\times 2}  $", "for $u_i \\in \\RR^2$.", "", "", "Regularization parameter $\\epsilon$ for the TV norm"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["epsilon = 1e-2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Define the normalization operator."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Amplitude = lambda u: sqrt(sum(u.^2, 3) + epsilon^2)", "Neps = lambda u: u./ repmat(Amplitude(u), [1 1 2])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The step size $\\tau$, should satisfy", "$$ \\tau<\\frac{\\epsilon}{4}. $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["tau = .9*epsilon/ 4"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Define the gradient of $J$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["G = lambda f: -div(Neps(grad(f)))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Perform the projected gradient descent.", "Record in a variable |J| the evolution of the TV energy $J_\\epsilon$."]}, {"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 the result."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(clamp(f), strcat(['SNR = ' num2str(snr(f0, f), 3) 'dB']))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the evolution of the TV norm $J_\\epsilon$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["plot(J)", "axis('tight')", "set_label('Iteration #', 'J_\\epsilon')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Inpainting with non-random mask", "-------------------------------", "Inpainting can be used to remove objects in pictures.", "", "", "Load an image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["n = 256", "f0 = load_image('parrot', n)", "f0 = rescale(sum(f0, 3))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display it."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Load the mask."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Gamma = load_image('parrot-mask', n)", "Gamma = double(rescale(Gamma) >.5)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Masking operator $\\Phi$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["Phi = lambda f: f.*Gamma"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Observation $y=\\Phi(f_0)$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["y = Phi(f0)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display it."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(y)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 3__", "", "Perform Sobolev inpainting on this image."]}, {"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": ["__Exercise 4__", "", "Try other methods to solve this inpainting problem.", "You can for instance have a look on the numerical on sparsity for", "deconvolution and inpainting."]}, {"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."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}