Shape From Shading¶
$\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}$
This tour explores the resolution of the shape from shading inverse problem using Fast Marching.
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('toolbox_graph')
addpath('solutions/graphics_7_shape_shading')
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)$.
name = 'mozart';
f = load_image(name);
n = size(f,1);
Height $h>0$ of the surface.
h = n*.4;
Final rescale surface $(x,y,f(x,y)) \in \RR^3$.
f = rescale(f,0,h);
Display the image.
clf;
imageplot(f);
Display as a 3-D surface.
clf;
surf(f);
colormap(gray(256));
shading interp;
axis equal;
view(110,45);
axis('off');
camlight;
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.
options.order = 2;
N0 = cat(3, -grad(f,options), ones(n));
Compute the normalized normal $N$.
s = sqrt( sum(N0.^2,3) );
N = N0 ./ repmat( s, [1 1 3] );
Display the normal map as a color image.
clf;
imageplot(N);
We compute the shaded image obtained by an infinite light coming from a given direction $d \in \RR^3$.
d = [0 0 1];
We use a Labertian reflectance model to determine the luminance $$ L(x,y) = (\dotp{N(x,y)}{ d })_+ \qwhereq (\al)_+ = \max(\al,0). $$
L = max(0, sum( N .* repmat(reshape(d,[1 1 3]), [n n 1]),3 ) );
Display the lit surface.
vmin = .3;
clf;
imageplot(max(L,vmin));
Exercise 1
Display the surface with several light directions $d \in \RR^3$.
exo1()