{"worksheets": [{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["Shape From Shading", "==================", "", "*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 tour explores the resolution of the shape from shading inverse", "problem using Fast Marching."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["from __future__ import division", "import nt_toolbox as nt", "from nt_solutions import graphics_7_shape_shading as solutions", "%matplotlib inline", "%load_ext autoreload", "%autoreload 2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Forward Image Formation", "-----------------------", "We consider here a simplified imaging system where the camera performs an", "orthogonal projection and the lighting is orthogonal to the camera focal", "plane and is at infinity. We also assume that the surface is perfectly", "diffuse (Lambertian).", "", "", "For more information about the shape from shading problem, including more", "complicated settings, one can read the review paper", "", "", "_Shape From Shading_, Emmanuel Prados, Olivier Faugeras", "Handbook of Mathematical Models in Computer Vision,", "Springer, page 375-388, 2006.", "", "", "The shape use in this tour is taken from", "", "", "_Shape from Shading: A Survey_, Ruo Zhang, Ping-Sing Tsai, James Edwin, Cryer", "and Mubarak Shah, IEEE Trans. on PAMI, Vol. 21, No. 8, 1999.", "", "", "", "Load a surface as a 2D height field, i.e. it is defined using an image $f(x,y)$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["name = 'mozart'", "f = load_image(name)", "n = size(f, 1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Height $h>0$ of the surface."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["h = n*.4"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Final rescale surface $(x,y,f(x,y)) \\in \\RR^3$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f = rescale(f, 0, h)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(f)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display as a 3-D surface."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["surf(f)", "colormap(gray(256))", "shading interp", "axis equal", "view(110, 45)", "axis('off')", "camlight"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the normal to the surface $ (x,y,f(x,y)) $.", "The tangent vectors are $ (1,0,\\pd{f}{x}(x,y)) $ and $", "(0,1,\\pd{f}{y}(x,y)) $", "and the normal is thus", "$$ N(x,y) = \\frac{N_0(x,y)}{\\norm{N_0(x,y)}} \\in \\RR^3", "  \\qwhereq N_0(x,y) = \\pa{\\pd{f}{x}(x,y), \\pd{f}{y}(x,y), 1}.   $$", "", "", "Compute $N_0$, the un-normalized normal, using centered finite differences."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["options.order = 2", "N0 = cat(3, -grad(f, options), ones(n))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the normalized normal $N$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["s = sqrt(sum(N0.^2, 3))", "N = N0 ./ repmat(s, [1 1 3])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the normal map as a color image."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(N)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We compute the shaded image obtained by an infinite light coming from", "a given direction $d \\in \\RR^3$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["d = [0 0 1]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We use a Labertian reflectance model to determine the luminance", "$$ L(x,y) = (\\dotp{N(x,y)}{ d })_+", "      \\qwhereq (\\al)_+ = \\max(\\al,0). $$"]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["L = max(0, sum(N .* repmat(reshape(d, [1 1 3]), [n n 1]), 3))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the lit surface."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["vmin = .3", "", "imageplot(max(L, vmin))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 1__", "", "Display the surface with several light directions $d \\in \\RR^3$."]}, {"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": ["Eikonal Shape From Shading", "--------------------------", "For a vertical ligthing direction $d=(0,0,1)$, the forward imaging", "operator reads", "$$ L(x,y) = \\frac{1}{ \\sqrt{ \\norm{ \\nabla f(x,y) }^2+1 } }. $$", "", "", "Compute the luminance map for $d=(0,0,1)$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["d = [0 0 1]", "L = sum(N .* repmat(reshape(d, [1 1 3]), [n n 1]), 3)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The height field $f$ can thus, in theory, be recovered by solving the", "following Eikonal equation:", "$$ \\norm{ \\nabla f(x,y) } = W(x,y)", "      \\qwhereq W(x,y) = \\sqrt{ 1/L(x,y)^2 - 1 } $$", "With additional boundary conditon", "$$ \\forall i \\in I, \\quad f( p_i) = f_i $$", "For a set of well chosen grid points $ p_i \\in \\RR^2  $.", "", "", "The issue is that this equation is ill-posed (several solution) due to", "the singular points $ (q_j)_{j \\in J} $ where $ L(q_j) $", "is close to 1, so that $W(q_j)$ is close to 0.", "A way to regularize the inversion is thus", "to choose the boundary condition points $p_i$ to be the singular points", "$q_j$.", "", "", "Set a tolerance $\\epsilon>0$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["epsilon = 1e-9"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We show in white the location of singular points $q_j$,", "that satisfies $L(q_j)>1-\\epsilon$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(L >1-epsilon)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Compute the Eikonal speed $W$ (right hand side of the equation)."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = sqrt(1./ L.^2-1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["To avoid too much ill-posedness, we threshold it."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["W = max(W, epsilon)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display the speed."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["imageplot(min(W, 3))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We use here a single point $ p $.", "Select the tip of the nose as base point $p$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["p = [140; 120]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Solve the Eikonal equation, assuming $f(p)=0$."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["[f1, S] = perform_fast_marching(1./ W, p)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Rescale the height."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["f1 = -f1*n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Display as a 3D surface."]}, {"cell_type": "code", "language": "python", "metadata": {}, "outputs": [], "collapsed": false, "input": ["hold on", "surf(f1)", "h = plot3(p(2), p(1), f1(p(1), p(2)), 'r.')", "set(h, 'MarkerSize', 30)", "colormap(gray(256))", "shading interp", "axis('equal')", "view(110, 45)", "axis('off')", "camlight"]}, {"cell_type": "markdown", "metadata": {}, "source": ["__Exercise 2__", "", "Try to reconstruct the image starting from other base points.", "What do you observe ?"]}, {"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": ["__Exercise 3__", "", "Try to improve the quality of the reconstruction by selecting several", "points, and imposing their height."]}, {"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."]}]}], "nbformat": 3, "metadata": {"name": ""}, "nbformat_minor": 0}