```SHAPE FROM SHADING
Mica Arie-Nachimson
May 2011


We can’t reconstruct shape from one shadow…
Image from www.moolf.com

Variable levels of darkness

Gives a cue for the actual 3D shape
There is a relation between intensity and shape

Images from S. Narasimhan, Carnegie Mellon; www.erco.com; www.alesmav.com;
H. Wang, University of California

These circles differ only in grayscale intensity

Intensities give a strong “feeling” of scene
structure
Talk Outline

Introduction and Basics
Notations
 Reflectance map
 Photometric stereo



Main Approaches
Propagation Solutions
Basic Solution by Horn [85]
 Using Fast Marching, Kimmel & Sethian [2001]


Global energy minimization solution using Belief
Propagation – only if time permits
Introduction and Basics
Problem definition
Notations
Reflectance Map
Photometric Stereo
What determines scene radiance?



The amount of light that falls on the surface
The fraction of light that is reflected (albedo)
Geometry of light reflection
 Shape
of surface
 Viewpoint
n

=  ∙
( ,  )
( ,  )

Albedo

A property of the surface



Can be modeled as a single scalar for the entire surface
or dependent on location
Usually assumed to be a known scalar


Usually normalized 1 ≥  ≥ 0
Is this a ball with uniform albedo or is this a 2D circle with
varying albedo?
We’re not going to talk about the albedo
Surface Normal
Convenient notation for surface orientation
 A smooth surface has a tangent plane at every
point

We can model the surface using the normal at
every point
The Shape From Shading Problem

Given a grayscale image
 And
albedo
 And light source direction

Reconstruct scene geometry
 Can
be modeled by surface normals
Lambertian Surface




Appears equally bright from all viewing directions
Reflects all light without absorbing
Matte surface, no “shiny” spots
Brightness of the surface as seen from camera is
linearly correlated to the amount of light falling on
the surface
Today we will discuss only
Lambertian surfaces under
point-source illumination
n
( ,  )
( ,  )
Some Notations:
Surface Orientation

A smooth surface has a tangent plane at every
point
 Mark

=

; =

Parametrize surface orientation by first partial
derivatives of

Some Notations:
Surface Orientation

Surface normal
  = , 0,1 ,  = 0, , 1



=  ×  = , , −1
Normalize  =

=

(,,−1)

=
; =

2 +2 +1

Reflectance Map
Relationship between surface orientation and
brightness
Lambertian surface
1 +   +
,  =  ∙  =
1 + 2 + 2 1 +  2 +  2


Image irradiance (brightness) is proportional to
,  = (, )
Reflectance Map Example

Brightness as a function of surface orientation
Lambertian surface
iso-brightness
q
contour
0 .8
0 .9
R p, q   0.7
1 .0
 pS , qS 
p
i  90
 ppS  qqS 1  0
,  =
1 +   +
0 .3
1 + 2 +  2 1 +  2 +  2
Illustration from S. Narasimhan, Carnegie Mellon
0 .0
Reflectance Map of a Glossy Surface

Brightness as a function of surface orientation
Surface with diffuse and
glossy components
Reflectance Map Examples

Brightness as a function of surface orientation
Graphics with a 3D Feel



Given a 3D surface  ,  , lighting and viewing
direction, we can compute the gray level of pixel
(, ) of the surface
Find the gradient of the surface (, )
Use I x, y =  ,  =  ∙  =
determine the gray level
1+ +
1+2 +2
1+ 2 + 2
to

Can we reconstruct shape from a single image?
q
p

Two variables ,  with one equation… So what do
we do?
Illustration from S. Narasimhan, Carnegie Mellon

Use more images
 Photometric

stereo
 Introduce
constraints
 Solve locally
 Linearize problem
Images from
R. Basri et al. “Photometric Stereo with General, Unknown Lighting“, IJCV 2006
Photometric Stereo

Take several pictures of same object under same
viewpoint with different lighting
q
p , q 
1
1
S
S
p
Photometric Stereo

Take several pictures of same object under same
viewpoint with different lighting
q
p , q 
p , q 
2
2
S
S
1
1
S
S
p
Photometric Stereo

Take several pictures of same object under same
viewpoint with different lighting
q
p , q 
p , q 
2
2
S
S
1
1
S
S
p
p , q 
3
3
S
S
Main Approaches for
A Quick Review
Main Approaches

Minimization: Solve for a global error function while
introducing constraint(s):
Brightness ( − )2
 Smoothness
( 2 + 2 )
(( −  )2 +( −  )2 )
…

Main Approaches

Minimization: Solve for a global error function while
introducing constraint(s):
Brightness ( − )2
 Smoothness
( 2 + 2 )
(( −  )2 +( −  )2 )
…




Propagation: Grow a solution around an initial known
point or boundary
Local: Assume local surface type
Linear: Make problem linear by a linearization of
reflectance function
Main Approaches

Minimization: Solve for a global error function while
introducing constraint(s):
Brightness ( − )2
 Smoothness
( 2 + 2 )
(( −  )2 +( −  )2 )
…




Propagation: Grow a solution around an initial known
point or boundary
Local: Assume local surface type
Linear: Make problem linear by a linearization of
reflectance function
Main Approaches

Minimization: Solve for a global error function while
introducing constraint(s):
Brightness ( − )2
 Smoothness
( 2 + 2 )
(( −  )2 +( −  )2 )
…




Propagation: Grow a solution around an initial known
point or boundary
Local: Assume local surface type
Linear: Make problem linear by a linearization of
reflectance function
Basic Propagation Solution
Horn [85]
Solution by Characteristic Curves
Propagating Solution


Suppose that we know the depth of some point on
the surface
Extend the solution by taking a small step  =
+
and  are unknown, and the image irradiance
equation gives only one constraint
 But

If we knew  and , can compute changes in ,
using second partial derivatives


=
2
2
+
2

=

and  =
=
2
2
+ 2

2
2
2

2

2
2
Propagating Solution

=
, what else can we use?

 Image irradiance equation  ,  =  ,

=
= =
+

=

= =

=

+

Propagating Solution

=
;
=

 We note that for a small step

=

 For these specific δ,

=
=
=


If you take a small step in the image plane parallel to
the direction of the gradient of , then can compute
change in (, )
Propagating Solution

=


;
Get ODEs:

=  ;
=

;
=  +  ;



=

=  ;
=

These equations describe a characteristic curve
(, , , , )
(  + ,  + ,  + ,  + ,  +  )
Propagating Characteristic Curve

Characteristic curve
(, , , , )
( + ,  + ,  + ,  + ,  + )

Need to initialize every curve at some known point
 Singular
points
 Occluding boundaries
Image by Dejan Todorović
Propagating Characteristic Curve

Characteristic curve
(, , , , )
( + ,  + ,  + ,  + ,  + )

Need to initialize every curve at some known point
 Singular
points
 Occluding boundaries

Curves are “grown” independently, very instable
Image by Dejan Todorović
Shape From Shading by Fast
Marching
“Optimal Algorithm for Shape from Shading and
Path Planning”, R. Kimmel and J. A. Sethian
[2001]
Vertical Light Source Case

Recall reflectance map
1 +   +
,  =  ,  =  ∙  =



1 + 2 + 2 1 +  2 +  2
Assume light source located near the viewer

=(0,0,1)

,  =  ,  =
=
1
2
1
1+2 +2
=
1
1+2
−1
This is an Eikonal equation

Can be solved by an   numerical algorithm
Fast Marching
Fast Marching

Expanding Dijkstra’s shortest path algorithm for
general surfaces
 Represented
 Flat

as triangulated mesh
domains
Many computer vision problems can be set into
flat domains
 Every
pixel is a node, edges
Dijkstra’s Shortest Paths
Connected graphs
9
6
2
14
11
10
9
s
7
15
Dijkstra’s Shortest Paths
Connected graphs
 Update all neighbors
9
14
14
6
2
9
9
s
11
10
15
7
7
Dijkstra’s Shortest Paths
Connected graphs
 Update all neighbors
 Go to the closest neighbor
 Set
9
14
14
6
2
9
9
its value
s
11
10
15
7
7
Dijkstra’s Shortest Paths
Connected graphs
 Update all neighbors
 Go to the closest neighbor
its value
 Compute its neighbors
 Update smaller scores
9
14
14
6
2
9
9
 Set
s
23
11
10
15
7
7
Dijkstra’s Shortest Paths
Connected graphs
 Update all neighbors
 Go to the closest neighbor
9
14
14
its value
 Compute its neighbors
 Update smaller scores
6
2
9
9
 Set

Continue with smallest value node
 Remember
path
s
23
11
10
15
7
7
Flat Domains: Why Does Dijkstra Fail?


Dijkstra will not find the
diagonal path
Need to examine triangles
Fast Marching: Problem Definition




Suppose there is a forest fire with multiple sources
Every point that was touched by the fire is burnt
and will not be visited by the fire again
Firemen register the time T(x) at which the fire
arrives to location x.
The fast marching algorithm simulates this
scenario
Fast Marching: Problem Definition




Multiple sources
Advances either at a constant rate (,)
= 1 or at

varying rate (,)
= (, )

What is the arrival time at every location?
Image from G. Rosman, Technion
Fast Marching Algorithm






Set  0 = 0 and mark the
points as done
Set the rest of the points as
∞ and mark them as far
All far points adjacent to
done points become verify
points
Update all verify points using
of the done set
The verify with the smallest
becomes done
Continue until all points are
done
Illustrations from selected publications of J. A. Sethian
Fast Marching on a Grid
Update step
Tij
?
Ti,j-1
Ti-1,j
Slide based on slides by R. Kimmel, Technion
Ti+1,j
i+1,j
i,j-1
ij
i-1,j
i,j+1
Ti,j+1
Fast Marching on a Grid
Update step
− min −1, , +1, ,
Tij
2
+  − min ,−1 , ,+1 ,
2
=
ℎ2  2
Ti,j-1
Ti-1,j
Slide based on slides by R. Kimmel, Technion
?
Ti+1,j
i+1,j
i,j-1
ij
i-1,j
i,j+1
Ti,j+1
Fast Marching on a Grid
Update step
1
− min −1, , +1, ,
Tij
2
+  − min ,−1 , ,+1 ,
2
Slide based on slides by R. Kimmel, Technion
2
=
ℎ2  2
?
Ti,j-1
Ti-1,j
Ti+1,j
i+1,j
i,j-1
ij
i-1,j
i,j+1
Ti,j+1
Fast Marching on a Grid

Initialization: all  = ∞ or known initial value

Update step:
− 1
 Fit
2
+  − 2
2
= ℎ2  2
Tij
a triangle with gradient
and two determined values at
neighboring grid points
Ti,j-1
Ti-1,j
Slide based on slides by R. Kimmel, Technion
Ti+1,j
i+1,j
i,j-1
ij
i-1,j
i,j+1
Ti,j+1
Update Step Details

Update step is a quadratic equation:
− 1

Solution:  =
1
2
2
+  − 2
2
= ℎ2  2
2
1 + 2 + 2ℎ2  − 1 − 2
Disregard the “minus” solution
which yields  < 1 , 2
2
Update Step Details

Update step is a quadratic equation:
− 1
1
2
2
+  − 2
Solution:  =

If 1 − 2 > ℎ , then  <
1 > 2 and 1 > ℎ + 2

< 1
= ℎ2  2
2
1 + 2 + 2ℎ2  − 1 − 2


2
1
2
1 + 2 + ℎ
2
Update Step Details

Update step is a quadratic equation:
− 1
1
2
2
+  − 2
Solution:  =

If 1 − 2 > ℎ , then  <
1 > 2 and 1 > ℎ + 2

< 1
= ℎ2  2
2
1 + 2 + 2ℎ2  − 1 − 2


2
1
2
1 + 2 + ℎ
2
Update Step Details

Update step is a quadratic equation:
− 1
1
2
2
+  − 2
Solution:  =

If 1 − 2 > ℎ , then  <
1 > 2 and 1 > ℎ + 2


2
1
2
2
1 + 2 + ℎ
< 1
2 > 1 and 2 > ℎ + 1


= ℎ2  2
1 + 2 + 2ℎ2  − 1 − 2


2
< 2
This means that the wave front propagation comes from outside
the triangle

Only choice is to update  = min 1 , 2 + ℎ
Update Step Details

Update step is a quadratic equation:
− 1

+  − 2
2
= ℎ2  2
If 1 − 2 < ℎ


2
=
1
2
Else

2
1 + 2 + 2ℎ2  − 1 − 2
= min 1 , 2 + ℎ
2
Vertical Light Source Case

Eikonal equation
=

− 1 = (, )
Numerical approximation of this equation
max
−−1,  −+1,
, ∆ , 0
∆
+ max

1
(,)2
2
−,−1  −,+1
, ∆ , 0
∆
Assume w.l.o.g. that ∆ = ∆ = 1
2
=  2
Vertical Light Source Case

Eikonal equation
=

1
2
− 1 = (, )
Numerical approximation of this equation
max  − −1, ,  − +1, , 0
+ max  − ,−1 ,  − ,+1 , 0

Assume w.l.o.g. that ∆ = ∆ = 1
2
2
=  2
Vertical Light Source Case
Derive Update Step
2
max  − −1, ,  − +1, , 0
+ max  − ,−1 ,  − ,+1 , 0
1
2
2
− min −1, , +1, ,
+  − min ,−1 , ,+1 ,
2
=
2
− 1
2
+  − 2
=  2
2
=
2
2
Vertical Light Source Case
Derive Update Step
− 1
2
+  − 2
2
=  2
Where:
1 = min −1, , +1,
2 = min ,−1 , ,+1

If  > 1 − 2
=
1 + 2 + 2 2 − 1 − 2
Otherwise  = min 1 , 2
2
+
2
Vertical Light Source Case
Fast Marching



Fast marching introduces order to these update
steps
Points are updated from small to large
Total complexity: ( log )
 log
for the selection of smallest point and update of
neighboring points
  points (pixles)
General Light Source

Recall reflectance map
,  =  ∙  =


1 +   +
1 + 2 + 2 1 +  2 +  2
In this case  =  , , ,  =  , ,  ,
 The right hand side depends on
This is not an Eikonal equation anymore
 How can we solve this?
Solve in the light source coordinates, then
change the variables back
General Light Source
Change Coordinate System

Before moving to light source coordinate system,
choose w.l.o.g light source direction to be
2
2
= 1 , 0, 3 , where 1 +3 = 1
 This will simplify our equations a bit
General Light Source
Change Coordinate System


3

= 0

−1
0 1
1 0
0 3

=
− 1 = (, )
Eikonal equation!

1

,  =  ,  = (3  + 1 , )
1
(,)2
=
3

Light source coordinates marked with ~


Change coordinate system

=

−
General Light Source
Solution in New Coordinate System

But…
,  =  ,  = (3  + 1 , )
 We need to find  in order to find the intensity
Use the smallest  value from all neighbors of
General Light Source Solution
Update Step
− 1
2
+  − 2
2
=  2
Where:
1 = min −1, , +1, ,
2 = min ,−1 , ,+1 ,

= 3  + 1 min 1 , 2

If  > 1 − 2
=
1 + 2 + 2 2 − 1 − 2
2
Otherwise  = min 1 , 2 +
2
Results


Results shown on
image
= (0.2,0,0.96)
Minimization Using Efficient
Belief Propagation
“Efﬁcient Belief Propagation for Vision Using
Linear Constraint Nodes”, B. Potetz [2007]
Belief Propagation

Estimate the marginals of a multivariate probability
distribution

=
( ),
⊂
  ( ) potential function
  consists of several elements

Usually represented as a factor graph
Factor Graph




=  ( ),
⊂
Bipartite graph
Every element of  is represented by a variable
node
Every potential function   represented by a
factor node , connected to all variable nodes of
circles: variable nodes
squares: factor nodes
Belief Propagation on a Factor Graph


=  ( ),
⊂
Estimate the marginals


=
∖ ()
Iteratively compute messages along edges of the
factor graph
 Every
variable node computes its marginal probability
according to the probabilities its neighbors sent

Usually converges to a good enough local
minimum
Shape From Shading Using a
Factor Graph

Recall I x, y =  ,  =  ∙
1+ +
=
2
2
2
2
1+ +

1+ +
Constraint I: Consistent integral
−  −  −  +  +  −  +  = 0
Where  −  =  ,  −   − 1,
Integral of the surface gradient along a closed curve
should be zero
Illustration from Klaus & Horn, “Robot Vision”, MIT Press, 1986
Shape From Shading Using a
Factor Graph

Constraint II: Lambertian constraint
,  =  ,
Where (, |) =
,  = (, )
0
ℎ

,  could be changed for other lighting
conditions

Constraint III: Shape prior
 Modeled
by a Laplace distribution
(− ∆ ), (− ∆ )
1
2
Shape From Shading Using a
Factor Graph

Constraint I: Consistent integral
−  −  −  +  +  −  +  = 0

Constraint II: Lambertian constraint
,  =  ,
Constraint III: Shape prior

(− ∆ ), (− ∆ )
1
2
Variable nodes
Factor node, consistent integral
Factor node, Lambertian constraint
Factor node, shape prior
one pixel one pixel
one pixel one pixel one pixel
Shape From Shading Using a
Factor Graph

Constraint I: Consistent integral
−  −  −  +  +  −  +  = 0

Constraint II: Lambertian constraint
,  =  ,
Constraint III: Shape prior

(− ∆ ), (− ∆ )
1
2
Variable nodes
Factor node, consistent integral
Factor node, Lambertian constraint
Factor node, shape prior
Results

Input:
original surface

Rendering of original
surface with  = (1,0,1)
Input image
Result:
recovered surface
Rendering of recovered
surface with  = (1,0,1)
Conclusion

Shape from shading problem
 Definition
 Main
difficulties
 Main approaches

Propagating solutions
 Horn
[85]
 Kimmel [2001]

Minimization with BP while introducing constraints
Thank You!
```