A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE
FROM SHADING PROBLEM
Markus Louw, Fred Nicolls
Dept. Electrical Engineering, University of Cape Town, South Africa
Dee Bradshaw
Dept. Chemical Engineering, University of Cape Town, South Africa
Keywords:
Loopy Belief Propagation, Bayesian network, Shape from shading, photoclinometry.
Abstract:
This paper describes a new approach to the shape from shading problem, using loopy belief propagation
which is simple and intuitive. The algorithm is called Loopy Belief Propagation Shape-From-Shading (LBP-
SFS). It produces reasonable results on real and synthetic data, and surface information from sources other
than the image (eg range or stereo data) can be readily incorporated as prior information about the surface
elevation at any point, using this framework. In addition, this algorithm proves the use of linear interpolation
at the message passing level within a loopy Bayesian network, which to the authors’ knowledge has not been
previously explored.
1 INTRODUCTION
The interested reader is referred to two surveys,
(R. Zhang and Shah, 1999), and (Jean-Denis Durou,
2004). In (R. Zhang and Shah, 1999), SFS approaches
are classified into minimization e.g. (Szeliski, 1994),
propagation e.g. (S. Osher, 1988), local e.g. (Pent-
land, 1984), or linear e.g. (P.S. Tsai, 1994) ap-
proaches. In (Jean-Denis Durou, 2004), SFS meth-
ods are classified into methods based on partial dif-
ferential equations (characteristic strips (B.K.P.Horn,
1975), power series expansion (Bruss, 1982), and
viscosity solutions e.g. (M. G. Crandall, 1983)),
minimization methods (P. Daniel, 2000), and ”meth-
ods approximating the image irradiance equation”,
which contain the local and linear methods surveyed
in (R. Zhang and Shah, 1999).
These surveys describe the development of shape
from shading methods, in which researchers have
tried to mimic the way the human brain and eyes ex-
tract shape information from shading on the object,
as well as trying to find analytical solutions based on
geometry and reflectance characteristics. This paper
describes the casting of the SFS problem into the be-
lief propagation paradigm, which would place it in
the minimization and also the propagation class of
method. Our method is algorithmically similar to
(Jian Sun, 2002), in which Sun et al. use Loopy Be-
lief Propagation (LBP) to solve the dense stereo cor-
respondence problem. In (Jian Sun, 2002), each pixel
in the left image is probabilistically assigned dispar-
ities for matching to a pixel in the right image, and
belief propagation is performed on the nodes which
are connected to their immediate (Ising) neighbours.
Our method uses a more complicated energy function
to approximate the correct elevation map for the sur-
face given the irradiance map, and to enforce surface
smoothness conditions. It also incorporates a multi-
resolution interpolation based approach which to our
knowledge has not used before in Loopy Belief Prop-
agation.
2 LAMBERTIAN LIGHTING
MODEL
This algorithm calculates a surface on the Lambertian
assumption that the intensity of a pixel is proportional
to the inner product of the direction vector of the in-
cident light and the surface normal at that point. We
may follow the notation of (Jean-Denis Durou, 2004),
to formulate this. The image irradiance equation is
R(
−→
n (x)) = I(x) (1)
376
Louw M., Nicolls F. and Bradshaw D. (2007).
A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE FROM SHADING PROBLEM.
In Proceedings of the Second International Conference on Computer Vision Theory and Applications - IFP/IA, pages 376-381
Copyright
c
 SciTePress

where I(x) is the image irradiance (usually the inten-
sity) measured at location x, and R(
−→
n (x)) is the re-
flectance function on the surface which takes the nor-
mal at point x as an argument. The surface normal
may be calculated as
1
√
(1 + p(x)
2
+ q(x)
2
)
(−p(x), −q(x),1) (2)
where
p = ∂u/∂x
1
(3)
and
q = ∂u/∂x
2
(4)
where u is the height of the surface. If there is a
unique light source at infinity, and shining in direc-
tion
−→
w = (w1,w2,w3), the pixel intensity is the inner
product
R(
−→
n (x)) =
−→
w ·
−→
n (x) (5)
Hereafter, without loss of generality (but assuming
all surface points are visible to both camera and light
source) we assume the light source is in the same di-
rection as the camera, which produces an orthogonal
projection.
3 LOOPY BELIEF
PROPAGATION
The Loopy Belief propagation algorithm using factor
nodes may be expressed in the following equations
(following loosely the notation of (Murphy, 2002)):
µ
x→f
(x) = o(x)
∏
g6= f
µ
g→x
(x) (6)
µ
f →x
(x) =
∑
u\x
f (u)
∏
y6=x
µ
y→f
(y) (7)
where x and y are the variable nodes, f and g are fac-
tor nodes, o(x) is the prior probability (observation)
on the variable node x, and where it is assumed µ is in
the domain of f . Our algorithm uses a parallel updat-
ing scheme, with the max product algorithm. After
the specified number of update iterations, the poste-
rior distribution on each corner vertex node may be
given as:
p
x
(x) = o(x)
∏
g
µ
g→x
(x), (8)
where g is the set of factor node neighbours of x.
4 FORMULATION OF LBP TO
SOLVE SFS
This algorithm calculates a posterior for the height
at each corner vertex on the image. A corner ver-
tex occurs at the corner of a pixel; at the intersec-
tion of four pixels, one corner vertex represents the
Figure 1: Depiction of the connectivity between vertex
nodes (round) and factor nodes (square). These four ver-
tex nodes correspond to the heights of the four corners of a
single pixel.
height of the surface at that location. Our energy func-
tional is evolved using Factor nodes which represent
the probability of triplets of these corner vertex nodes.
Each triplet of vertices creates a unique plane, and the
orientation of that plane relative to the direction of
the light source allows a probability to be assigned
to that configuration for that triplet. If the dimen-
sions of the image are width w and height h and if
only simple right angled triangles are used, with the
topology shown in Fig. 1, the number of vertices is
(h + 1)(w+ 1) and the number of factor nodes is 4wh.
If smoothing using point triplets (described later) is
incorporated, the number of Factor nodes is 8wh. A
diagram of the topology for this scheme is shown in
Fig.1. The plane generated by each triplet of corner
nodes forms an angle against the incident light, giving
an illumination for that pixel. This is shown in Fig. 3.
Next we define a Markov Random Field (MRF) on
this set of vertex nodes X, given the image data Y and
explicit range data (which gives a prior probability for
the height of the surface at a particular location on the
surface):
P(X|Y,Z) ∝
∏
i, j,k:k> j>i
exp−ψ
t
i jk
(x
i
,x
j
,x
k
,y
i jk
)···
∏
i
exp−ψ
i
(x
i
,z
i
) (9)
Each element of the state vector of a vertex node
corresponds to the vertex node taking on a particular
height. At each iteration, our implementation of Eqn.
7 is the maximum product (aka max. prod. algorithm)
of the input messages with the elements of factor node
u. A factor node is therefore a 3D array which con-
tains probabilities, each element is derived from the
energy
ψ
t
i jk
(x
i
,x
j
,x
k
,y
i
) ∝ |y
i jk
−|
−→
n
i jk
(x
i
,x
j
,x
k
)·
−→
L ||, (10)
A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE FROM SHADING PROBLEM
377

where N is the number of images (one for each light-
source), n
i jk
(x
i
,x
j
,x
k
) is the normal of the surface in
the triangle between i, j, k, y
i jk
is the average image
intensity in the image region enclosing vertex nodes
i, j,k.
−→
L is the lightsource direction, x
i
is the height
of the vertex at right angles to the other two, x
j
is the
height of the vertex horizontal to that vertex, and x
k
is
the height of the remaining vertex. We thus have the
following equations for the partial derivatives in the
height (with respect to change in position in the hor-
izontal and vertical directions on the image) in terms
of the three heights of the surface at the points on
which the three corner vertex nodes lie:
p = ∂u/∂x
1
= (x
j
−x
i
)/∂x
1
(11)
q = ∂u/∂x
2
= (x
j
−x
i
)/∂x
2
. (12)
Assuming an affine camera with square pixels and
an overall scale of one unit per pixel width we set ∂x
1
and ∂x
2
to 1. We can then use Eqn. 2 to calculate the
plane.
Note that it is simple at this stage to extend the
energy function to include static scene/moving light-
source information (on the assumption that all points
on the surface are always visible to both the camera
and to all lightsources). We therefore adjust Eqn. 10
above to be:
ψ
t
i jk
(x
i
,x
j
,x
k
,
−→
y
i jk
) ∝
N
∑
s=1
|y
i jks
−|
−→
n
i jk
(x
i
,x
j
,x
k
)·
−→
L
s
||,
(13)
Where N is the number of images (one for each light-
source), n
i jk
is the normal of the surface in the trian-
gle between i, j, k, y
i jks
is the average image intensity
in the image region enclosing vertex nodes i, j,k, in
image s (note that we have vectorized
−→
y
i jk
to show
that it contains average intensity values over all inten-
sity images, indexed with the subscript s. So s iterates
over each of the images, and
−→
L
s
is the lightsource di-
rection in image s.
4.1 Boundary Conditions and Range
Data
In Shape from Shading algorithms, it is usually neces-
sary to establish some boundary conditions, since all
surface heights calculated (if only shape from shad-
ing be used) are relative to each other, but not against
any fixed frame of reference. In addition, the specifi-
cation of boundary conditions may solve some of the
ambiguities, since there are usually a number of sur-
faces which may generate a particular intensity map
under particular lighting conditions. LBP allows such
boundary conditions and range data to take on the
form of either hard or soft constraints, and this is han-
dled naturally within the LBP paradigm. Each corner
vertex node x
i
may be given a prior probability on the
heights of its state vector, such that
x
i
(X = h) ∝ exp(−|h −u(x
i
)|), (14)
where u(x
i
) is the specified range or depth of the
point. Whether the point is given a value because it
lies on a known boundary, or because we have range
data about the point, the point is treated the same way.
4.2 Multi-Resolution
The LBP method allows us to approach the problem
in a multi-resolution manner: after N LBP iterations,
we may iterate through each of the M corner vertex
nodes and adjust the heights which each element in
the vertex node’s state vector corresponds to, so that
we may home in on a closer approximation of the
correct value. The multi-resolution algorithm is de-
scribed as follows:
1. for i = 1 to N
2. for j = 1 to M
3. calculate MAP estimate for all corner vertex node
4. end
5. increase height resolution of each corner vertex node
6. end
In step 5 in the above algorithm, for each corner
vertex node, we calculate the posterior on that node
using Eqn. 8. Then, we find the entry with the highest
probability. The new set of heights for that node are
then calculated using a finer height resolution about
the MAP value. The algorithm for calculating the new
set of heights is
1. Calculate new height resolution: h
new
= kh
old
2. n ← 1; a ← 1; b ← 1
3. while (n < numlabels)
4. if (LB < H
x
+ a ·h
new
< UB)
5. vertex height(n) ← H
x
+ a ·h
new
6. n ← n + 1; a ← a + 1
7. end
8. if (n < numlabels)
9. if (LB < H
x
−b ·h
new
< UB)
10. vertex height(n) ← H
x
−b ·h
new
11. n ← n + 1; b ← b + 1
12. end
13. end
14. end
where in the above algorithm, H
x
is the MAP esti-
mate for the height at x, h
new
is the height resolution
for the current LBP resolution, k is a compression ra-
tio applied to the previous height resolution h
old
at
each resolution cycle, LB and U B are the lower and
upper bounds respectively. We then update the mes-
sages from the corner vertex nodes to the factor nodes
VISAPP 2007 - International Conference on Computer Vision Theory and Applications
378

the reconstruction (e.g. Fig. 7). Even after the bound-
ary conditions and some range data have been given,
there are many different surfaces which could produce
the observed image intensity data. Ways of overcom-
ing this within the LBP-SFS paradigm include adding
a smoothing energy term on collinear point triplets
and the use of surface triangles of varying scales to
eliminate high frequency error in the resultant digital
elevation map. Other methods e.g. (P.L. Lions, 1993)
usually remove some ambiguity by defining a sin-
gle highest point or characteristic curve (M. G. Cran-
dall, 1983),(M. G. Crandall, 1984),(E. Rony, 1992),
(S. Osher, 1988). Our algorithm may be made more
reliable by increasing the number of labels used to
represent corner vertex heights, although each extra
label greatly increases the computation time. We have
also noted that elevation map results using this al-
gorithm are more reliable for smaller images (less
than 50×50). As the image size increases, the algo-
rithm is more prone to fall into bad local minima. We
are experimenting with a wavelet-style spatial multi-
resolution approach to overcome this. The algorithm
results will improve with an increase in the number
of soft range data points supplied, and with increased
certainty on those points. This method should re-
ally be seen as a data fusion method for incorporating
range data with intensity information.
7 CONCLUSION
The algorithm has the advantage of being adjustable
in terms of the height resolution required per vertex.
Unlike other SFS methods, LBP-SFS can incorporate
both hard and soft constraints on the boundary con-
ditions of the surface and range data at points on the
surface. Different reflectance models per surface can
be easily accounted for in the energy term. This algo-
rithm requires long computation time and large stor-
age for images with large depth variation (such im-
ages would require larger vertex node state vectors
given an initial height resolution).
8 CURRENT AND FUTURE
WORK
Preliminary results show that minimizing the same
MRF formulations using simulated annealing with
Gibbs sampling is faster and more reliable. We are
also investigating integrating SFS information into a
dense stereo formulation.
ACKNOWLEDGEMENTS
The authors are grateful for the financial support
given by the National Research Foundation of South
Africa, and Anglo American via the Minerals Pro-
cessing Research Unit at the University of Cape
Town.
REFERENCES
B.K.P.Horn (1975). Obtaining Shape from Shading Infor-
mation. McGraw-Hill.
Bruss, A. (1982). The eikonal equation: Some results ap-
plicable to computer vision. Journal of Mathematical
Physics, pages 890–896.
E. Rony, A. T. (1992). A viscosity solutions approach to
shape-from-shading. SIAM. J. Numer. Anal., pages
867–884.
Jean-Denis Durou, Maurizio Falcone, M. S. (2004). Nu-
merical methods for shape from shading: A survey
with benchmarks. CVIU.
Jian Sun, Heung-Yeung Shum, N.-N. Z. (2002). Stereo
matching using belief propagation. ECCV.
M. G. Crandall, P. L. (1983). Viscosity solution of hamilton-
jacobi equationssfs. Trans. Amer. Math. Soc., pages
1–42.
M. G. Crandall, P. L. (1984). Two approximations of so-
lutions of hamilto n-jacobi equations. Math. Comput.,
pages 907–922.
Murphy, K. (2002). Dynamic Bayesian Networks: Repre-
sentation, Inference and Learning. PhD thesis, Uni-
versity of California, Berkeley.
P. Daniel, J.-D. D. (2000). From Deterministic to Stochastic
Methods for Shape from Shading. In Proc. 4th Asian
Conf. on Comp. Vis.
Pentland, A. (1984). Local shading analysis. IEEE. trans-
actions on Pattern Analysis and Machine Intelligence,
pages 170–187.
P.L. Lions, E. Rouy, A. T. (1993). Shape-from-shading, vis-
cosity solutions and edges. Numerische Mathematik
64, pages 323–353.
P.S. Tsai, M. S. (1994). Shape from shading using linear
approximation. Image and Vision Computing Journal,
pages 187–198.
R. Petrovic, I. Cohen, B. F. R. K. and Huang, T. (2001). En-
forcing integrability for surface reconstruction algo-
rithms using belief propagation in graphical models.
CVPR.
R. Zhang, P. Tsai, J. C. and Shah, M. (1999). Shape from
shading: A survey. PAMI.
S. Osher, J. S. (1988). Fronts propagating with curvature-
dependent speed: algorithms based on hamiltonian-
jacobi formulations. J. Comput. Phys., pages 12–49.
Szeliski, R. (1994). Fast shape from shading. Computer
Vision, Graphics, Image Processing: Image Under-
standing, pages 129–153.
A LOOPY BELIEF PROPAGATION APPROACH TO THE SHAPE FROM SHADING PROBLEM
381