IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 20131

Photometric Stereo using Sparse Bayesian
Regression for General Diffuse Surfaces
Satoshi Ikehata, Student Member, IEEE, David Wipf, Member, IEEE, Yasuyuki Matsushita, Senior
Member, IEEE and Kiyoharu Aizawa, Member, IEEE
Abstract—Most conventional algorithms for non-Lambertian photometric stereo can be partitioned into two categories. The firstcategory is built upon stable outlier rejection techniques while assuming a dense Lambertian structure for the inliers, and thus
performance degrades when general diffuse regions are present. The second utilizes complex reflectance representations and nonlinear optimization over pixels to handle non-Lambertian surfaces, but does not explicitly account for shadows or other forms ofcorrupting outliers. In this paper, we present a purely pixel-wise photometric stereo method that stably and efficiently handles various
non-Lambertian effects by assuming that appearances can be decomposed into a sparse, non-diffuse component (e.g., shadows,
specularities, etc.) and a diffuse component represented by a monotonic function of the surface normal and lighting dot-product. Thisfunction is constructed using a piecewise linear approximation to the inverse diffuse model, leading to closed-form estimates of thesurface normals and model parameters in the absence of non-diffuse corruptions. The latter are modeled as latent variables embedded
within a hierarchical Bayesian model such that we may accurately compute the unknown surface normals while simultaneously
separating diffuse from non-diffuse components. Extensive evaluations are performed that show state-of-the-art performance using
both synthetic and real-world images.
Index Terms—Photometric Stereo, Sparse Regression, Piecewise Linear Regression, Sparse Bayesian Learning.

F

1

I NTRODUCTIONPHOTOMETRIC stereo involves estimating the surface normals of an object given appearance variations in multiple images taken under different lighting conditions. While
more than 20 years have passed since the original workon Lambertian photometric stereo [1], recovering the surface
orientations of a natural scene from appearance variations
remains challenging, largely because real-world objects exhibit diverse non-Lambertian effects. The design of practicalphotometric stereo algorithms involves complicated trade-offs
related to robustness, efficiency, and versatility, with mostapproaches falling into one of two categories. First, there
are many techniques that assume a basic Lambertian model
but augmented with outlier detection for handling all nonLambertian regions of the scene [2], [3], [4], [5], [6]. While
this strategy is numerically stable and relatively insensitive
to initializations, it may be computationally expensive since alarge number of images is required for robust outlier rejection.
Additionally, complex reflections such as rough specularities
or non-Lambertian diffuse reflections can be highly disruptive.
In contrast, a second basic approach relies on sophisticatedreflectance models or the sharing of properties across materials
to explicitly account for non-Lambertian reflections [7], [8],[9], [10], [11], [12]. While these methods are more capable ofhandling a wide variety of objects including rough surfaces,
they are often more sensitive to initialization, numerical instabilities derived from the complex nonlinear optimization, and• S. Ikehata and K. Aizawa are with University of Tokyo.
• D. Wipf and Y. Matsushita are with Microsoft Research Asia.
Manuscript received **********, 2012; revised **********, 2013.

the effects of shadows and other corrupting outliers. Further
details regarding existing approaches will be discussed in
Section 2.
This paper presents an alternative photometric stereo algorithm for stably and accurately estimating the surface normals
of a scene in the presence of various non-Lambertian effects.
For this purpose, a hierarchical Bayesian model is developedin Section 3 that automatically decomposes observed appearances into a continuous piecewise linear diffuse component
and a sparse, non-diffuse component for capturing shadows,
specularities, and other corruptions. Optimization and inference is accomplished using a robust majorization-minimization
technique akin to the popular EM algorithm, with desirable
convergence properties and quantifiable advantages over standard convex estimators.
For mathematical convenience, we develop our diffuse component by assuming, in the absence of non-diffuse corruptions,that pixelwise appearances are well-approximated by a monotonic (and therefore invertible) function of the dot-productbetween the surface normal and the lighting direction. We maythen consider the inverse representation of the image formation
process, where the unknown normal vector is now separated
from the unknown monotonic inverse reflectance function. By
parameterizing the latter using a piecewise linear approximation, we obtain a set of linear equations in both the surface
normals and diffuse parameters, leading to simple, closedform estimators. We note that the piecewise linear function
has been widely used for regression purposes [13] because of
its simplicity and representational capability. In the context ofphotometric stereo, this function is a the natural extension of
the Lambertian reflectance model, which directly corresponds
to the case where there is only one linear segment. Later weIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

will empirically demonstrate that additional piecewise linear
diffuse segments (e.g., three) can effectively represent many
complex non-Lambertian reflections.
Of course we do not know a priori where non-diffuse
appearances will be located, and thus the overall estimation
problem is underconstrained even with many available images.
Thus, our hierarchical Bayesian model attempts to maximize
the number of observations that can be adequately explained
via the piecewise linear inverse diffuse function (the inlier
model) while treating the remaining observations, including
specularities, shadows, and sensor saturations, as non-diffuse
elements with separate, unknown variances learned from thedata (the outlier model). The partitioning between inliers andoutliers is estimated simultaneously with the normal vectors
and piecewise linear parameters using a principled variational
Bayesian technique.
The proposed framework benefits from simple, efficientpixelwise optimization, which is easily amenable to parallel
processing. Moreover it does not require the pre-processingof specularities/shadows, careful initialization strategies, or
typical smoothness constraints for both object structure andreflectance, which can disrupt the recovery of fine details.
Consequently, we do not suffer from numerical instability evenwith relatively few input images (e.g., 20) and we do nothave difficulty handling spatially varying materials with highfrequency structures. Extensive experiments in Section 4 showthat our implementation produces more stable, accurate, andefficient surface normal estimates than other robust algorithms
such as [3], [6], even when applied to non-Lambertian materials that previously required complex model-based approaches
such as [11]. A preliminary version of this work, without
the piecewise linear diffuse model, many algorithmic details,
or comprehensive empirical testing, has appeared in recent
conference proceedings [14].

2

R ELATED W ORK

Since Woodham [1] first introduced photometric stereo for
Lambertian scenes, the extension to non-Lambertian surfaces
has drawn significant interest. These extensions are largely
categorized by whether the non-Lambertian reflections are
considered as outliers for removal, or as inliers using nonlinear reflectance modeling.
In the former class, surface normals are recovered via a
simple Lambertian reflectance model while non-Lambertian
corruptions are robustly neglected. The most traditional yet
still widely used technique is intensity thresholding for avoiding the effects of shadows and specularities because of its
simplicity. To handle spatially varying albedo or sensor saturations, some recent approaches incorporate classification
methods for detecting shadows and specularties [15], [16],[17], or exploit graph-based visibility estimation for shadow
removal [18]. A color-based technique has also been used for
specular removal [2]. More recently, techniques for handlingother corruptions, such as inter-reflections [19] and calibration
errors [20], have been proposed.
For robustly identifying outliers, the RANdom SAmple
Consensus (RANSAC) scheme has been used in [3], [21], [5].

2While RANSAC is effective when the number of observations
and unknown inlier parameters are relatively small, it becomes computationally intractable with higher dimensionality.
Alternatively, some works cast photometric stereo as a matrix decomposition or factorization problem considering nonLambertian corruptions as missing entries. For example, Wu etal. [6] propose a rank-minimization approach that decomposes
observations into a low-rank Lambertian structure and sparse
non-Lambertian corruptions/outliers, while del Bue et al. [22]solve photometric stereo as a manifold constrained bilinear
factorization problem in the presence of missing entries.
Spatial constraints are also exploited to effectively remove
outliers. Tang et al. [23] use an Markov random field (MRF)
for imposing discontinuity preserving smoothness using belief
propagation. A similar setting has also been used in [24] wheregraph cuts is used for deriving the optimal solution. While the
MRF-based approaches preserve discontinuities, they tend toover-smooth the surface normal map. Spatial information is
also incorporated in filtering approaches [25], [4]. Miyazaki etal. [25] use a filter-based method where an inaccurate surface
normal is refined via the median estimate of neighboring
pixels. Yu et al. [4] propose a pixelwise scheme to find a
maximum feasible subset of Lambertian observations via a
Big-M algorithm.
In contrast to the former outlier removal techniques, asecond class of methods treat non-Lambertian reflections as inliers using non-linear reflectance modeling. Various physicallyinspired parametric reflectance models of the bidirectionalreflectance distribution function (BRDF) have been used inphotometric stereo. For example, Torrance-Sparrow [26] andWard [27] models have been used to account for specularities [7], [8]. To handle spatially varying materials, Goldman etal. [9] approximate basis materials by the Ward model, andreflectance parameters and surface normals are simultaneously
estimated in a non-linear manner. More recently, Alldrin etal. [10] use a non-parametric representation of BRDFs for
simultaneously recovering shape and and BRDFs, which are
represented by bi-variate form [28]. The model and optimization technique are further simplified in Shi et al. [11], wherea compact biquadratic representation of isotropic BRDF is
introduced.
Instead of explicitly modeling the reflectances, there aremethods that use general reflectance properties that are shared
among diverse real reflectances. The properties of isotropy
and reciprocity symmetry [29], [30] have been exploited
for surface normal reconstruction in [31], [32], [33]. The
monotonicity property of general diffuse reflectances is also
used in [33], [12]. Other than these properties, radiancesimilarity [34], image deviates [35] and diffuse maxima [36]are also used in the photometric stereo context. Examplebased approaches [37], [38] also fall in this category, in whichsurface normals are estimated using a reference object whose
shape is known and BRDF is the same as the target object.
Our method tries to capitalize on both of the above approaches in that it relaxes the often restrictive Lambertianreflectance model by using a piecewise linear representationof reflectances, while simultaneously performing robust estimation to avoid over-fitting. More specifically, our sparsity

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

penalty effectively rejects outliers while the inverse piecewise
linear representation of diffuse reflectances enables us tohandle a wide variety of materials.

3 P HOTOMETRIC S TEREO V IA I NVERSEP IECEWISE L INEAR S PARSE R EGRESSION
In this section, we formulate the estimation of surface normals
using photometric stereo as a piecewise-linear sparse Bayesian
regression problem. Henceforth, we rely on the following
assumptions:
1) Relative position between the camera and the object isfixed across all images.
2) Object is illuminated by a point light source at infinityfrom varying and known directions.
3) Camera view is orthographic, and the radiometric response function is linear.
3.1

Problem Statement

Diverse appearances of real world objects can be encoded by a
BRDF, which relates the observed intensity I at a given pointon the object to the associated surface normal n ∈ R3 , theincoming lighting direction l ∈ R3 , and the outgoing viewing
direction v ∈ R3 via
I = f (n, l, v).

(1)

The dichromatic reflectance model [39] states that if thescene is illuminated by a single dominant point light source,
the radiance is a linear combination of diffuse and specularreflectance via
I = fd (n, l) + fs (n, l, v).

(2)

In practice, various additional effects are observed, including
attached/cast shadows, image noise, inter-reflections, and soon. We can interpret these effects as additive corruptions e
applied to the ideal scene leading to the image formation
model giving
I = fd (n, l) + fs (n, l, v) + e.

(3)

The photometric stereo is a problem to recover surface normaln of a scene by inversely solving Eq. (3) from a collectionof m observations under the unknown set of parameters(fd , fs , e). Note that except for uncalibrated photometricstereo problems such as [36], l and v are usually known.
Early photometric stereo works assumed that the diffuse
component of observation is represented by Lambertian reflectance model (i.e., fd = ρnT l, where ρ is a surface
albedo) and discarded the non-Lambertian component fs + e
as outliers [3], [6], [14]. While the Lambertian assumption
is effective for a certain materials, this strong assumption on
the reflectance substantially limits the target objects. Instead,
we introduce the general representation of the material diffuse
function to handle non-Lambertian diffuse materials as follow,
fd (n, l) = f (nT l),(4)where we assume f (nT l) passes through the origin, that is
f (0) = 0. Note that Lambertian image formation is a special

3case when f is a linear function. Then, we also merge fs and e
as deviations from the diffuse reflection and use the following
imaging model,
I = I˜ + e = f (nT l) + e.

(5)Here, I˜ is the diffuse component of I. Our goal is to recoverunknown surface normal n, diffuse reflectance function fand non-diffuse corruptions e from a collection of lighting
directions l and associated appearances I. However, there are
two critical issues which must be solved: (a) the coincidenceof unknown parameters n and f in the same term, and (b)an under-constrained problem since the number of unknowns(equal to m + 2 plus however many degrees of freedom
are needed to describe f ) always exceeds the number of
equations (equal to the number of images m). We overcome
these difficulties by combining a convenient, piece-wise linear
inverse representation of the imaging model with a sparsity
penalty applied through latent variables embedded in a robust
hierarchical Bayesian framework.
3.2 Inverse Diffuse Reflectance Model
For simplicity, we first neglect the non-diffuse corruptions
e in Eq. (5). Then, we assume the monotonicity of the
diffuse reflectance function like [33], [12] which provides
the following constraint on the function under two different
lighting directions l1 and l2 ,nT li > nT lj ↔ f (nT li ) > f (nT lj ).

(6)Under this assumption, the unique existence of the inverse
function of f is guaranteed giving
˜ = g(I)˜ = nT l.
f −1 (I)(7)

Now that we assume only diffuse reflections appear in thescene, i.e., I˜ = I, the following equation is acquired,nT l = g(I).

(8)

We call Eq. (8) the inverse diffuse reflectance model. The fundamental advantage of this model is that an unknown function
g(I) and a surface normal n are separated, which contributes
to the simplicity of the problem. Eq. (8) suggests that the perpixel collection of 2-D plot (Ii , nT li ) corresponding to i-th
image must draw a monotonic inverse diffuse function g(I).
While this relationship limits the solution space of both n andg, the problem is that there are still multiple feasible solutionsof a pair of n and g(I) as illustrated in Fig. 1, especiallywhen m is small. To reduce inherent ambiguity of the problem,
we further assume a parametric model of the general inverse
diffuse function g(I).
Given that the left-hand-side of Eq. (8) is linear in theunknown normal vector n, for computational simplicity wewould like to impose similar linearity on the right-hand-side
in our parameterized representation of g(I). For this purpose,
we then choose to express g(I) as a summation over p fixedand known, non-linear basis functions gk (I) weighted by an
unknown coefficient vector a, leading to the representationg(I) =

pX
k=1

ak gk (I).

(9)

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

𝑔𝑎 (𝐼)𝒏𝑇 𝒍

𝑔𝑏 (𝐼)𝒏𝑇 𝒍(𝐼𝑖 , 𝒏𝑇𝑎 𝒍𝑖 )(𝐼𝑖 , 𝒏𝑇𝑏 𝒍𝑖 )
𝐼

𝐼

(a)

(b)

Fig. 1. We can draw different monotonic curves whichinterpolate 2-D plots (Ii , nT li ) derived from (na ) in (a)and (nb ) in (b), which illustrates that there are multiple
feasible solutions of a pair of n and g(I) which satisfy the
inverse diffuse reflectance model in Eq. (8).
𝑔1 (𝐼)

𝑔(𝐼)

𝑏1

𝑏2

𝑔2 (𝐼)+𝑏0
𝑏1
𝑔3 (𝐼)

𝑏0

𝐼

𝑏2

𝑏1

𝑏2

(a)

𝐼

𝐼

𝑏0

𝑏1

𝑏2

𝐼

(b)

Fig. 2. The illustration of the piecewise linear function.
We show the case when the number of basis functions
is three. (a) Each basis function is defined as a polylinear
function which has a breaking point, and (b) the piecewise
linear function is defined as the summation of these basis
functions.

4

when all elements of a are set to the same positive value. We
thus preserve preferable properties of the Lambertian model
where appropriate, while still allowing us to handle more general non-linear diffuse reflections where needed. In fact, evenwith p small, we may nonetheless approximate a wide variety
of non-linear functionals, with monotonicity ensured whenevera has all non-negative elements (while we do not strictly
enforce non-negativity, our learning procedure described in
the following sections strongly disfavors any ak < 0). Here
we should note that our model relates to the recent bilinear
BRDF model proposed by Romeiro and Zickler [40], whichis represented as a linear combination of non-negative basis
functions learned through non-negative matrix factorization of
100 materials in the MERL BRDF database [41]. However,there are two different points. Firstly, our basis functions arelearned directly from the data itself (i.e., linear segment ends
in Eq. (10) are decided from samples at each pixel) while[40] learns basis functions from the external database. More
importantly, we only model the diffusive component by asimple linear function which contributes to the computationally efficiency and stability, while [40] adapts basis functions
represented as discretized two-dimensional matrices.
By substituting Eq. (9) into the inverse diffuse reflectancemodel, Eq. (8) becomesnT l =pX

ak gk (I).

Collecting variations of observation at the same pixel under
different lighting directions, the equations can be merged into
following linear problem,Ãx = 0,where each bk−1 denotes the point where the (k − 1)-th linear
segment ends and the k-th segment begins. See Fig. 2 fordetails. By construction, g(I) will be a continuous piecewise linear function, meaning each adjoining linear segment
connects to one another at the corresponding point bk−1 ,regardless of the coefficients a. We also assume that b0 = 0,
meaning g(I) will intersect the origin. Remaining values ofbk−1 are chosen such that each piecewise linear segment spansan equal-sized bin over the range of I.
This choice of basis function leads to a natural generalization of the Lambertian reflectance model, which is obtained(12)

where x , [nx , ny , nz , a1 , a2 , . . . , ap ]T ∈ Rp+3 , andnx , ny , nz are the three elements of the surface normal.
Ã ∈ Rm×(3+p) is a data matrix whose j-th row is given byÃj =[−lxj , −lyj , −lzj , b1 − b0 ,While non-linear in I, g(I) is clearly linear in a =[a1 , . . . , ap ]T . Choices for each gk include polynomial, Gaussian, logistic, and spline functions as well as many others.
We choose to adopt a piecewise linear representation [13]which is composed of multiple polylinear functions of the form
(0 ≤ I < bk−1 )0
gk (I) = I − bk−1 (bk−1 ≤ I < bk )(10)
bk − bk−1 (bk ≤ I)

(11)

k=1

. . . , bk−1 − bk−2 , Ij − bk−1 , 0, . . . , 0].

(13)Here we assume bk−1 ≤ Ij < bk and lxj , lyj , lzj are threeelements of the j-th lighting direction.
Without loss of generality, we may avoid P
the degeneratepx = 0 solution to Eq. (12) by constraining
k=1 ak = 1.
For this purpose we replace Ã by A where [0, 0, 0, 1, . . . , 1] isappended as the (m+1)-th row of Ã and we define y ∈ Rm+1
as a vector of all zeros except for a one as the last element.
Given the appearance variations (I1 , I2 , . . . , Im ) under different known lighting conditions (l1 , l2 , . . . , lm ), the optimal
surface normal (n) and model parameters (a1 , a2 , . . . , ap ) arerecovered by solving the linear problem
min kW (Ax − y)k22 .
x(14)Here W , diag[1, · · · , 1, ∞] ∈ R(m+1)×(m+1) is a diagonal
weightingmatrix designed to strictly enforce the constraint
Ppa=
1. Note that 2 + p linearly independent images
k=1 k
are sufficient for producing a unique solution to Eq. (14). We
will refer to this photometric stereo method as piecewise linear

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013𝐼

least-squares regression (PL-LS). In spite of its simplicity, PLLS works for a wide variety of non-Lambertian diffuse materials (see our experimental results in Section 4.3). The problem
of course is that real images are frequently contaminated withvarious non-diffuse effects as modeled in Eq. (5). The next
section will focus on how to handle these corruptions within
our inverse diffuse reflectance model.

𝐼𝑖

n l=

p
X

ak gk (I) + e0 .

(15)

k=1

Note that although the value of e0 in Eq. (15) is different
from the e in Eq. (5), both factors are directed at the sameobservations and thus serve the same overall purpose (seeFig. 3). The standard least-squares based photometric stereo
problem corresponding to Eq. (15) would involve solving
min0 kW (Ax + e0 − y)k22
x,es.t. e0(m+1) = 0,(16)where e0 , [e01 , e02 , . . . , e0m+1 ]T ∈ Rm+1 . Note that e0(m+1)
Ppmust be set to zero to maintain the constraint k=1 ak = 1(see details in Section 3.4).
Given observations I and lighting directions l, the ultimate
goal is to recover surface normals and model parameters x
and non-diffuse corruptions e0 . However, this is an underconstrained problem since the number of unknowns p + 2 + m
will always exceeds the number of independent equations m.
One solution to this ambiguity is to apply simple
shadow/specular thresholding [9], [11] or a color channeltransformation [2] as a preprocessing step, to obtain an estimate of e0 and/or discard outliers. However, these types ofheuristics may discard useful information at times and comewith an additional computational expense. Moreover, graphbased approaches [24], [23], [42] and robust algorithms [3],[6] do not naturally embed within our framework since they
may conflict with our inverse piecewise diffuse model and/or
degrade the numerical stability.
In this paper, we instead introduce a special sparsity penalty
applied to e0 , whose minimization disambiguates the infinityof feasible solutions to Eq. (15). This penalty quantifies thereasonable observation that objects in the natural world exhibit
dominant diffuse reflections while non-diffuse effects emerge
primarily in limited areas of its appearance. For example,
specularities surround the spot where the surface normal isoriented halfway between lighting and viewing directions,while shadows are created only when lT n ≤ 0 (attachedshadow) or when a non-convex surface blocks the light (castshadow). Strictly speaking, we assume that the optimal feasible solution to Eq. (16) is acquired when the largest possible
number of observations are lying on the piecewise linear
diffuse reflectance function. Reflecting this assumption, ourestimation problem can be formulated asmin0 ke0 k0
x,e

s.t.

y = Ax + e0 , e0(m+1) = 0.

(17)

𝑓 𝒏𝑇 𝒍 = 𝑔−1 (𝐼)
𝑔(𝐼𝑖 )𝑒𝑖 ′𝑒𝑖𝑓(𝒏𝑇 𝒍𝑖 )3.3 Piecewise Linear Sparse Regression with NonDiffuse Corruptions
By introducing additive non-diffuse corruptions e0 in Eq. (11),the image formation model becomes
T

5

𝑇

𝒏 𝒍𝑖

𝒏𝑇 𝒍

Fig. 3. A 2-D point (nT li , Ii ) can be represented by boththe forward image formation model (Ii = f (nT li ) + ei )and the inverse image formation model (nT li = g(Ii )+e0i ).
The illustration indicates that they are uniquely convertible
if the reflectance function is monotonic.
Here, k · k0 represents the `0 -norm, which counts the number of non-zero entries in a vector. To reiterate, Eq. (17)builds on the assumption that images are captured underknown lighting conditions and any non-diffuse corruptions
have sparse structure. If these assumptions are not true (e.g.,
because of imperfect lighting calibration, no dominant diffuse
structure, etc.), then the hard constraint in Eq. (17) is no longer
appropriate. To compensate for more modeling errors, we relax
the hard constraint via an additional model mismatch penalty
givingmin0 kW (Ax + e0 − y)k2 + λke0 k0x,es.t. e0(m+1) = 0, (18)where λ is a nonnegative trade-off parameter balancing data
fit with sparsity. Note that in the limit as λ → 0, problems(17) and (18) are equivalent (the limit must be taken outsideof the minimization). Eq. (18) entails a difficult, combinatorial
optimization problem that must be efficiently solved at every
pixel. Here we consider two alternatives to brute force exhaustive search. First, in the machine learning and statistics
literature, it is common to replace the discontinuous, nonconvex `0 norm with the convex surrogate`1 norm. The `1Pnorm of a vector z is given by
|z
|,
which
constitutes
i ithe tightest convex approximation to the `0 norm. In certain
situations the resulting estimate will closely match the solution
to Eq. (18); however, in the context of photometric stereo thissubstitution may not always be adequate (see Section 3.5 formore details). Secondly, we can apply a simple hierarchical
Bayesian approximation to estimate x while simultaneously
accounting for e0 . This formulation, a generalized version ofsparse Bayesian learning (SBL) [43], is described in detail
next.
3.4

Recovery of Normals and Corruptions Via SBL

We assume the standard Gaussian likelihood function for the
first-level, diffuse errors givingp(y|x, e0 ) = N (y; Ax + e0 , λW −1 ),(19)IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

Note that we define W −1 , diag[1, · · · , 1, 0] ∈R(m+1)×(m+1) . We next apply an independent, zero-mean
Gaussian prior distributions on both x and e0 :p(x) = N (x; 0, Σx ),p(e0 ) = N (e0 ; 0, Γ).

While some details are omitted for brevity, using resultsfrom convex analysis it can be shown that for all u ≥ 0,
log |Σy | =log |Γ| + log AΣx AT + λW −1+ log Γ−1 + AΣx AT + λW −1(20)

Σx describes the prior variance of the unknown normal vector
n and model parameters a as Σx = diag(σn2 I(3) , σa2 I(p) )where I(k) ∈ Rk×k is the identity matrix; they are fixed
to convey our lack of a priori certainty about x. Thus
the prior on x will be relatively uninformative. The valuesof σn2 and σa2 will be discussed further below. In contrast,Γ , diag[γ] is a fully-parameterized, diagonal matrix, whereγ , [γ1 , . . . , γm+1 ]T is a non-negative vector of variancesin one-to-one correspondence with elements of e0 . A large
variance γi indicates that the corresponding e0i is free to reflect
the data, compensating for non-diffuse effects (outliers), whilea small or zero-valued variance implies that the associatederror term is constrained near zero (inliers). While we are
ignorant of which observations are outliers, γ(m+1) is fixed to
be zero because of the constraint regarding e0(m+1) in Eq. (18).
Combining the likelihood and prior using Bayes’ rule leads
to the posterior distribution p(x, e0 |y) ∝ p(y|x, e0 )p(x)p(e0 ).
To estimate x, we may further marginalize over e0 to giveZp(x|y) = p(x, e0 |y)de0 = N (x; µ, Σ),
(21)

6

−1

log |Γ| + log AΣx AT + λW −1
X ui− h∗ (u)+γiiX  ui≡+ log γi − h∗ (u),
γi
i

≤

(24)

where h∗ (z) denotes the concave conjugate function [45] of
−1and we haveh(β) , log diag[β] + AΣx AT + λW −1
removed irrelevant factors independent from u or γ. It can beshown that equality (and therefore the minimum of the righthand side) is obtained in Eq. (24) if and only ifh−1 i−1
u = diag Γ−1 + AΣx AT + λW −1.
(25)
In a somewhat related fashion, the second term in L(Γ) can
be upper-bounded via
−1Ty T Σ−1
AΣx AT + λW −1(y − z)y y ≤ (y − z)
X z2
i
,
(26)
+
γ
i
i
which holds for all z, with equality if and only ifwith mean and covariance defined as
−1µ = ΣAT Γ + λW −1
y,
h
i−1

T
−1 −1
.
Σ = Σ−1
+
A
Γ
+
λW
A
x

z = ΓΣ−1y y.
(22)

We now have a convenient closed-form estimator for x givenby the posterior mean. The only issue then is the values forthe unknown parameters Γ. Without prior knowledge as to thelocations of the sparse errors, the empirical Bayesian approachto learning Γ is to marginalize the full joint distribution over
all unobserved random variables, in this case x and e0 , andthen maximize the resulting likelihood function with respectto Γ [43]. Equivalently, we will minimizeZ
L(Γ) , − log p(y|x, e0 )p(x)p(e0 )dxde0≡ log |Σy | + y T Σ−1
y y
with

(23)
T

Σy , AΣx A + Γ + λW

−1,

with respect to Γ. While L(Γ) is non-convex, optimization
can be accomplished by adapting a majorization-minimization
approach from [44] to the photometric stereo problem. This
technique essentially involves the construction of rigorous
upper bounds (see below) on each of the two terms in
Eq. (23) using auxiliary variables z , [z1 , . . . , zm+1 ]T andu , [u1 , . . . , um+1 ]T (the EM algorithm can be viewed asa special case). For fixed values of z and u, a closed form
solution for Γ exists. Likewise, for a fixed value of Γ, theauxiliary variables can be updated in closed form to tighten
the upper bound around the current Γ estimate.

(27)

Now with z and u fixed, the overall upper bound on L(Γ)
decouples and we can solve for each γi individually bycollecting the γ-dependent terms from Eq. (24) and Eq. (26),leading to the problem
zi2 + ui
+ log γi ,
γi ≥0
γi
min

(28)

which has a simple closed-form solution.
Combining all of the above, and using matrix inversion formulae to produce numerically efficient computations, produces
update rules for the (k + 1)-th iteration given by

2
(k)
(k+1)
(k)
γi
←
zi
+ ui , ∀i, Γ(k+1) = diag[γ (k+1) ]−1z (k+1) ← Γ(k+1) S (k+1)
y
(29)


2 
−1 
u(k+1) ← diag Γ(k+1) − Γ(k+1)
S (k+1),
where S (k+1) is computed via

−1 T
T
S (k+1) = D − DA Σ−1
A D
x + A DAand

D , (Γ(k+1) + λW −1 )−1 .

(30)

These expressions only require O(n) computations and are
guaranteed to reduce L(Γ) until a fixed point Γ∗ is reached.
This value can then be plugged into Eq. (22) to estimate theunknown normal vector and model parameters. We denote thispoint estimator as xsbl . If the variances Γ∗ reflect the true
profile of the sparse errors, then xsbl will closely approximate

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

the true surface normal. This claim will be quantified more
explicitly in the next section.
We have thus far omitted details regarding the choice ofλ, σn2 and σa2 . The first one can be reasonably set according
to our prior expectations regarding the magnitudes of diffuse
modeling errors, but in practice there is considerable flexibility
here since some diffuse errors will be absorbed into e0 . In
contrast, we can realistically set σn2 → ∞, which implies zeroprecision or no prior information about the normal vectors
and yet still leads to stable, convergent update rules. However
we have observed that on certain problems a smaller selection
for σn2 can lead to a modest improvement in performance,
presumably because it has a regularizing effect that smoothes
out the cost function and improves the convergence path of
the update rules from Eq. (29). Interestingly and perhaps
counter-intuitively, in certain situations it does not alter the
globally optimal solution as discussed below. It is also possible
to learn σn2 using similar updates to those used for Γ, but
this introduces additional complexity and does not improve
performance.
Finally, we discuss the effects of σa2 on our algorithm. Since
the value of σa2 regularizes the shape of the piecewise linear
function Eq. (10), smaller values will prevent the reverse of
the sign of each linear segment and enforce the monotonicityof the function. While too small of a σa2 may limit thegenerality of g(I), it can be proven that as σa2 → 0 the model
naturally reduces to a simple, Lambertian form for the diffuse
component. We will empirically determine an appropriate
value for σa2 in Section 4.3.
3.5

Analytical Evaluation

Previously we discussed two tractable methods for solving
Eq. (18): a convex `1 -norm-based relaxation and a hierarchical
Bayesian model called SBL. This section briefly discusses
comparative theoretical properties of these approaches relevantto the photometric stereo problem. To facilitate the analysis,
here we consider the idealized case where there are no diffuse
modeling errors, or that λ is small. In this situation, the basic
estimation problem reverts back to Eq. (17).
If the lighting directions and sparse errors are in general
position, meaning they are not arranged in an adversarial configuration with zero Lebesgue measure, then it can be shownthat the minimizer of Eq. (17), denoted x0 , is guaranteed tobe the correct normal vector as long as the associated feasible
error component e0 = y − Ax0 satisfies ke0 k0 < m − p − 2.
Therefore, a relevant benchmark for comparing photometricstereo algorithms involves quantifying conditions whereby acandidate algorithm can correctly compute x0 .
In this context, recent theoretical results have demonstrated
that any minimizer x1 of the `1 relaxation approach will
equivalently be a minimizer of Eq. (17) provided ke0 k0 is
sufficiently small relative to a measure of the structure incolumns of the matrix A [46]. Unfortunately however, fortypical photometric stereo problems the requisite equivalency
conditions often do not hold (i.e., ke0 k0 is required to beprohibitively small) both because of structure imposed by thelighting geometry and implicit structure that emerges from the7

relatively small dimensionality of the problem (meaning we donot benefit from asymptotic large deviation bounds that apply
as m becomes large). Fortunately, SBL offers the potential for
improvement over `1 via the following result.
Theorem: For all σn2 , σa2 > 0 (and assuming λ → 0), if Γ∗ isa global minimum of Eq. (23), then the associated estimator
xsbl will be a global minimum of Eq. (17). Moreover, for σn2and σa2 sufficiently large it follows that: (i) Any analogous
locally minimizing SBL solution is achieved at an estimate
xsbl satisfying ky − Axsbl k0 ≤ m − p − 2, (ii) SBL canbe implemented with a tractable decent method such thatconvergence to a minimum (possibly local) that produces anxsbl estimator as good or better than the global `1 solution is
guaranteed, meaning ky − Axsbl k0 ≤ ky − Ax1 k0 .
The proof, which has been omitted for brevity, is relatively
straightforward and uses block-matrix inverse and determinant
identities, as well as ideas from [46], to extend SBL properties
derived in [47] to problems in the form of Eq. (17). We
may thus conclude that SBL can enjoy the same theoreticalguarantees as the `1 solution yet boosted by a huge potential
advantage assuming that we are able to find the global minimum of Eq. (23), which will always produce an xsbl = x0 ,unlike `1 .
There are at least two reasons why we might expect thisto be possible based on insights drawn from [47]. First, as
discussed previously, A will necessarily have some structure
unlike, for example, high dimensional random matrices. In
this environment SBL performance is often vastly superior to`1 because it can be shown to be implicitly based on an Adependent sparsity penalty that can compensate, at least inpart, for structure in A. Secondly, the sparse errors e0 will
likely have substantially different magnitudes depending on
image and object properties (meaning the non-zero elementsof e0 will not all have the same magnitude), and it has beenshown that in this condition SBL is more likely to converge
to the global minimum [47].

4

E XPERIMENTAL R ESULTS

In this section, we quantitatively evaluate our method on
synthetic and real image data. All experiments were performedon an Intel Core2 Duo E6400 (2.13GHz, single thread) machine with 4GB RAM and were implemented in MATLAB.
For the SBL- and `1 -based methods we used λ = 1.0−6 inthe synthetic experiments with no additive image noise andperfect lighting calibrations (Section 4.2 and Section 4.3), and
λ = 10−2 for the other cases (Section 4.4). We set σn2 = 106
for all experiments. As for σa2 , which can affect the solution
when p > 1, we experimentally find the optimal value inSection 4.3-(a).
4.1Datasets

For quantitatively evaluating our method, four different
datasets are used. We generate 32-bit HDR images of twotarget scenes, Bunny (256 × 256) and Caesar (300 × 400),
illuminated under random directional lightings. We use a few

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013(a) Ground Truth(c) SBL(d) L18(e) R-PCA(f) LS1.00.0(g) SBL(b) Input(h) L1(i) R-PCA(j) LS

Fig. 4. Recovery of surface normals from 40 images of Caesar (300 × 400) in dataset (A) with explicit shadow removal.
(a) Input, (b) Ground truth, (c)-(j) Recovered surface normals and Error maps (in degrees).
TABLE 1
Experimental results of Bunny (left) / Caesar (right) in dataset (A) with varying number of images.
No. ofimagesMean error (in degrees)
L1 R-PCA LSSBLMedian error (in degrees)
Elapsed time (sec)
L1 R-PCA LS SBL L1 R-PCA LS

SBL

Mean error (in degrees)
L1 R-PCA LSSBLMedian error (in degrees)
L1 R-PCA LS

SBL

Elapsed time (sec)
L1 R-PCA LS

SBL

6.2
6.0 26.2 7.4 4.73 4.63 31.0 4.97 106.4 34.2 45.8 15.9
5
6.0
6.0 15.3 7.0 4.7
4.3 10.7 4.8 46.5 13.6 15.7 4.6
10 0.09 0.61 3.8 1.9 0.27 0.58 0.81 1.8 36.3 13.6 37.8 5.9 0.24 0.40 10.7 0.94 0.19 0.26 1.8 0.93 97.2 34.1 93.7 19.2
15 0.076 0.16 0.21 1.6 0.052 0.13 0.19 1.6 26.8 13.1 55.1 6.3 0.044 0.11 2.6 0.77 0.047 0.083 0.14 0.76 67.0 31.4 153.3 22.0
20 0.033 0.080 0.11 1.6 0.022 0.078 0.11 1.6 24.2 13.5 70.5 6.9 0.018 0.051 0.079 0.76 0.015 0.035 0.065 0.72 60.9 32.7 177.5 23.3
25 0.018 0.055 0.084 1.6 0.010 0.048 0.069 1.6 23.1 14.1 86.0 7.6 0.011 0.034 0.068 0.76 0.0081 0.023 0.059 0.76 57.9 34.3 196.9 25.1
30 0.012 0.037 0.080 1.7 0.0048 0.032 0.065 1.7 23.1 14.2 121.0 8.4 0.0063 0.018 0.043 0.77 0.0082 0.018 0.043 0.77 58.1 33.2 231.7 27.7
35 0.0057 0.023 0.098 1.6 0.0029 0.019 0.093 1.6 22.7 14.6 161.3 8.5 0.0045 0.012 0.036 0.78 0.0031 0.0084 0.033 0.80 58.4 34.5 259.2 29.4
40 0.0039 0.019 0.12 1.6 0.0020 0.015 0.12 1.6 22.6 15.0 200.7 9.4 0.0031 0.0094 0.037 0.76 0.0019 0.0063 0.034 0.78 59.6 35.2 281.2 31.5

different BRDF settings for rendering; (A) combination ofLambertian diffuse reflection and specularity of the CookTorrance reflectance model [48], (B) combination of Lafortune
diffuse reflection [49] and Cook-Torrance specularity, and (C)MERL BRDF database [41]. Additionally, as the forth dataset,
denoted (D), we record real images for qualitatively evaluating
our method in a practical scenario. For the datasets (A), (B),and (C), both cast shadows and attached shadows are also
rendered. Note that though light calibration is noise-free, andshadowed pixel intensities are exactly zero in the synthesizeddatasets, some calibration errors and non-zero shadowed pixels
exist in the real data (D).
4.2

Performance Evaluation with p = 1

We begin with the simplest case where the number of basis
functions in the piecewise linear model is one i.e., the diffuse
function is assumed to be Lambertian. We evaluate performance using the synthesized images in dataset (A), which
has Lambertian diffuse reflections along with non-Lambertian

specularities and shadows. Here we compare our SBL-basedand convex `1 -norm-based sparse regression models from
Sections 3.3 and 3.4 with the recent R-PCA-based method
proposed by Wu et al. [6], which attempts the simultaneousrecovery of a low-rank Lambertian factor and sparse corruptions. We also compare with the standard least squares(LS) Lambertian photometric stereo estimator [1], which isequivalent to solving Eq. (14) with p = 1.
We change experimental conditions with regard to thenumber of images, surface roughness (i.e., the ratio ofspecularities), shadow removal (i.e., whether or not a shadow
mask is used to remove zero-valued elements from theobserved images), and the presence of additional Gaussian
noise. Note that when in use as defined for each experiment,
the shadow mask is applied equivalently to all algorithms.
To increase statistical reliability, all experimental results areaveraged over 20 different sets of 40 input images. While
the amount of cast/attached shadows are fixed to 24.0%
and 27.8% of the number of all the pixels, respectively, theIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013(a) Ground Truth(c) SBL(d) L1(e) R-PCA

9(f) LS
1.0(b) Input(g) SBL(h) L1(i) R-PCA(j) LS

0.0

Fig. 5. Recovery of surface normals from 40 images of Bunny(256×256) in dataset (A) without explicit shadow removal.
(a) Input, (b) Ground truth, (c)-(j) Recovered surface normals and Error maps (in degrees).

(a) Ground Truth(c) SBL(d) L1(e) R-PCA(f) LS
1.0(b) Input(g) SBL(h) L1(i) R-PCA(j) LS

0.0

Fig. 6. Recovery of surface normals from 40 images of Bunny in dataset (A) with explicit shadow removal and additive
Gaussian noises (30%). (a) Input, (b) Ground truth, (c)-(j) Recovered surface normals and Error maps (in degrees).

percentage of specularities depends on the surface roughness
parameter that is independently selected in each experiment.
To quantitatively evaluate the performance, we compute theangular error between the recovered normal map and theground truth.
Valid number of images for efficient recovery in thepresence of specularities
In this experiment, we vary the number of images to estimate
the minimum number required for effective recovery whenusing the shadow mask with fixed surface roughness. The
percentage of specular pixels in Bunny and Caesar are 8.4%
and 11.6%, respectively.
Once 40 images are generated for each dataset, the image
subset is randomly subsampled without replacement. The
results are displayed in Table 1 and Fig. 4. We observe thatthe sparse-regression-based methods are significantly moreaccurate than both R-PCA and LS. It is also clear that SBL
is more accurate than `1 , although somewhat more expensive
computationally.1 Note that, although not feasible in general,
when the number of images is only 5, the most accurate
and efficient implementation for regression could be to justsystematically test every subset of 3 images (i.e., brute force
search only requires 10 iterations at each pixel).
Robustness to shadows and image noise
We now evaluate the robustness of our method against
corruptions; shadows and image noise. We set two conditions
for evaluating the effects of (i) shadows (fixed specularities, no
shadow removal, no image noise), (ii) additive Gaussian image
1. The SBL convergence rate can be slower with fewer images because of
the increased problem difficulty. This explains why the computation time may
actually be shorter with more images.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

TABLE 3
Experimental results of Bunny in dataset (A) with varying
amount of additive Gaussian noises.

TABLE 2
Results of Bunny in dataset (A) without shadow removal.
Mean error (in degrees)
No. ofimages SBL L1 R-PCA LS

5
10
15
20
25
30
35
40

5.2 11.9
2.8 5.6
1.9 4.0
1.2 2.7
0.81 1.9
0.62 1.6
0.59 1.5
0.53 1.2

Median error (in degrees)
L1 R-PCA LS

SBL

10

Elapsed time (sec)
L1 R-PCA LS

SBL

Mean error (in degrees)
Dens. ofnoises (%) SBL
L1 R-PCA LS
100.0079 0.040 0.16 3.3
20
0.021
0.79 4.4
0.11
30
0.068
3.6
5.3
0.29
40
0.21
9.8
6.2
0.70
50
0.58
11.7 7.0
1.5

12.1 12.1 5.0 12.3 12.5 12.5 213.0 37.0 45.8 5.1
10.9 10.9 2.3 5.6 11.3 11.3 98.9 33.0 93.7 6.0
9.9 10.0 2.3 4.0 10.1 10.2 66.8 32.5 153.3 7.4
9.4 9.5 1.0 2.7 9.6 9.6 52.9 30.0 177.5 7.6
8.9 9.0 0.69 1.8 8.9 9.0 46.2 31.0 196.9 9.1
9.0 9.1 0.61 1.5 8.9 8.9 41.1 32.0 231.7 9.4
9.1 9.1 0.58 1.4 9.3 9.3 41.1 34.4 259.2 11.0
8.8 8.9 0.58 1.2 9.0 9.1 39.4 33.3 281.2 10.7

Median error (in degrees)SBL
L1 R-PCA LS
0.0060 0.039 0.16 3.3
0.019 0.099 0.80 4.3
0.060
0.25
3.2
5.2
0.18
0.63
9.9
6.1
0.53
1.4
11.7 6.9

5

10

0

0
(a)

(b)

(c)

(d)

(e)

Fig. 7. Comparison between SBL and `1 -based method. Error maps of (a) SBL and (b) L1 (in degrees). The per-pixel
number of (c) specularities, (d) shadows, (e) corruptions (The maximum is 5). These maps are illustrated based on
one of twenty datasets.

The results are illustrated in Fig. 5, Fig. 6, Table 2 andTable 3. While performance of each method degrades when additional corrupted pixels (outliers) contaminate the estimation
process, our sparse regression methods outperform both RPCA and LS in accuracy and outperform R-PCA in efficiency
in the presence of shadows and noise. We observe that theperformance of R-PCA degrades when a shadow mask isunavailable (i.e., the positions of missing entries are unknown)
while our sparse-regression based method automatically compensates for missing entries in the estimation process.
We also observe that SBL is more accurate than `1 in allconditions with slight more expense of computational effort.
We further compare the results of SBL and `1 using the casewhere the number of images is 5 without removing shadows.
The error maps and the numbers of corruptions per-pixel are
displayed in Fig. 7. We observe that the `1 method typically
suffer from shadows while SBL can find the correct solution
in most pixels as long as the number of corruptions is lessthan 3, which is the theoretical limit when only 5 observations
are given.
Sparseness of corruptions
While it does not often occur, our assumption that corruptions

10

Mean angular errors(in degrees)

noise (fixed specularities, explicit shadow removal, varying
amount of image noise). The ability to estimate surface
normals without an explicit shadow mask is important, sincein practical situations shadow locations are not always easy
to be determined a priori. The number of images is 40 in (ii)and varying from 5 to 40 in (i). We use Bunny for evaluation,and the ratio of specularities is 8.4% in (i) and (ii), and image
noise is 10% to 50% in (ii). Image noise obeys a zero-mean
iid Gaussian distribution with σ 2 = 0.1.

SBL
L1
R-PCA
LS

8
6
4

2
0
0

10

20

30

40

50

60

% of specularities
Fig. 8. Experimental results of Bunny in dataset (A) withvarying amount of specularities. The x-axis and y-axis
indicate the ratio of specularities and the mean angularerror of normal map.

appear sparsely over images can be violated in somesituations, e.g., observations include wide-lobe specularities.
In this experiment, we examine how our approach is affectedby the sparseness of corruptions by changing the width of
the specular lobe (i.e., percentage of specularities appearedin observations). Here, we use Bunny dataset whose amountof specularities is varying from 10% to 60%, respectively, by
appropriately choosing the model parameter of Cook-Torrance
reflectance model. In this experiment, we remove shadows
explicitly to factor out their effect.
The result is illustrated in Fig. 8. As expected, performanceof both SBL and `1 degrades as the ratio of specular
corruptions increases. However, even when the sparsenessof outliers is difficult to be held, we observe that ourIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013sparse-regression based method still outperforms R-PCA andLS.
Comparison with RANSAC-based approaches
We also compared our approach with two different RANSACbased methods. Here we use 40 images of Bunny dataset
whose percentage of specular pixels is 8.4%. Shadows are
removed and images noises are not added in this experiment.
The first one is a photometric linearization approach [3],which converts an input image into the ideal Lambertian
image. It uses RANSAC for robustly identifying the basis
images in the presence of non-Lambertian corruptions. Oncelinearized images are acquired, the standard Lambertian-based
approach [1] is applied with known lightings to estimate
surface normals. The second one uses RANSAC more directly in Lambertian photometric stereo method like [21]. In
this approach, three images are randomly sampled for eachposition independently, and the surface normal and albedo are
estimated from them using Lambertian photometric stereo [1].
Then the Euclidean distances between observations and predicted intensities are computed, and the number of inliers
whose distances are less than a threshold are counted. Afterthe sampling process, the surface normal and albedo with the
maximum number of inliers are adopted.
The results are summarized in Table 4 and Table 5, respectively. We have also included standard deviations in the tables
for observing the estimation stability. We observe in Table 4that although we use a large number of samples for RANSAC(e.g., 2000) for the photometric linearization, this approachcannot always stably find the solution especially when thenumber of images is large. On the other hand, our method
is able to successfully and efficiently finds the solution. We
also observe in Table 5 that a large number of samples (e.g.,200) gives very accurate estimation while it takes much time
to compute it. The smaller number of samples gives more
efficient but less accurate estimation. When the number ofsamples is 20, the computational cost is almost same as our
SBL-based method, but the mean error is 50 times larger
than ours, which also demonstrates the effectiveness of the
proposed method.
It is also crucial to emphasize that the RANSAC algorithm
scales very poorly as the number of inlier model parametersis increased. Simply put, a larger number of inlier random
samples is required to robustly estimate a larger numberof model parameters, and obtaining such inlier samples becomes combinatorially more difficult in higher dimensions.
Consequently, while RANSAC may work reasonably wellhere recovering surface normals under a Lambertian diffusivemodel (where the number of unknowns is effectively only 3), it
will become intractable when using more complex reflectancemodels, which is a central concern herein. In contrast, ourmethod scales linearly in the number of inlier parameters andcan therefore be robustly expanded to handle non-Lambertianparameterized reflectance functions with outliers, e.g., the p >1 model detailed in Section 4.3 or other more sophisticatedextensions.

11

TABLE 4Comparison with RANSAC based approach [3] with
Bunny in dataset (A) (Number of samples is 2000).
Mean error
[3]

SBL

6.7
0.74
0.61
0.70
1.0
2.3
3.2
2.1

4.7
0.27
0.052
0.022
0.010
0.00480.00290.0020No. ofImages

SBL

5
10
15
20
25
30
35
40

6.0
0.09
0.076
0.033
0.018
0.012
0.00570.0039Median error
[3]

5.4
0.38
0.12
0.058
0.063
0.042
0.051
0.046

Elapsed time[3]Standard deviation
[3]
SBL

SBL

4.4
1.4
1.9
1.2
2.6
4.0
9.2
4.3

46.5
36.3
26.8
24.2
23.1
23.1
22.7
22.6

4.1
0.35
0.059
0.027
0.020
0.012
0.00640.0048

52.3
544.0
958.4
1048.91141.81227.81327.91430.4

TABLE 5
Performance of RANSAC based approach [21] with
Bunny in dataset (A) (Number of images is 40).

No. of Samples Mean error

200
100
80
60
40
20
10

4.3

0.000019
0.00350.0099
0.019
0.034
0.22
0.51

Median error

7.9E-06
7.6E-06
7.6E-06
7.6E-06
7.3E-06
7.4E-06
7.5E-06

Standard
deviation
0.000027
0.31
0.70
0.98
1.2
3.4
5.0

Elapsed
time
283.0
142.4
115.2
86.0
57.9
29.8
15.6

Performance Evaluation with p > 1

If the majority of observations are represented by non-linear
diffuse reflections, then piecewise linear sparse regression withp > 1 basis functions or segments is expected to be considerably more effective. In this experiment, we evaluate the pfunctions piecewise sparse linear regression to those complex
objects by using the dataset (B) rendered with non-linear
Lafortune diffuse reflectance model [49] and Cook-Torrance
Reflectance model, and (C) rendered with one hundred BRDF
functions from the MERL BRDF database [41].
We generate 32-bit HDR images of Bunny (256 × 256)
scene with Lafortune diffuse reflectance and Torrance sparrow
specularities under 40 different lightings in dataset (B), and100 different lightings for 100 different materials in dataset(C). Note that the average percentage of specularities in dataset(B) is 8.4% and 24.0% of cast/attached shadows are renderedin both dataset (B) and dataset (C).
Here, in addition to the SBL-based piecewise linear sparse
regression (PL-SBL), we further implement the method which
is solving Eq. (14) (PL-LS). We compare our methods with
the R-PCA-based method [6] and a recent parametric nonLambertian photometric stereo method with biquadratic reflectance model [11], which reasonably represents the lowfrequency of non-linear reflectance though shadows and highfrequency observations must be removed in advance. The

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

I = α1 (nT l)2 (lT h)2 + α2 (nT l)2 (lT h) + α3 (nT l)2(31)+ α7 (lT h)2 + α8 (lT h) + α9 ,where h is a half vector directed at l + v andα , [α1 , α2 , . . . , α9 ] are model parameters. In theimplementation of [11], we used the non-shadowed pixels
whose intensities are ranked below the 25%, and both
surface normals and model parameters are iteratively updated
initialized by the Lambertian photometric stereo method. Note
that, shadowed pixels are also rejected, but all observations
are taken into account in other methods.
Choosing σa2 and p
Here we examine the choices for σa2 and for the numberof basis functions p in piecewise linear model by using
dataset (C). Note that if p is large enough, we can essentiallyrepresent any complex non-linear function, although we run
the risk of over-fitting when m is too small. However, we canalways compensate to some extent via σa2 , since regardlessof p, for σa2 sufficiently small we approximate the standard
Lambertian model with all linear segments having equal
slope, and therefore equal diffuse albedo.
In this experiment, we vary individually σa2 and p to findthe optimal parameters. The results are illustrated in Fig. 9.
Here, average angular errors in each normal map are further
averaged over 100 materials. We observe that as expected, too
small σa2 deteriorates the performance, yet does not affect theperformance when the value is sufficiently large. Therefore,
in the following experiments, we fix σa2 by 1.0. As for thenumber of basis functions, it appears that p = 3 is optimal in
the case of dataset (C).
We also illustrate per-material angular errors with different
p in Fig. 10. We observe that piecewise linear function withmany basis functions, e.g., 6, works very well under thematerial with narrow specular peaks, e.g., for (1) specularwhite-phenolic and (4) gray-plastic, while it suffers fromover-fitting in polished metal with the broad specularity, e.g.,for (82) silver-metallic paint or materials with complex 2-lobe
BRDF, e.g., (53) natural-209 and (83) ipswich-pine-221),
where we can hardly see the underlying monotone diffuse
structure in observations.
Quantitative comparison with other methods
We evaluate of the performance of our method by a numericalcomparison with other methods using datasets (B) and (C).
Here, our methods (PL-SBL, PL-LS) are compared with thebiquadratic photometric stereo method by Shi et al. [11](Biquadratic), LS and R-PCA (R-PCA is included only fordataset (C)).
First, we use dataset (B) to verify that our inverse piecewise linear diffuse model reasonably handles non-Lambertiandiffuse reflections. In the Lafortune model [49], a generalrotationally symmetric diffuse component ρd is written as
ρd = (nT l)k (nT v)k . Here k is a model parameter which
is fixed to 3.0 in our experiment. Note that specularities ofMean angular errors(in degrees)

10biquadratic model is a approximated BRDF as follow,+ α4 (nT l)(lT h)2 + α5 (nT l)(lT h) + α6 (nT l)

12

8
6
4

2

p=1
p=2
p=3
p=4
p=5
p=6

0
1.00E-06 1.00E-04 1.00E-02 1.00E+00 1.00E+02 1.00E+04 1.00E+06

The value of 𝜎𝑎Fig. 9. Experimental results of dataset (B) with different
values of p and σa2 . The results are averaged over 100
different materials.

dataset (B) are also rendered with Cook-Torrance model inthe same manner with dataset (A). Therefore we can consider
dataset (B) as a combination of dominant general diffuse
reflections and sparse specular reflections.
The results are illustrated in Fig. 11 and Table 6. We present
the results of PL-SBL with p = 1 and p = 3, and PL-LS with
p = 3. Note that PL-LS with p = 1 is exactly same with
LS. We observe that while non-Lambertian diffuse reflections
degrade the performances of LS and PL-SBL(p = 1), Biquadratic, PL-LS(p = 3) and PL-SBL(p = 3) work better since
they are potentially capable of non-Lambertian reflections.
We also observe that PL-SBL(p = 3) performs best since
the estimation of Biquadratic and PL-LS are disrupted byspecularities which are not included in each model while thesparsity penalty in PL-SBL(p = 3) reasonably neglects themas the model outlier.
Secondly, we evaluated the performance of our method byusing datsaet (C). The results are illustrated in Fig. 12. We
observe that while SBL(p = 1) outperforms R-PCA and LS
in most of materials, the average angular errors are large forvarious complex materials. In contrast, PL-SBL(p = 3) works
better for those kind of objects since the piecewise linear function with several basis functions can capture the nonlinearityof the non-Lambertian diffuse structure. We also observe thatthe Biquadratic model is the most effective method of all fordataset (C), however we emphasize that our method achieves
competitive performance even without heuristically filteringhigh-frequency specularities. Consequently, on real images ourapproach produces qualitatively superior results (see Section
4.4 below). Finally, we can see from the comparison between
PL-LS(p = 3) and PL-SBL(p = 3) that the sparsity penalty
is also helpful for large p in our method, which can reject
specularities efficiently.
We also compare our method with others in the presenceof shadows and image noises. The results are illustrated inTable 7. We observe that our SBL-based methods still work
under attached/cast shadows and image noises in contrast
that Biquadratic and least-squares-based methods are easily
disrupted by those corruptions.

35

30

25

20

15

10

Mean errors (in degrees)yellow-paintspecular-white-phenolicorange-paintwhite-paintgreen-plastic
red-specular-plastic
specular-yellow-phenolicpink-jasper
white-fabricpure-rubber
white-diffuse-bballwhite-acrylic
specular-orange-phenolic
specular-red-phenolic
specular-maroon-phenoliclight-red-paintspecular-violet-phenolic
alumina-oxide
polyurethane-foam
purple-paintpink-fabricgray-plastic
dark-red-paintneoprene-rubber
white-marble
specular-blue-phenolicmaroon-plastic
specular-green-phenolicpink-plastic
teflongreen-acrylic
fruitwood-241
red-plastic
violet-rubber
blue-acrylic
yellow-matte-plastic
white-fabric2
red-fabric2
blue-rubberdelringreen-fabricred-phenolic
aventurine
yellow-plastic
nylon
pearl-paintpvc
pink-felt
red-fabriccolonial-maple-223
polyethylene
gold-paintdark-blue-paintnatural-209
cherry-235
pickled-oak-260
silicon-nitrade
light-brown-fabric
blue-fabricviolet-acrylic
yellow-phenolic
green-latex
silver-paintcolor-changing-paint2
gold-metallic-paint2
special-walnut-224
dark-specular-fabricgreen-metallic-paint2
two-layer-gold
gold-metallic-paint3
black-phenolic
black-fabricspecular-black-phenolicpink-fabric2
black-oxidized-steel
red-metallic-paintblack-soft-plastic
green-metallic-paintblue-metallic-paint2
two-layer-silver
black-obsidianalum-bronze
silver-metallic-paintipswich-pine-221
gold-metallic-paintcolor-changing-paint3
chromecolor-changing-paint1
hematitesilver-metallic-paint2
beige-fabricsteel
aluminum
nickelchrome-steel
blue-metallic-painttungsten-carbide
ss440
brass
grease-covered-steel

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 20130(b) Input

p=1
p=2
p=3
p=4
p=5
p=6

10

(a) Ground Truth

20(c) LS(h) LS

0
10
20
30
40

50

60
70
80
90

30
40

(d) BQ(i) BQ5060(e) PL-SBL (p=1)(j) PL-SBL (p=1)

13

40

5

0
70

(f) PL-LS (p=3)(k) PL-LS (p=3)

8090Fig. 10. Experimental results of dataset (C) under the different number of basis functions [41]. We aligned results inascending order of mean angular error of PL-SBL(p = 3). Some rendered images corresponding to each material ID
are also illustrated.

(g) PL-SBL (p=3)10.0(l) PL-SBL (p=3)

0.0

Fig. 11. Recovery of surface normals from 40 images of Bunny in dataset (B) with explicit shadow removal. (a) Input,(b) Ground truth, (c)-(l) Recovered surface normals and Error maps (in degrees).

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

14

40
LS

Mean errors (in degrees)

35
30
25
20
15
10

0

10

20

30

40

50

60

70

80

90

SBL (p=1)
R-PCA
Biquadratic
PL-LS (p=3)
PL-SBL (p=3)

5

0
10

20

30

40

50

60

70

80

90

yellow-paintspecular-white-phenolicorange-paintwhite-paintgreen-plastic
red-specular-plastic
specular-yellow-phenolicpink-jasper
white-fabricpure-rubber
white-diffuse-bballwhite-acrylic
specular-orange-phenolic
specular-red-phenolic
specular-maroon-phenoliclight-red-paintspecular-violet-phenolic
alumina-oxide
polyurethane-foam
purple-paintpink-fabricgray-plastic
dark-red-paintneoprene-rubber
white-marble
specular-blue-phenolicmaroon-plastic
specular-green-phenolicpink-plastic
teflongreen-acrylic
fruitwood-241
red-plastic
violet-rubber
blue-acrylic
yellow-matte-plastic
white-fabric2
red-fabric2
blue-rubberdelringreen-fabricred-phenolic
aventurine
yellow-plastic
nylon
pearl-paintpvc
pink-felt
red-fabriccolonial-maple-223
polyethylene
gold-paintdark-blue-paintnatural-209
cherry-235
pickled-oak-260
silicon-nitrade
light-brown-fabric
blue-fabricviolet-acrylic
yellow-phenolic
green-latex
silver-paintcolor-changing-paint2
gold-metallic-paint2
special-walnut-224
dark-specular-fabricgreen-metallic-paint2
two-layer-gold
gold-metallic-paint3
black-phenolic
black-fabricspecular-black-phenolicpink-fabric2
black-oxidized-steel
red-metallic-paintblack-soft-plastic
green-metallic-paintblue-metallic-paint2
two-layer-silver
black-obsidianalum-bronze
silver-metallic-paintipswich-pine-221
gold-metallic-paintcolor-changing-paint3
chromecolor-changing-paint1
hematitesilver-metallic-paint2
beige-fabricsteel
aluminum
nickelchrome-steel
blue-metallic-painttungsten-carbide
ss440
brass
grease-covered-steel0

Fig. 12. Comparison among different methods with dataset (C). We aligned results in ascending order of mean angularerror of PL-SBL(p = 3). Some rendered images corresponding to each material ID are also illustrated.

TABLE 6Comparison among different methods with dataset (B).
No. ofImages

LS

5
10
15
20
25
30
35
40

4.9
3.7
3.7
3.8
3.8
3.7
3.8
3.8

Mean error (in degrees)
PL-SBL PL-LS
BQ
(p=1)
(p=3)
6.9
5.4
25.1
4.4
3.3
2.6
3.5
3.5
1.4
2.8
3.6
1.3
2.2
3.6
1.3
1.9
3.5
1.2
1.7
3.6
1.2
1.6
3.7
1.2

PL-SBL
(p=3)
6.0
1.5
0.80
0.67
0.58
0.49
0.47
0.49

TABLE 7
Experimental results of dataset (C) under the different
kind of corruptions. This table illustrates average angularerrors of surface normals for 100 materials.
Mean error (in degrees)
PL-LS
LS SBL (p=1) Biquadratic(p=3)
w/o shadow and noise 13.9
w/ shadow
w/ noise

14.3
35.8

8.9
7.7
8.2

2.9
44.0
11.0

4.4
28.0
14.1

PL-SBL
(p=3)

4.1
10.5
4.8

4.4 Qualitative Evaluation with Real Images
We also evaluate our algorithm (PL-SBL) using real images(dataset (D)). We captured RAW images without gamma correction by Canon 30D camera with a 200[mm] tele-photo lensand set it 1.5[m] far from target object. Lighting conditions
are randomly selected from a hemisphere whose radius is
1.5[m] with the object placed at the center. For calibrating
light sources, a glossy sphere was placed in the scene. We
use a set of 25 images of Chocolate bear (261 × 421), and 40
images each of Doraemon (269×420) and Fat guy (293×344).
Note that in this experiment, we did not remove shadows fromimages by zero-intensity thresholding since even shadowed
pixels have non-zero values due to the presence of slight ambient illumination, sensor satiations, low signal-to-noise ratio,
inter-reflections between the object and the floor and so on.
Therefore, we sort all observations in increasing order and use
lowest 25% of observations for the low-intensity condition of
Biquadratic (low). We also applied the biquadratic reflectancemodel using all observations denoted as Biquadratic (all). We
evaluate the performance of PL-SBL by visual inspection of
the output normal maps and recovered surface meshes by apoisson solver [50], which recovers surface meshes from a
gradient map.
The results are illustrated in Fig. 13. First, we observe that
PL-SBL(p = 1, p = 3) succeeds in efficiently rejecting specularities and estimating smoother and more reasonable normal
maps. While Biquadratic (low) has produced the highest
performance in Section 4.4(b), it has substantial difficulty with

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

real images where calibration errors, shadowed pixels, andsensor saturations are often mis-classified as low-frequency
reflectance, leading to unpredictable errors in practice (see Table 7). We note that increasing the threshold for shadow
removal might solve this problem, but the optimal threshold selection remains hard and scene-dependent. Although
including higher frequency observations may sometimes give
better results by diluting low-frequency corruptions with alarger number of samples (see Biquadratic (all)), accurate
estimation still remains hard since the simple biquadraticmodel has difficulty representing complex, non-linear highfrequency observations (as is mentioned in [11]).
Second, we observe that reconstructed surface meshes by
PL-SBL(p = 3) are more reasonable than PL-SBL(p = 1) e.g.,
the stomach of fatguy recovered by PL-SBL(p = 1) is shaper
than that of PL-SBL(p = 3), though we can not make thefurther quantitative comparison due to the lack of the groundtruth. From those observations, we can say that our method iseffective, especially when p > 1, in the practical situation.

5

C ONCLUSION

In this paper, we have proposed the inverse piecewise linear
sparse regression based photometric stereo method whichworks under various kind of non-diffuse corruptions. At first,
we presented a special case of Lambertian-based sparse regression which consider all kind of non-Lambertian corruptions
as outliers. Our detailed discussion and experimental results
has shown that our method are not only robust to thosekind of corruptions, but also computationally more efficientand stable comparing than other robust approaches. However,since Lambertian sparse regression could not handle nonLambertian materials, we have secondly discussed the generalcase where the piecewise linear function has several basis
functions. Our extensive evaluation with both synthetic andreal images indicated that our method works not only forvarious non-diffuse corruptions, but also real images where
the shadow removal is not easy.
Our current limitation is that we assume the existence of theunique dominant diffuse structure in observations. Therefore
our method cannot potentially handle materials with unknown
dominant structure e.g., materials with unknown BRDF where
dominant structure is specularities or 2-lobe BRDF which isrepresented by a sum of several functions whose dominantdirections are unknown polyethylene and natural ones likefabrics. We are very interested in extending our work to handle
those complex materials. The other interesting direction of our
future work is applying our method to the uncalibrated casesince we currently assume that the lighting environment iscalibrated. Since the simple structure of our image formation
model have the possibility to be combined with other known
constraints such as isotropy symmetry, it seems that the realization of the uncalibrated photometric stereo method forgeneral diffuse materials is not impossible.

ACKNOWLEDGMENTS
This work was supported by the Grants-in-Aid for JSPS
Fellows (248615).

15

R EFERENCES[1]

P. Woodham, “Photometric method for determining surface orientationfrom multiple images,” Opt. Engg, vol. 19, no. 1, pp. 139–144, 1980.
[2] S. P. Mallick, T. E. Zickler, D. J. Kriegman, and P. N. Belhumeur,“Beyond lambert: reconstructing specular surfaces using color,” Proc.
CVPR, 2005.
[3] Y. Mukaigawa, Y. Ishii, and T. Shakunaga, “Analysis of photometricfactors based on photometric linearization,” Journal of the Optical
Soceity of America, vol. 24, no. 10, pp. 3326–3334, 2007.
[4] C. Yu, Y. Seo, and S. Lee, “Photometric stereo from maximum feasible
lambertian reflections,” in Proc. ECCV, 2010.
[5] K. Sunkavalli, T. Zickler, and H. Pfister, “Visibility subspaces: Uncalibrated photometric stereo with shadows,” Proc. ECCV, 2010.
[6] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma, “Robust
photometric stereo via low-rank matrix completion and recovery,” inProc. ACCV, 2010.
[7] A. Georghiades, “Incorporating the torrance and sparrow model ofreflectance in uncalibrated photometric stereo,” Proc. ICCV, 2003.
[8] H. Chung and J. Jia, “Efficient photometric stereo on glossy surfaces
with wide specular lobes,” Proc. CVPR, 2008.
[9] D. Goldman, B. Curless, A. Hertzmann, and S. Seitz, “Shape andspatially-varying brdfs from photometric stereo,” in Proc. ICCV, October
2005.
[10] N. Alldrin, T. Zickler, and D. Kriegman, “Photometric stereo with nonparametric and spatially-varying reflectance,” in Proc. CVPR, 2008.
[11] B. Shi, P. Tan, Y. Matsushita, and K. Ikeuchi, “A biquadratic reflectancemodel for radiometric image analysis,” Proc. CVPR, 2012.
[12] ——, “Elevation angle from reflectance monotonicity,” in Proc. ECCV,
2012.
[13] E. Vieth, “Fitting piecewise linear regression functions to biologicalresponses,” Journal of Applied Physiology, vol. 67, no. 1, pp. 390–396,
1989.
[14] S. Ikehata, D. Wipf, Y. Matsushita, and K. Aizawa, “Robust photometricstereo using sparse regression,” Proc. CVPR, 2012.
[15] S. Barsky, “The 4-source photometri stereo technique for threedimensional surfaces in the presence of highlights and shadows,” IEEE
Trans. Pattern Anal. Mach. Intell., vol. 25, no. 10, pp. 1239–1252, 2003.
[16] V. Argyriou and M. Petrou, “Recursive photometric stereo when multiple
shadows and highlights are present,” Proc. CVPR, 2008.
[17] J. Ackermann, F. Langguth, S. Fuhrmann, and M. Goesele, “Photometricstereo for outdoor webcams,” Proc. CVPR, 2012.
[18] M. Chandraker, S. Agarwal, and D. Kriegman, “Shadowcuts: Photometric stereo with shadows,” Proc. CVPR, 2007.
[19] M. Liao, X. Huang, and R. Yang, “Interreflection removal for photometric stereo by using spectrum-dependent albedo,” Proc. CVPR, 2011.
[20] T. Okatani and K. Deguchi, “Optimal integration of photometric and geometric surface measurements using inaccurate reflectance/illumination
knowledge,” Proc. CVPR, 2012.
[21] C. Hernández, G. Vogiatzis, and R. Cipolla, “Multi-view photometricstereo,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 30, no. 3, pp.
548–554, 2008.
[22] A. del Bue, J. Xavier, L. Agapito, and M. Paladini, “Bilinear modeling
via augmented lagrange multipliers (balm),” IEEE Trans. Pattern Anal.
Mach. Intell., vol. 34, no. 8, pp. 1496–1508, 2012.
[23] K. Tang, C. Tang, and T. Wong, “Dense photometric stereo using
tensorial belief propagation,” in Proc. CVPR, 2005.
[24] T. Wu and C. Tang, “Dense photometric stereo using a mirror sphereand graph cut,” Proc. CVPR, 2005.
[25] D. Miyazaki, K. Hara, and K. Ikeuchi, “Median photometric stereo
as applied to the segonko tumulus and museum objects,” International
Journal of Computer Vision, vol. 86, no. 2, pp. 229–242, 2010.
[26] K. Torrance and E. Sparrow, “Theory for off-specular reflection fromroughened surfaces,” Journal of the Optical Society of America, vol. 57,
no. 9, pp. 1105–1112, 1967.
[27] G. Ward, “Measuring and modeling anisotropic reflection,” Computer
Graphics, vol. 26, no. 2, pp. 265–272, 1992.
[28] S. Rusinkiewicz, “A new change of variables for efficient brdf representation,” in Proc. of Eurographics Symposium on Rendering, 1998.
[29] P. Tan, L. Quan, and T. Zickler, “The geometry of reflectance symmetries,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 33, no. 12, pp.
2506–2520, 2011.
[30] M. Chandraker and R. Ramamoorthi, “What an image reveals aboutmaterial reflectance,” Proc. ICCV, 2012.
[31] P. Tan, S. P. Mallick, L. Quan, D. Kriegman, and T. Zickler, “Proc.
CVPR,” Isotropy, Reciprocity and the Generalized Bas-Relief Ambiguity,
2007.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

Input

LS

R-PCA

Biquadratic (low)Biquadratic (all)16

PL-SBL (p = 1)

PL-SBL (p = 3)

Fig. 13. Experimental results with real datasets (D). We used three kind of datasets called Chocolate bear (25 images
with 261 × 421), Fat guy (40 images with 293 × 344) and Doraemon (40 images with 269 × 420). We show example of
input images, recovered surface normals and surface meshes (only for PL-SBL(p = 3)).

[32] N. Alldrin, “Toward reconstructing surfaces with arbitrary isotropic
reflectance: A stratified photometric stereo approach,” Proc. ICCV, 2007.
[33] T. Higo, Y. Matsushita, and K. Ikeuchi, “Consensus photometric stereo,”Proc. CVPR, 2010.
[34] I. Sato, T. Okabe, Q. Yu, and Y. Sato, “Shape reconstruction based on
similarity in radiance changes under varying illumination,” Proc. ICCV,2007.
[35] M. Chandraker, J. Bai, and R. Ramamoorthi, “A theory of differential
photometric stereo for unknown brdfs,” Proc. CVPR, 2011.
[36] P. Favaro and T. Papadhimitri, “A closed-form solution to uncalibrated
photometric stereo via diffuse maxima,” Proc. CVPR, 2012.
[37] W. M. Silver, Determining shape and reflectance using multiple images.
Master’s thesis, MIT, 1980.
[38] A. Hertzmann and S. M. Seitz, “Example-based photometric stereo:
shape reconstruction with general, varying brdfs,” IEEE Trans. PatternAnal. Mach. Intell., vol. 27, no. 8, pp. 1254–1264, 2005.
[39] S. Shafer, “Using color to seperate reflection components,” COLOR
research and application, vol. 10, no. 4, pp. 210–218, 1985.
[40] F. Romeiro and T. Zickler, “Blind reflectometry,” in Proc. ECCV, 2010.
[41] W. Matusik, H. Pfister, M. Brand, and L. McMillan, “A data-driven
reflectance model,” ACM Trans. on Graph., vol. 22, no. 3, pp. 759–769,
2003.
[42] F. Verbiest and L. V. Gool, “Photometric stereo with coherent outlier
handling and confidence estimation,” Proc. CVPR, 2008.
[43] M. Tipping, “Sparse bayesian learning and the relevance vector machine,” J. Mach. Learns. Res., vol. 1, pp. 211–244, 2001.
[44] D. Wipf and S. Nagarajan, “Iterative reweighted `1 and `2 methods
for finding sparse solutions,” IEEE Journal of Selected Topics in Signal
Processing, vol. 4, no. 2, pp. 317–329, 2010.
[45] S. Boyd and L. Vandenberghe, Convex Optimization.
Cambridge
University Press, 2004.
[46] E. Candés and T. Tao, “Decoding by linear programming,” IEEE
Transactions on Information Theory, vol. 51, no. 12, pp. 4203–4215,
2005.
[47] D. Wipf, B. D. Rao, and S. Nagarajan, “Latent variable bayesian models
for promoting sparsity,” IEEE Transactions on Information Theory,
vol. 57, no. 9, 2011.
[48] R. Cook and K. Torrance, “A reflectance model for computer graphics,”
ACM Trans. on Graph., vol. 15, no. 4, pp. 307–316, 1981.

[49] E. Lafortune, S.-C. Foo, K. Torrance, and D. Greenberg, “Non-linear
approximation of reflectance functions,” in Proc. ACM SIGGRAPH,
1997.
[50] A. Agrawal, R. Raskar, and R. Chellappa, “What is the range of surface
reconstructions from a gradient field ?” Proc. ECCV, 2006.

Satoshi Ikehata received the B.A. of Psychology in 2009 and the M.S. of Information Studies
in 2011 from the University of Tokyo respectively
and am working for Ph.D at the University ofTokyo supervised by Prof. Kiyoharu Aizawa. He
worked as an intern at the Visual Computing
Group at Microsoft Research Asia as a visiting
student in 2011. His main interests lie in Imagebased 3D Scene Reconstruction (photometricstereo, reflectance recovery, multi-view stereo),
Time-of-flight Imaging, Image Super-resolution,
Human Cognition and Perception.

David Wipf received the B.S. degree in electrical engineering from the University of Virginia,
Charlottesville, and the M.S. and Ph.D. degrees
in electrical and computer engineering from the
University of California, San Diego. He was later
an NIH Postdoctoral Fellow in the Biomagnetic
Imaging Lab at the University of California, San
Francisco where his research involved the development and analysis of Bayesian learning algorithms for functional brain imaging and sparse
coding. Currently he is with the Visual Computing Group at Microsoft Research Asia.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. *, NO. *, *** 2013

Yasuyuki Matsushita received his B.S., M.S.
and Ph.D. degrees in EECS from the Universityof Tokyo in 1998, 2000, and 2003, respectively.
He joined Microsoft Research Asia in April 2003.
He is a Senior Researcher in Visual Computing
Group. His areas of research are computer vision (photometric techniques, such as radiometric calibration, photometric stereo, shape-fromshading), computer graphics (image relighting,video analysis and synthesis). He is on the editorial board member of IEEE Transactions on
Pattern Analysis and Machine Intelligence (TPAMI), International Journal of Computer Vision (IJCV), IPSJ Journal of Computer Vision andApplications (CVA), The Visual Computer Journal, and Encyclopedia of
Computer Vision. He served/is serving as a Program Co-Chair of PSIVT
2010, 3DIMPVT 2011, and ACCV 2012 and a General Co-Chair for
ACCV 2014. He is appointed as a Guest Associate Professor at Osaka
University (April 2010-), Visiting Associate Professor at National Instituteof Informatics (April 2011-) and Tohoku University (April 2012 - ), Japan.
He is a senior member of IEEE.

Kiyoharu Aizawa received the B.E., the M.E.,
and the Dr.Eng. degrees in Electrical Engineering all from the University of Tokyo, in 1983,
1985, 1988, respectively. He is currently a Professor at the Department of Information andCommunication Engineering of the University ofTokyo. He was a Visiting Assistant Professor
at University of Illinois from 1990 to 1992. His
current research interest is in image processing
and multimedia applications, and he is currently
engaged in multimedia life log and 3D video etc.
He received the 1987 Young Engineer Award and the 1990, 1998 Best
Paper Awards, the 1991 Achievement Award, 1999 Electronics SocietyAward from IEICE Japan, and the 1998 Fujio Frontier Award, the 2002
and 2009 Best Paper Award from ITE Japan. He received the IBM Japan
Science Prize in 2002. He is currently on Editorial Board of ACM TOMCCAP, Journal of Visual Communications and Image Processing and
APSIPA Transactions on Signal and Information Processing. He served
as the Editor in Chief of Journal of ITE Japan, an Associate Editor ofIEEE Trans. Image Processing, IEEE Trans. CSVT and IEEE Trans.
Multimedia. He has served a number of international and domestic
conferences; he was a General co-Chair of ACM Multimedia 2012 andIEEE VCIP2012.

17