PREPRINT OF A REGULAR PAPER @PAMI2015

1

From Intensity Profile to Surface Normal:
Photometric Stereo for Unknown Light Sources
and Isotropic Reflectances
Feng Lu, Yasuyuki Matsushita, Imari Sato, Takahiro Okabe, and Yoichi Sato
Abstract—We propose an uncalibrated photometric stereo method that works with general and unknown isotropic reflectances. Our
method uses a pixel intensity profile, which is a sequence of radiance intensities recorded at a pixel under unknown varying directional
illumination. We show that for general isotropic materials and uniformly distributed light directions, the geodesic distance between
intensity profiles is linearly related to the angular difference of their corresponding surface normals, and that the intensity distribution of
the intensity profile reveals reflectance properties. Based on these observations, we develop two methods for surface normal estimation;
one for a general setting that uses only the recorded intensity profiles, the other for the case where a BRDF database is available while
the exact BRDF of the target scene is still unknown. Quantitative and qualitative evaluations are conducted using both synthetic andreal-world scenes, which show the state-of-the-art accuracy of smaller than 10◦ without using reference data and 5◦ with reference
data for all 100 materials in MERL database.
Index Terms—Uncalibrated photometric stereo, general reflectance, BRDF, intensity profileF

1

I NTRODUCTIONPHOTOMETRIC stereo recovers surface normals of ascene from a set of images recorded under varying
lighting conditions. A variety of methods have been
proposed in the literature under the assumptions ofLambertian reflectance and known lighting directions.
These assumptions, however, greatly limit the applicability of photometric stereo, because real-world scenes
exhibit more general reflectances and precise calibrationof lighting directions is laborious in practice. Therefore,
it is desired to develop a technique that works with general surface reflectance and unknown lighting directions
simultaneously with good accuracy.
In most photometric stereo approaches, observed intensities are treated as independent measurements rather
than ordered measurements. The ordered measurements
offer additional information that cannot be found in the
un-ordered measurements as demonstrated by clusteringof surface normals [1]. The strength of the ordered
measurements in photometric stereo has been shown bysome previous approaches [2], [3], [4], [5]. These previous methods make use of an intensity profile, which is asequence of radiance intensities recorded at a pixel undervarying lightings. These methods are shown effective• F. Lu and Y. Sato are with the Institute of Industrial Science, theUniversity of Tokyo, Tokyo 153-8505, Japan. E-mail: {lufeng,ysato}@iis.utokyo.ac.jp• Y. Matsushita is with Microsoft Research Asia, Beijing 100080, China.
E-mail: yasumat@microsoft.com• I. Sato is with National Institute of Informatics, Tokyo 101-0003, Japan.
E-mail: imarik@nii.ac.jp• T. Okabe is with Kyushu Institute of Technology, Fukuoka 820-8502,
Japan. E-mail: okabe@ai.kyutech.ac.jp
Manuscript received xxx.

under individual assumptions: calibrated light directionsand Lambertian reflectance assumption for [4], known
occluding boundaries [3], and calibrated reference object
with the same reflectance property [2].
In this paper, we further exploit the observations on
intensity profiles to develop two methods for accurate
surface normal estimation from a scene with generaland unknown isotropic reflectances observed under unknown lighting directions. Our methods are built upon
the notion of conditional linearity . In particular, we showthat there exist two strong relations: one between thegeodesic distance of intensity profiles and the angulardifference of surface normals, and the other between
the intensity distribution of intensity profiles and surface reflectance property. We investigate these relations
quantitatively and develop a technique that recovers
surface normals purely from intensity profiles (Sec. 3).
We also develop the second method in the case wherea database of bi-directional reflectance distribution functions (BRDFs) is available. For this setting, we present atechnique that recovers not only surface normals but also
unknown light sources with improved accuracy (Sec. 4).
Our methods make the following assumptions: (1) directional light sources of the same intensity are uniformly positioned around the target surface, (2) the test
surface is made of a uniform reflectance up to scaling.
We begin with these assumptions for development of
the method; however, later in this paper, we show thatthey can be relaxed to some extent. In addition tothese assumptions, we assume a common orthographic
camera and exclude extreme cases where the test surface
contains only a few planar regions. With these assumptions, our method is able to recover surface normals up
to a binary convex/concave ambiguity.

PREPRINT OF A REGULAR PAPER @PAMI2015

1.1

Previous work

A variety of previous studies have been conducted to relax the assumptions of the Lambertian reflectance modeland known lighting directions in photometric stereo.
Non-Lambertian reflectances have been handled either
by regarding non-Lambertian components as outliers, orusing more general reflectance models. The former classof methods finds non-Lambertian observations in a robust estimation framework using color-cues [8], medianfilter [9], rank-minimization [10], sparse regression [11],or Markov random field [12]. The latter class of methods
studies reflectance properties, such as bilateral symmetry [13], reflective symmetry about the halfway vector [14], isotropy and monotonicity [15], [16], and otherreflectance symmetries [17]. There are methods that use
complex parametric models, such as ones that use theWard model [18], [19], Lambertian+specular models [20]and bi-variate BRDF representation [21]. These methods
assume that the lighting directions are known.
There are uncalibrated photometric stereo techniques
assuming unknown lighting directions. Most of the existing methods attempt to resolve the Generalized BasRelief (GBR) ambiguity [7] with the Lambertian model. Various surface properties are used such as diffuse
maxima [22], specularity [23], [24], low-dimensional space [25], minimum entropy [26], interreflections [27],color profiles [28], reflectance symmetry [17], holisticreflectance symmetry [29], and certain configurationsof the light sources [30]. If considering the perspective
projection case, camera calibration can solve the ambiguity [31]. These methods rely on the assumption that the
diffuse reflectance component obeys the Lambertian law.
Therefore, the applicability to real-world surfaces thathave non-Lambertian diffuse components is restricted.
Handling both non-Lambertian reflectances and uncalibrated light sources simultaneously is far more challenging and has been less studied. Silver [32] andHertzmann and Seitz [2] capture reference objects withknown and proper reflectances together with the target

0.08

A

B

C

D

Pixel intensity

Our methods have the following advantages over theprevious approaches. First, our methods handle both uncalibrated light sources and general isotropic reflectances
simultaneously; therefore, both of the two proposed
methods remove the conventional assumptions about
known lighting directions and reflectances. Second, byusing a MERL BRDF database [6], our second method
is able to estimate lighting directions as well as material and surface normals, which has been difficult withprevious methods. Third, derived surface normals have
only a binary convex/concave ambiguity, while general
uncalibrated photometric stereo methods suffer fromgeneralized bas-relief (GBR) ambiguity [7]. Our average
normal estimation error evaluated for 100 real-world
materials is smaller than 10◦ without BRDF database (thefirst method) and smaller than 5◦ with BRDF database(the second method).

2

A&C
B
D

0.06
0.04
0.02
0
0

200

400

Light source position

600

Fig. 1. Illustration of intensity profiles. Surface points A,
B, and C have the same reflectance, but D is different. A,
C, and D have the same surface normal, while B has adifferent normal.
scene, while Ren et al. [33] captures a well-designed
calibration board with multiple reference materials andshapes. These methods use intensity values of reference
objects to determine surface normals of the test scene.
Because reference objects are captured with the test
scene, the effect of uncalibrated light sources is eliminated. Georghiades [34] uses the Torrance-Sparrow modeland proposes to optimize over a large set of variables.
Chandraker et al. [35] recover surface iso-contours fromimages taken under light sources positioned in a ring. An
initial normal needs to be assigned to determine surface
normals. Sato et al. [3] use the similarity of intensity
profiles to recover the normal difference with the assumptions about uniform light sources and Lambertianor Torrance-Sparrow reflectance models. Okabe et al. [36]use attached shadow codes to deal with general BRDFs
under dense and uniform light sources. Both of the two
methods require visible occluding contours for resolving
the ambiguity. Unlike these methods, our method does
not require any reference objects being captured under
the same lights with the test scene, and assumes noneof initial normals, certain reflectance models, visible
occluding contours, or strictly uniform light sources.

2 I NTENSITY P ROFILES FOR P HOTOMETRICS TEREO
An intensity profile is a sequence of radiance intensities
recorded at a pixel across varying lighting conditions. It
has been used in various contexts due to its useful properties when assuming distant lightings, no cast shadow,
and no inter-reflections.
Orientation-consistency: Intensity profiles become exactly the same, if and only if they correspond to thesame surface normal and material (A and C in Fig. 1).
Using this simple observation, surface normals can be
determined by looking up a pre-stored table indexedwith surface brightness values [32], or match the intensity profiles to those from a reference object that is madeof the same material [2].
Geometry-extrema: For isotropic BRDFs, intensity profiles reach the extremas synchronously, if and only ifthey correspond to the same surface normal (A and D in
Fig. 1) even for scene points with different materials. This

111
76
41
6
0.00

0.29

0.57

0.86

1.14

Euclidean distance

where I˜pl is the recorded intensity at the p-th pixel for ascene point (p = 1, . . . , P ), under the l-th lighting direction (l = 1, . . . , L), and Ipl is the normalized intensity.
Previous methods have shown that the similarity
of intensity profiles and surface normals are strongly
correlated [1], [3]. Similarity can be defined using the
Euclidean distance of two intensity profiles as kIp − Iq k2
in a straightforward manner. It indeed correlates withthe angular difference cos−1 (nT
p nq ) of surface normals
np and nq at scene points p and q; however, this linear
relationship only holds in a limited range as depicted
in Fig. 2 (left). To extend the range, our method uses ageodesic distance dG (Ip , Iq ) instead of kIp − Iq k2 as it
is used in [3] to measure the similarity of more distinct
normals (Fig. 2 (right)).
The geodesic distance corresponds to the distance ofthe shortest path between two nodes in a graph, whichis computed by adding up the Euclidean distances ofneighboring nodes along the path [37]. Here we refer
to “neighbors” in the space of intensity profiles rather

Angular difference [deg.] Angular difference [deg.]

180

Green-latex

145
111
76
41
6
0.00
0.00

0.76
2.81

1.53
5.61

2.29
8.42

3.05
11.22

Geodesic distance180

Delrin

145
111
76
41
6
0.00

1.58

3.17

4.75

6.33

Geodesic distanceAngular difference of Surface normals

Let us now assume evenly distributed but unknown
light directions and a scene with a uniform material; weshow later that these assumptions can be relaxed to someextent. Let np , nq ∈ R3×1 be a surface normal pair at pixel
locations p and q, and {Ip , Iq } be the corresponding pixel
intensity profiles in a normalized form as:
h
iT q
 1

L T
1
L
˜
˜
Σl (I˜pl )2 , (1)Ip = Ip , . . . , Ip
= Ip , . . . , Ip

Angular difference [deg.]

145

Estimation
error [deg.]

2.1

180

180
145
111
76
41
6
0.00

0.80

1.59

2.39

3.19

Geodesic distanceFig. 2. Using geodesic distance (right) preserves a linear
relationship over a greater range of angular differences incomparison with using Euclidean distance (left).
Angular difference [deg.] Angular difference [deg.]

observation is used for clustering surface orientations [1]without determining the exact orientations.
Similarity: Similarity of two intensity profiles is strongly
related to the angular difference between two surface
normals for the same material (A and B in Fig. 1).
Sato et al. [3] analyze this relationship and recoversurface normals in the cases of Lambertian and TorranceSparrow reflectances illuminated under evenly distributed light directions with an assumption of known occluding boundaries.
This paper makes a further observation about intensity profiles and introduces the property of conditional
linearity. Unlike [3], we do not restrict our analysis tocertain reflectance models. Instead, we take into accountmore general isotropic reflectances such as materials in
the MERL BRDF database [6].
Conditional linearity: For various real-world isotropic
reflectances, we observe a strong linear relation between
the distance among intensity profiles seen under evenly
distributed lightings and the angular difference of surface normals, up to a certain normal angular difference.
We also observe that the linear coefficient is materialdependent and closely related to the intensity distribution of the intensity profiles. These observations allow
us to develop uncalibrated photometric stereo methods
that work with unknown isotropic reflectances.

3

Angular difference [deg.]

PREPRINT OF A REGULAR PAPER @PAMI2015

Materials tested on the right180

Alum-bronze

145
111
76
41
6
0.00

3.77

7.53

11.30

15.07

Geodesic distance
180
145

Specular-orange-phenolic

111
76
41
6
0.00

1.27

2.55

3.82

5.10

Geodesic distance

10
5
0

20

50

80

Angular threshold δ [deg.]

Fig. 3. Examples of the linear relationship. Top: four
typical distributions and their partial linear fittings (solidlines) are shown under uniform lightings. Dotted lines
show non-uniform light cases. Bottom: setting different
angular threshold δ varies the estimation errors.

than pixel coordinates. In our case, we first compute
Euclidean distances between all possible combinationsof intensity profiles Ip and Iq as(
kIp − Iq k2 if kIp − Iq k2 < εpd(Ip , Iq ) =(2)+∞otherwise,
where each Ip determines a threshold εp ; it is empirically
set to the 10-th shortest distance from this Ip to allother Iq . Then geodesic distance between Iu and Iv iscomputed using the pre-computed Euclidean distances{d(Ip , Iq )} asdG (Iu , Iv ) = Dsp ({d(Ip , Iq )}),

(3)

where function Dsp (·) applies a Dijkstra’s algorithm tofind the shortest path that connects Iu and Iv , and adds
up multiple d(Ip , Iq ) along the path. In this manner, inthe geodesic distance, linearity is locally well preserved
and extended to a greater range of angular differences
as shown in Fig. 2 (right).

−1
cos−1 (nT(nT
u nv ) ≈ αm dG (Iu , Iv ) if cos
u nv ) < δ.

Normal A 0.3

0.2

0.2

0.1

0.1

200

Lighting positions

300

0.3

Normal A 0.3

0.2

0.2

0.1

0.1

0
0

0

100

200

Lighting positions

300

Material reflectance property

The linear coefficient αm described in the previous
section is material-dependent, i.e., it is related to thesurface reflectance property of a material. To characterizesuch a reflectance property, we begin with showing thatthe intensity distribution observed in an intensity profile conveys information about the reflectance property.
Fig. 4 shows four intensity profiles, where each row
corresponds to the same material but different surface
normals. These figures indicate that an intensity profile’s
shape depends on both reflectance and surface normal.
On the other hand, its intensity distribution appears
stable against surface normal changes while it varies
with materials as shown in Fig. 4.
Based upon this observation, we develop a new measure for characterizing the relationship between theintensity distribution and the corresponding material.
The intuitive examples in Fig. 4 suggest that a shinnymaterial shows a biased distribution due to sporadic
specular observations compared to a more symmetric

0.3

Normal B 0.3

0.2

0.2

0.1

0.1

0
0

0

Intensities

100

100

200

Lighting positions

300

0

0.3

Normal B 0.3

0.2

0.2

0.1

0.1

0
0

100

200

Lighting positions

300

0

Fig. 4. Intensity profiles w.r.t.material and surface normal.
The top row shows intensity profiles of a specular surface
recorded at two distinct surface normals. The bottom row
shows those of a diffuse surface. The intensity distributions are plotted on their right-hand side.

distribution from a diffused material. We therefore use
skewness for characterizing the relationship. The skewness γ of an intensity distribution, which is invariant to
the intensity/lighting order, is calculated as
 X
 32
X
1
(I l )3(I l )2
γ(I) = L 2
,
(5)
l

(4)

While the result is not sensitive to the choice of threshold δ, extreme values should be avoided. We examine
the error across varying δ in the bottom of Fig. 3. It
shows that for non-diffused materials, δ should be set
no larger than π/4 to guarantee a good estimation. As
for the material-dependent slope αm , the inference ofαm enables us to eliminate the necessity of the known
occluding boundary assumption that is used in [3] as we
will see later in this paper.
2.2

0.3

0
0

Intensities

We examine this linearity for all 100 materials in the
MERL BRDF database [6] by plotting the angular differences cos−1 (nTu nv ) and geodesic distances dG (Iu , Iv ).
Fig. 3 shows four typical plots, from which we can makethe following observations. First, the use of geodesic
distance generally shows the greater linear relationshipto the angular difference of normals in a larger range.
Second, for most materials, even those having complex
reflectances, the linearity still largely holds in the rangeof [0◦ , 70◦ ] of the angular difference, while it does notspan the entire range for many materials. Finally, thelinear coefficient, or slope, varies with the material (solidlines). The slope is insensitive to random noise in lighting directions as shown by the dotted lines in Fig. 3,
which are the line fittings to the plots produced withfluctuating light directions by 7◦ standard deviations.
A more thorough analysis on the linearity under nonuniform light distributions is given in Sec. 5.1.
Interestingly, the first case in Fig. 3 is actually the
problem solved in [3], which corresponds to an ideal
subset of the more general reflectances. Based on theabove observations, we define a partial linear conversion
from the geodesic distance dG (Iu , Iv ) to the angulardifference of surface normals cos−1 (nTu nv ) in a certain
range bounded by a threshold δ as

Intensities

4

Intensities

PREPRINT OF A REGULAR PAPER @PAMI2015l

where I is an intensity profile, and I l is its l-th elementthat corresponds to the l-th observation under the L
lighting direction. If we look at examples in Fig. 4 again,
the skewness value is large for the shiny material andis small for the matte one. Therefore, we quantitatively
examine the relationship between the skewness of the−1
intensity distributions and the linear coefficient αm.
To verify the effectiveness of the skewness measure,
we render a spherical surface under varying numbers oflight directions (162 and 42) using all the 100 materials in
MERL BRDF database. The skewness values γ computed−1from all pixels and inverse slopes αmare plotted inFig. 5. As we see in the figure, they have a strong
linear relationship, which in fact yields the correlation
coefficient of 0.99. Our experiment shows that varying
numbers of light directions show nearly identical linear
relations. From this result, we use the skewness of theintensity distribution for estimating αm of the unknown
material using the linear relationship in Fig. 5, and it
results in constraining Eq. (4).
2.3

Photometric stereo using intensity profiles

The two linear relationships described in the previous
sections 2.1 and 2.2 are our key findings in this paper.
Based on them, we develop two uncalibrated photometric stereo methods that work for unknown isotropic
BRDFs. We begin with the case where only intensity
profiles are available in Sec. 3, and extend the method
for the case where a database of BRDFs is available
in Sec. 4. We denote the former method as a method
without reference data, and the latter as a method with

PREPRINT OF A REGULAR PAPER @PAMI2015

510

N = [n1 | n2 | . . . | nP ] ∈ R3×P that we solve for. We define the observation matrix O as

m

α-1

8
6

P ×P
O = [ou,v = nT,
u nv ]P ×P ∈ R

4
2
1

2

3

4

5

6

skewness

7

8

5
m

α-1

4
3
2
1
1

1.5

2

2.5

3

skewness

3.5

4

−1
Fig. 5. Relationship between skewness values and αmvalues of 100 materials. The correlation coefficients are
0.99 (Top/bottom: 162/42 light. directions)reference data. The method without reference data uses
the two linear relationships to quantitatively convert theobserved intensity profiles into surface normals, while
the method with reference data casts the original problem into a matching of intensity profiles, where its search
space is bounded by the two linear relationships.
In practice, they can be applied in accordance withthe scenarios: quick implementation without a pre-stored
database, or simultaneous light source recovery and agenerally better guaranteed accuracy (Sec. 5).
In what follows, we begin with describing the method
without reference data to explain how surface normals
can be determined only from intensity profiles. Then we
describe the method with reference data, in which the
problem can be more constrained.

3

S URFACE NORMAL ESTIMATION WITHOUTREFERENCE DATA

We describe our first method for recovering surface
normals only from the intensity profiles without reference data. From the observed intensity profiles, we
compute geodesic distance dG (Iu , Iv ) using Eqs. (2) and(3). We then convert dG (Iu , Iv ) into a normal angulardifference cos−1 (nTu nv ) using Eq. (4). As discussed in
the previous section, since the conversion is valid only
when cos−1 (nTu nv ) < δ holds, we rewrite Eq. (4) as thefollowing:(
αm dG (Iu , Iv ) if αm dG (Iu , Iv ) ≤ δ−1Tcos (nu nv ) =Undefined
otherwise,(6)
where αm is obtained using the skewness measure, andthreshold δ is set to π/4 to ensure good linear regions
for diverse classes of materials as discussed in Sec. 2.1.
3.1

Formulation

We wish to recover surface normals of scene points
that correspond to P pixels in the observed image from a set of images taken under varying unknown lightings. Let the surface normal matrix be(7)whose elements ou,v (= nT
u nv ) are readily obtained from
Eq. (6). The diagonal elements ou,u in O are all ones,
which ensures the unit normal length for {np }. Equation (6) has undefined cases for some u and v; therefore,
not all O’s elements have well-defined values. In otherwords, the observation matrix O has missing elements.
With a matrix E that accounts for the errors due to themissing entries, the relationship between the observation
matrix O and surface normal N can be written as
O = N T N + E.

(8)

We wish to solve for surface normal N by using theincomplete observation matrix O and unknown but
sparse error matrix E.
3.2

Matrix decomposition with missing data

Solving Eq. (8) for N involves recovering and decomposing the incomplete observation matrix O. We use
a matrix A for A = N T N , where we know thatrank(A) = 3 since rank(N ) = 3. Let Ω be a set ofindices where ou,v has a defined value in O, and let its
complement set be Ωc . By restricting the error matrix E
only to account for the missing entries Ωc , the original
problem of Eq. (8) can be written as
argmin kA − O + Ek2F

s.t. rank(A) = 3, ΦΩ (E) = 0,A(9)where ΦΩ (E) is an operator that only keeps E’s entriesin Ω unchanged and sets others zero. We solve the problem of Eq. (9) by alternately estimating A and E. The
optimization begins by initializing E = 0 and setting themissing entries of O to zeros so that ΦΩc (O) = 0.
At the k-th iteration, we update Ak+1 by singularvalue decomposition (SVD) as
T U SV ← SVD(O − E k ),
 

(10)
 Ak+1 ← U S I (3) 0 V T = U S (3) V T ,00where Ak+1 is reconstructed using only the first three
singular values of S to ensure rank(Ak+1 ) = 3. I (3) is a
3 × 3 identity matrix. We then update E k+1 by
E k+1 ← ΦΩc (O − Ak+1 ).

(11)Equation (11) assigns values to E k+1 only for thoseentries in Ωc . This ensures the second constraint of Eq. (9)to hold.
The iteration stops when it converges using the criterion kE k − E k+1 kF < ξkOkF , where ξ is a small value (setto 10−4 ). Finally, the method obtains the surface normal
estimates N̂ as
1

1

2
2
U T = S (3)V T.
N̂ = S (3)(12)

PREPRINT OF A REGULAR PAPER @PAMI2015

3.3

6

Like any other uncalibrated approach, our solution contains ambiguity. In our case, any matrix Q ∈ R3×3 thatsatisfies QT Q = I can be multiplied with the solutionof Eq. (8) without violating the equality [38]:N T QT QN = (QN )T (QN ) = O − E.

Reference data

We produce a reference dataset containing intensity values that are sampled from a large variety of surface
normals, light directions, and materials. Let srk ∈ R3
(k = 1, . . . , S) and nri (i = 1, . . . , N ) denote the sampled light directions and surface normals, respectively,
and let m(= 1, . . . , M ) be the index for materials included in the reference dataset. Then I r (k, i, m) is them

3D reference data array Dr
Determinematerial m
ˆ

……

stl

2D test data matrix Dt

Iteration

Intensity
matching

Row/column
operations

……

nri

……

s

r
k

r
2D slice Dmˆ…

Section 3 recovers surface normals for general isotropic materials using only intensity profiles without any
reference data. In this section, we describe an extended
method by assuming that a database of BRDFs is additionally available.
A BRDF database includes dense measurements of realmaterial reflectances, and diverse appearances can besynthesized using varying surface normals and lightings.
If such a synthetic data is dense enough, any test data
can be considered as its subset, and the problem can bemade more bounded and tractable by using the known
synthetic data as a reference. In this extended method,
we first generate the reference dataset using the BRDF
database, and then use the reference data for deriving
the estimates of surface normals and lighting directions.

ntj

……

skr

stl

…

N ORMAL ESTIMATION BY USING REFER ENCE DATA

ntj

nri

…

4

m

skr

…

T

Satisfy Q̃ Q̃=I

This is a binary ambiguity where λ = ±1 corresponds
to convex/concave surfaces that are not distinguishable
without light calibration. Therefore, by using the integrability constraint, our method recovers surface normals up to only a binary convex/concave ambiguity. To
apply the integrability constraint, our method uses theapproach of [39].

4.1

nri

(13)

Therefore, QN ∈ R3×P is also a solution to the problem.
In a previous study, Belhumeur et al. [7] show that if
rank(N ) = 3, the ambiguity due to a general invertible
3 × 3 matrix can be reduced to the generalized bas-relief(GBR) ambiguity using the integrability constraint [39].
This is also true in our case since Q ∈ R3×3 . Therefore,
by enforcing integrability, we have a simpler ambiguity
described by a matrix Q̃. As a result, in our case, Q̃should take a form of a GBR transformation, and at theTsame time, satisfy Q̃ Q̃ = I. Therefore, it becomes




1 0 µ
1 0 0
1
1
Q̃ ⇒ 0 1 ν  ⇒ 0 1 0 , λ = ±1. (14)
λ
λ
0 0 λ
0 0 λ
|
{z
}
|
{z
}
GBR transformation

Captured test data

Synthesized reference data

Convex/Concave ambiguityr

Reformed 2D slice RDm
ˆCFig. 6. Matching-based photometric stereo: determining
surface normals, light directions, and material of the test
scene by matching intensity values between test data andreference data.

corresponding pixel intensity value for k-th light, ith surface normal, and m-th material. In other words,
the reference dataset {I r (k, i, m)} defines a mapping(srk , nri , BRDF(m)) 7→ I r (k, i, m). In the rest of the paper,
we use superscripts r and t for indicating reference andtest data, respectively.
We assume that all light sources have an identicalintensity; therefore, both light source srk and surface
normal nri are represented as unit vectors. To densely
sample the reference dataset I r , we use 1249 frontal
vertices on a uniformly subdivided icosahedron to define
light directions srk and surface normals nri visible to thecamera. Pixel intensity values are then rendered for eachcombination of light source srk , surface normal nri , andmaterial m selected from the MERL BRDF dataset [6]that contains 100 isotropic reflectances. In this manner,
we obtain the reference data {I r (k, i, m)} for all combinations of sampled light directions, surface normals, and100 different surface reflectances.
4.2

Formulation

We formulate the problem as intensity matching between
test data and reference dataset. For a test scene, we
obtain its test data by recording a sequence of images
from a fixed viewpoint with varying directional illumination. Let stl (l = 1, . . . , S t ) denote the unknown lightvectors and ntj (j = 1, . . . , N t ) be the unknown surface
normals for the recorded pixels. Their corresponding
pixel intensity values {I t (l, j)} are stored as the test data,
which form a test matrix D t = [I t (l, j)]S t ×N t .

PREPRINT OF A REGULAR PAPER @PAMI2015

7

Let D r = [I r (k, i, m)]S r ×N r ×M r denote the three
dimensional tensor containing all the intensity values
from the reference dataset. Then, our task is to match
intensity values between D r and D t . To reduce thecomputational complexity, our method first obtains theestimate of material m̂ using the skewness measuredescribed in Sec. 2.2. Then, its corresponding 2D matrix
D rm̂ = [I r (k, i, m = m̂)]S r ×N r , which is a slice of D r , isused for matching with D t to solve for surface normals
and light directions. An overview of this approach isillustrated in Fig. 6.
To match pixel intensities between D rm̂ and D t , lettr
us define two operators as two matrices R ∈ RS ×Srt
and C ∈ RN ×N . When the reference matrix is leftmultiplied by R, its q-th row is selected and placed in thep-th row of the resulting matrix. Similar to a permutation
matrix, this is achieved by setting Rp,q = 1 to be the only
non-zero element in the p-th row of R. Similarly, matrix
C applies column selection. When the reference matrix
is right-multiplied by C, its p-th column is selected andplaced in the q-th column of the resulting matrix. Thisis performed by setting the element Cp,q = 1 to be theonly non-zero element in the q-th column of C.
Using these operators R and C applied to the reference matrix D r , we can use pixel intensities of theselected rows (light directions) and columns (surface
normals) of the reference matrix D r to compare withthe test matrix D t . If m, R, and C are properly chosen,
the following objective function can be minimized:{R̂, Ĉ, m̂} = argmin kRD rm C − αD t kF ,(15)

R,C,m,αwhere α is a scalar value to account for exposure
differences. This matching procedure is illustrated inFig. 6 (right). In the optimal solution set {R̂, Ĉ, m̂}, m̂ isconsidered the index of the reference material that best
approximates the unknown test material mt . Besides,
positions of ‘1’s in R̂ relate each unknown light source
directions stl to their closest reference light source directions in {srk }. Similarly, Ĉ maps each unknown surface
normals ntj to their closest reference surface normals in{nri }. In this manner, the problem can be formulated
as finding the best matching between the test data andreference data.
Due to the discontinuous nature of the operators R
and C, where only binary values {0, 1} are contained,
directly finding the global optima via optimization iscomputationally difficult. We solve this problem in thesubsequent sections in an iterative refinement framework. Namely, our strategy is to first determine the material m̂, then infer light directions, and finally estimate
surface normals, and iterate these steps for reaching toa better solution.
4.3

Material estimation

To solve the problem of Eq. (15), we first obtain theinitial solution by inferring the material. By assuming

a uniform surface reflectance in the test scene, or asurface patch with a unique material after segmentation,
we infer the surface reflectance from the observed pixel
intensity values by computing the skewness measureusing Eq. (5). Let the skewness for the test materialr
and those of 100 reference materials be γ t and {γm},respectively. Then a best matching material m̂ is foundas:rm̂ = argmin kγ t − γm
k2 .
(16)
m

However, for highly specular materials, if the test surface
is lit by much fewer light directions than those in thereference data, γ t computed for a single pixel may beunstable due to the limited angular resolution of lightings. To avoid such a situation, we take multiple pixels
with similar surface normals (whose intensity profiles
have positive correlation coefficients), and use all ofthem to compute γ t . We also retain a list (empirically 20)
matching candidates {m̂} with their ascending matchingcosts, among which the final solution will be obtained
via optimization described later in this paper.
4.4

Initialization of light directions

In this section, we describe a method for initializing light
directions, i.e., corresponding to the estimation of R, by
assuming a roughly uniform light direction distribution.
In estimating the light directions, our method takes a
two-step approach, i.e., it first determines the azimuth
angles, then estimates the elevation angles.
To determine the azimuth angles of the light directionsin the test data, we use the angular difference calculation
described in Sec. 2.1. We begin with approximately determining two special elements in the test data: the light
direction st and surface normal nt , which are parallelto the viewing direction. To identify st , we take a simple
strategy by selecting the lighting st that illuminates thegreatest number of surface points. Specifically, we examine every image captured under a different light stl , andselect the one that contains the largest number of pixels
whose intensity values are greater than a threshold; thethreshold for detecting shadows is determined based on
the first quartile of all measured intensities, following[36]. The frontal surface normal nt is then determined by
selecting the surface normal that produces the brightest
pixel intensity under the frontal lighting direction st .
At the same time, we find one set of surface normals{nte } with similar elevation angles as illustrated in Fig. 7
after grouping similar intensity values under the frontal
light source sr . From this set, our method estimates theazimuth angles of {nte } using their intensity profiles {Ite };specifically, the linear relationship between the angulardifference of surface normals and geodesic distance of
their intensity profiles described in Sec. 2.1. For obtaining
angular differences among {nte }, we apply Multi Dimensional Scaling (MDS) [40], which is a common techniquefor distance-preserving dimensionality reduction using
the geodesic distances among {Ite } as input. The output

PREPRINT OF A REGULAR PAPER @PAMI2015

Viewing
direction

8

s//t
n //t
{net }

Same elevation
angle

Fig. 7. For light source estimation, we find light source stand surface normal nt parallel to the viewing direction,
also a set of normals {nte } with similar elevation angles.

while notation {stl } does not change for convenience. By
assuming a uniform distribution of {stl } in the frontal
hemisphere, we order uniform reference lights {srk } inaccordance with their elevation angles, sub-sample theirelevation angles to a total number equal to the numberof {stl } and assign them to {stl }.
In this procedure, we have assumed a uniform lightsource distribution for the purpose of initialization.
However, iteration in the next section removes this assumption and handles non-uniform lights for deriving
the final solution. The 2D rotation and flip ambiguity inthis section will also be addressed later.
4.5

Rotationambiguity

Colorized
azimuthanglesRotation & flip
ambiguity

Estimated azimuth angles via dimension reduction

Surface normal estimation

Previous sections estimate a list of candidate materialsm̂ and light selection matrix R̂. What remains unknown
in Eq. (15) is Ĉ that encodes the surface normals of thetest scene.
Let A = R̂D rm̂ , then Eq. (15) can be rewritten asĈ = argmin kAC − αD t kF ,

(17)

C,α
Rotationambiguity

Rotation & flip
ambiguity

Fig. 8. The estimated azimuth angles of {nte } for twomaterials are shown in two rows. The spots’ positions in
the plots are the dimensionality reduction results of {Ite }.
They correctly convey the relative azimuth angles of {nte }and the color shows the ground truth azimuth angles.
Note the 2D rotation and flip ambiguity.
are a set of 2D coordinates that preserve the distancesamong {Ite }. The geodesic distances among these 2D
coordinates are linearly related to the angular differenceof {nte }; in other words, they show a similar structurewith that of {nte } in the 2D space up to a 2D rotation andflip ambiguity as shown in the example in Fig. 8. Since{nte } have similar elevation angles, the 2D coordinates
represent the distribution of their azimuth angles.
For each light source stl , we examine its corresponding pixel intensity values at surface normals {nte }. The
brightest intensity observation should correspond to oneof {nte } that has the closest azimuth angle to stl assumingthat the reflectance is isotropic. Based on this, we assign
an azimuth angle obtained from {nte } to each light source{stl }. Note that, at this point, the 2D rotation and flip
ambiguity in {nte } is still embedded in stl .
Once the azimuth angles of {stl } are determined, theirelevation angles are ordered in accordance with thebrightness of the intensity observations at the frontal
surface normal nt , based on the fact that when the light
direction’s elevation angle gets closer to the frontal surface normal nt , it yields a brighter intensity observation.
In the meantime, we remove light sources from {stl }which cannot illuminate nt to only keep frontal lights,

where C contains a single ‘1’ in each column while
all the other elements are zero. Because of this discontinuous and non-differentiable constraint, Eq. (17) is
not computationally tractable by efficient convex/nonconvex optimization methods. Instead of directly optimizing Eq. (17), we use a nearest-neighbor searchin the normalized column vector domain. Namely, wesequentially select a column of D t and find the mostsimilar column in A. The similarity is measured using
the normalized column vectors to eliminate the effect ofα, and then computing the cosine of their inner product.
For instance, for the q-th column vector of D t , if its mostsimilar column in A is found to be the p-th one, Cp,q = 1
is set to be the only non-zero element in C’s q-th column.
By sequentially determining each column of C, correspondences between test surface normals {ntj } andreference normals {nri } are obtained to give us an initial
estimate of the test surface geometry. At this point, a 2D
rotation and flip ambiguity of light directions (Sec. 4.4)
are still embedded in the recovered surface normal {ntj }.
4.6

Iterative refinement

To further refine the estimates of light and surface normal selection matrices R̂ and Ĉ, we take an alternating
update procedure. For updating R̂, we fix the previous
estimate of Ĉ and compute α̂ using Eq. (17). Then R̂ is
updated byR̂ ← argmin kRD rm̂ Ĉ − α̂D t kF ,

(18)

R

which is solved in a similar manner as Eq. (17). We
then fix R̂ and update Ĉ using Eq. (17) as describedin Sec. 4.5. This procedure continues iteratively until thecost decrement of Eq. (17) becomes smaller than 0.1%.
We also define an outer loop aiming at finding the best
matching material. As described in Sec. 4.3, we have a

PREPRINT OF A REGULAR PAPER @PAMI2015

9

Algorithm 1 Reference-based method algorithm
Find candidate materials {m̂} via Eq. (16).
// Outer loop: test each candidate material
for Select a new candidate material m̂ do
Initial estimation of light directions R̂ (Sec. 4.4)
Initial estimation of surface normals Ĉ (Sec. 4.5)// Inner loop: update light directions and surface normals
while Not converge do
Update R̂ by using Eq. (18)

Fig. 9. Synthetic surfaces used in the experiments: hemisphere, Bunny, Dragon, Rabbit, Beethoven, Lion, andHappy Buddha.

Update Ĉ by using Eq. (17)Update the minimum cost value.
end while

TABLE 1
Average normal / light errors of 100 materials. Withreference data, only frontal lights are used (89 and 25).

Update the solution with minimum cost value
end for
Return the final solution {m̂, R̂, Ĉ}list of candidate materials obtained by evaluating theskewness. For each of them, we estimate the light andsurface normal selection matrices R̂ and Ĉ using theprocedure described above. Among those material-wise
estimates, the minimizer {R̂, Ĉ, m̂} of the cost of Eq. (15)is taken as the final solution. The entire procedure issummarized in Algorithm 1.
4.7

Convex/Concave ambiguity

As discussed in the end of Sec. 4.5, the solution is not yetcorrectly determined since it contains ambiguity in thesurface normal estimates. Namely, the derived surface
normals can rotate and flip about the viewing direction
together with light sources by still satisfying Eq. (15).
Similar to Sec. 3.3, the integrability constraint can beused to reduce this ambiguity to a binary one. The 2D
rotation and flipping ambiguity B for a surface normal
n is written as± cos φ − sin φ 0cos φ 0
Bn = 
(19)
 ± sin φ
 n,
0
0
1
where φ is the 2D rotation angle around the viewing
direction. Using the integrability constraint, the general
invertible 3 × 3 ambiguity can be reduced to a GBR ambiguity. Therefore, in our case, the GBR transformationthat satisfies our particular ambiguity B becomes




κ 0 0



B̂ = 
 0 κ 0 ,0 0 1

κ = ±1,(20)

which correspond to two surfaces whose surface normals
differ by a π rotation. Therefore, the ambiguity corresponds to a binary convex/concave ambiguity. Applying
integrability constraint in our case is simple, mainly because the original ambiguity in Eq. (19) only contains anangle φ and sign ‘±’ as unknowns. Compared to general
invertible 3×3 matrices solved by conventional methods,
our problem is much easier to solve numerically.

Mat–erialsAll 100
Best 67
Best 33

162 (89? ) lights
No ref. With ref.?
9.65◦ /– 4.36◦ /5.45◦
6.77◦ /– 2.58◦ /2.98◦
3.65◦ /– 1.94◦ /2.03◦

42 (25? ) lights
No ref. With ref.?
9.30◦ /– 4.88◦ /5.84◦
7.02◦ /– 3.13◦ /3.64◦
4.72◦ /– 2.18◦ /2.32◦

Since our iterative update scheme does not explicitly
enforce the ambiguity to be a rotation and flipping
ambiguity, for a Lambertian surface, the ambiguity may
become the GBR ambiguity during the iterative update.
This is due to the fact that Eq. (15) has minimizers up to
the GBR ambiguity in the Lambertian case. For general
non-Lambertian reflectances, the ambiguity will not be
generalized as such because of the information provided
by specular components. This is consistent with previous
studies that resolve the GBR ambiguity using specularobservations [23], [24], [34].

5

E XPERIMENTAL RESULTS

In this section, we evaluate our methods, one without
reference data (Sec. 3) and the other without reference
data (Sec. 4), using synthetic and real-world data for bothquantitative and qualitative assessments.
5.1

Synthetic scenes

We quantitatively assess the proposed methods usingsynthetic scenes. We use several 3D models include
a hemisphere, Bunny, Dragon, and Happy Buddha
from Stanford 3D model dataset [41], and also Rabbit,
Beethoven and Lion as shown in Fig. 9. These 3D models
are rendered using BRDFs of the MERL BRDF database
under varying lighting directions. The rendering settings
are detailed in the each test described below.
Estimation accuracy We first assess the estimation
accuracy of surface normal and light directions usinga hemisphere that contains all visible surface normals.
We render images using 100 different materials of the
MERL BRDF database illuminated under multiple lightPREPRINT OF A REGULAR PAPER @PAMI2015

Average error [deg.]

10
8

10

Angularerror

No ref.
With ref.

40

6

30

4

20

2
0

20

40

60

Normal angle range

80

10
0

TABLE 2Comparison of average angular errors of surface normal
estimates for 100 materials and 162/42 light directions.
Ours
Ours?
Wu et al. Hayakawa(no ref.) (with ref.)
[10]
[38]
162 (89? ) 9.60◦
4.36◦
18.52◦
22.09◦
42 (25? )
9.30◦
4.88◦
20.05◦
26.95◦
Lightsdirections uniformly distributed over a sphere. In particular, we first test 162 uniform light directions, and then
42. For the method with reference data (Sec. 4), we only
use the frontal light directions that span in the upper
hemisphere, resulting in 89 and 25 light directions.
Table 1 shows the average accuracy for 100 materials.
Without reference data, our method in Sec. 3 achieves
average errors smaller than 10◦ for all 100 materials,
while for the best 33 materials the average accuracy isaround 4◦ . On the other hand, if using reference data,
our method in Sec. 4 has better accuracies than 5◦ ,while the number for the best 33 materials is around 2◦ .
In addition, it also allows for light sources estimation,
whose accuracy is similar to those for surface normals.
We study the distribution of surface normal errors.
Fig. 10 shows the relationship between normal angle
range and average error. Without reference data, ourmethod in Sec. 3 accurately recovers surface normals
when they are more close to the viewing direction;
while errors become bigger for those normals in theboundary. On the other hand, our method in Sec. 4
using reference data performs more stably for variant
surface normals. Note that when using reference data,
the discrete sampling of srk and nri in Sec. 4.1 is one
source of the error. In particular, the average samplinginterval is about 4.1◦ in our setting, and thus it will cause
a 1.6◦ average error for arbitrary surface normal or lightsource. This error exists in the final results inevitably.
Comparison with previous methods In this section,
we further compare our results with those of previous
methods. The same dataset as above, i.e., hemisphericalsurfaces rendered with 100 materials and both 162 and42 light directions, are used by all methods. Strictly
speaking, it is not easy to find prior methods that cancompletely handle unknown reflectances and uncalibrated illuminations without additional assumptions to ours.
In our experiments, we choose the following ones thatat least separate the reflectance and illumination factors
without requiring known light directions:1) Hayakawa [38]: a baseline method that assumes
Lambertian reflectance and no shadows.
2) Wu et al. [10]: a state-of-the-art method that robustly handles non-Lambertian components and also
cast shadows.
During experiments, we find these methods difficult tosolve the GBR ambiguity for very specular materials.
Therefore we provide them with the ground truth light

Surface normal error [deg.]

Fig. 10. Average surface normal error w.r.t.normal range.
Spherical plots show the average error maps.

70
60
50
40
30
20
10
0

Ours (no ref.)Ours (with ref.)Wu et al. 2010
Hayakawa 1980

1

8 15 22 29 36 43 50 57 64 71 78 85 92 99

Materials (sorted by the results of Wu et al. 2010)

Fig. 11. Comparison of average normal errors for 100
materials with 42 lights. Large errors are truncated at 70◦ .

sources. In this sense, they actually work in a calibrated
mode and thus give their theoretically best results.
Average estimation errors of all methods are shownand compared in Table 2. Our method achieves significantly better estimation accuracy than others for both
162 and 42 light sources. Note that errors for methods[42] and [10] are around 20◦ in average, but it does notmean these methods cannot handle all materials. In fact,as per-material results shown in Fig. 11 for 42 lights, ourmethod clearly achieves good accuracies for all materials
while the other two methods only fail for those difficult
ones. This shows the significant merit of our method tohandle general and unknown reflectances.
Another significant merit of our method is that, if reference data is available, we can also recover light source
directions with high precision. This cannot be done byconventional methods without accurately knowing bothreflectance and shape.
Inexact reference material We test the method witha reference database where the exact reference material
does not exist. We perform the experiments for 100 materials by forcing the method not to choose the exact reference material. The resulting accuracy does not reducesignificantly, as shown in Table 3. It shows that mostcommon reflectances can be approximated by similar
ones in the database when only considering their grayscale values. In addition, we also try interpolating the
unknown material by using three best matched inexact
materials during optimization, where the weights are
set inversely proportional to their individual matchingcosts. The accuracy in Table 3 shows an even more stable
result with the interpolation method than when given
the exact reference material, since the non-interpolation
method may not always select that best material.
Convergence The problem in Algorithm 1 is clearlynon-convex, and we test its convergence experimentally.

PREPRINT OF A REGULAR PAPER @PAMI2015
Avg. 5°shift
full 42 lights

Avg. 13°shift
full 42 lights

Global skewing
full 42 lights

TABLE 4
Normal error increments due to non-uniformity of lights.

Lightsources

Frontal25 lights

11

5.8°

5.8°

6.8°

10.0°

Surface normal error [deg.]

Estimated Estimated lightnormal map
sources

Uniform Frontal 5◦ shift 13◦ shift Skewing
No ref.
9.30◦+1.67◦ +0.67◦ +1.74◦+6.17◦◦◦◦◦
W. ref.
4.88
+0
+0.51
+1.48
+3.70◦

Normalerrors

40
20
0

5.4°

4.9°

6.4°

8.6°

Recovery with reference data
Frontal25 lights

Avg. 5°shift
full 42 lights

Avg. 13°shift
full 42 lights

Global skewing
full 42 lights

Frontal, 25 lights

35

Avg. 5°shift, 24 lights

30

Avg. 13°shift, 24 lights

25

Global skewing, 21 lights

20
15
10
5
0
1

7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97

Materials (sorted by the results of uniform lights)

Fig. 14. Per-material surface normal estimation errors
under different non-uniform light sources.

40

Normalerrors

40

20
0

10.0°

11.0°

11.0°

Recovery without reference data

15.5°

Fig. 12. Average surface normal estimation results for
100 materials under non-uniform distributions of light directions. Values represent the average angular errors.
TABLE 3Results with/without exact reference material.
Lights Exact mat. Inexact mat. Interpolated mat.
89
4.36◦
4.65◦
4.25◦
◦
◦
25
4.88
5.31
4.86◦

Normal error [deg.]

We select eight materials which represent four typical reflectance types: exceptional, specular, glossy and diffuse.
We test the convergence with different initial states byadding Gaussian noise to the ground truth light directions to cause average biases from 0◦ to 50◦ . Results in
Fig. 13 show that: 1) convergence becomes bad for glossy
materials with sharp specular lobes when the initial lightbias > 20◦ ; 2) some exceptional materials have globallylarge errors since they are not ideally isotropic; 3) othertypes of materials show good convergence even withlarge initialization biases. For reference, the average biasof our light direction initialization is below 10◦ .
Non-uniform distribution of light directions To as-

30
20
10
0
0

10

20

30

40

Average initial light directions bias [deg.]50

Fig. 13. Results by distorting the initial light directions.
Each material on the left corresponds to an error curve.

sess the effects of non-uniform distributions of light
directions, we test our method with four representative
non-uniform light patterns. These light source directions
are generated by 1) only keeping the frontal ones ofthe original uniform lights; 2) randomly shifting theuniform lights with an average bias of 5◦ ; 3) randomly
shifting the uniform lights with an average bias of 13◦ ;and 4) globally skewing the uniform lights. These lightpatterns and the estimation results are summarized inFig. 12. Note that with reference data, our method inSec. 4 by default only uses frontal lights.
Intuitively, both surface normal map and light source
pattern can be correctly obtained. Randomly shifted lightsources do not affect the accuracy too much, while global
skewing causes relatively larger errors. This intuition
can be supported by the error increments given in
Table 4, where the errors increase by 10% ∼ 20% due
to random shifts, and the number is around 70% in
the cases of global skewing. These results indicate thatuniform light sources are preferred by our method, andnon-uniformity due to slight random shifts is also well
acceptable. Good accuracy is achievable by just assigning
light sources to be roughly uniform in practice.
Per-material results are further given in Fig. 14 for easycomparison. The results obtained by using reference data
are shown. The uniform-light case achieves best accuracies, while other light sources produce bigger errors but
maintain the same error-material correspondences.
Results for synthetic 3D surfaces We conduct experiments using synthetic 3D surfaces shown in Fig. 9. They
are rendered with six typical non-Lambertian materials
under 42 light sources with 5◦ random shifts. Their
normal maps are recovered by using our methods bothwith/without reference data. Representative results are
given in Fig. 15. The recovered normal maps by two
methods look highly similar, although error maps showthat using reference data improves accuracy especiallyfor boundary areas. Surface normal errors and light

12

40

0

Angular error

Estimated normal
map with ref. data

Angular error

Estimated normal
map w/o ref. data

One of input
imagesPREPRINT OF A REGULAR PAPER @PAMI2015

40

0

Fig. 15. Results for synthetic surfaces with isotropic reflectances. Top row shows one of the images for each scene.
Second and third rows respectively show the recovered surface normal and its error map of our method without
reference data (Sec. 3). Fourth and fifth rows show those of our method with reference data (Sec. 4).
TABLE 5Results for 3D models with six non-Lambertian materials.
No ref.
With ref. With ref.
3D Modelnormal err. normal err. lights err.
Bunny
8.91◦
4.41◦
6.00◦
◦
◦
Dragon
8.06
3.72
5.12◦
◦
◦
Rabbit
8.30
4.39
5.84◦
◦
◦
6.65
4.46
5.96◦
Beethoven
◦
◦
Lion
7.12
3.72
5.33◦
◦
◦
Happy Buddha
6.23
4.22
5.44◦
◦
◦
Average
7.55
4.15
5.62◦

source errors are further provided in Table 5. The average normal errors are 7.55◦ and 4.15◦ for all six nonLambertian reflectances with uncalibrated light sources.
5.2

Real-world scenes

We conduct experiments for real world scenes. Images
for real-world objects are captured by rotating the lightsource around the objects. Some objects with large concavity suffer from severe cast shadows produced bysidelights and back-lights. To avoid this, we use lights
mainly from the frontal hemisphere for these objects. In
most cases around 50 light sources are finally used.

Experimental results are summarized in Fig. 16 where
estimated surface normal maps are presented in color.
Results in the second row are obtained without usingreference data, and results in the last row are obtained byusing reference data. Both methods give nearly identical
outputs. Note that different colors on the same object are
well handled, because our methods use gray-scale pixel
intensities that are normalized. Besides, in all results,
there exists a convex/concave ambiguity as discussed.
The convex/concave surfaces are theoretically undistinguishable by our method. However, if use additional
information such as simple user interaction or boundary
assumption, it is feasible to resolve the ambiguity because our ambiguity is only a binary one. For example,
the rightmost object “red-leaf” in Fig. 16 has a concave
shape that is correctly and automatically recovered byusing the outward-looking boundary normals.
We make quantitative evaluations for these results.
First, the average angular error of “metallic hemisphere”(the leftmost case in Fig. 16), which has a known shape,is 3.45◦ in the second row, and is 9.99◦ in the third
row. The latter one has a larger error partially because
intensity matching is more sensitive to intensity loss insaturated areas, as we did not purposely use a highbit-depth camera. As for other objects, we do not haveground truth of their surface normals. However, we can13Est. normal Est. normal One of input
with ref. data w/o ref. data
imagesPREPRINT OF A REGULAR PAPER @PAMI2015

Estimated error:

8.27°

9.15 °

6.73 °

11.30 °

10.43 °

8.57 °

Fig. 16. Result of real objects using our methods without (middle row) and with (bottom row) reference data.

Estimated
error:
8.30°

Estimated
error:
10.63°

Fig. 17. Real scenes without occluding boundaries.
assess their light source estimation accuracy since weknow their ground truth. Because the average error oflight sources are similar to (usually larger than) that of
surface normals as shown in Table 1 and Table 5, they canbe used to estimate the surface normal accuracy. Theseestimated errors are shown in the bottom of Fig. 16.
We also test a few scenes that have no occluding
boundaries, as shown in Fig. 17. As the results show,
our methods, no matter with or without reference data,
can estimate surface normals purely from the images ofthe surface patches. This shows the advantage of the
proposed methods over the prior methods [3], [36].

6

D ISCUSSION AND CONCLUSION

In this paper, we propose photometric stereos methods
that handle both uncalibrated light sources and general
unknown reflectances. Our approach uses intensity profiles, and it is built upon our findings on the relations between the geodesic distance of intensity profiles and theangular distance of surface normals, and also between
the intensity distribution of intensity profiles and surface
reflectance property. We investigate these relations anddevelop methods that recover surface normals under
two conditions: purely from intensity profiles and by using reference data. Experimental results show impressive
results in surface normal estimation with unknown lightsources and general reflectances.
The assumptions about uniform lights and reflectanceare somehow restrictive. Although we have shown thatthey can be relaxed to some degree without significantly
affecting the results, further investigation on how to
remove them can be a good future research direction.

R EFERENCES[1]

S. Koppal and S. Narasimhan, “Clustering appearance for scene
analysis,” in Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), 2006, pp. 1323–1330.
[2] A. Hertzmann and S. Seitz, “Example-based photometric stereo:
shape reconstruction with general, varying BRDFs,” IEEE Trans.
on Pattern Analysis and Machine Intelligence, vol. 27, no. 8, pp. 1254–
1264, 2005.
[3] I. Sato, T. Okabe, Q. Yu, and Y. Sato, “Shape reconstruction basedon similarity in radiance changes under varying illumination,” inProc. of Int’l Conf. on Computer Vision (ICCV), 2007, pp. 1–8.
[4] G. Oxholm, P. Bariya, and K. Nishino, “The scale of geometrictexture,” in Proc. of European Conf. on Computer Vision (ECCV),
2012.
[5] F. Lu, Y. Matsushita, I. Sato, T. Okabe, and Y. Sato, “Uncalibrated
photometric stereo for unknown isotropic reflectances,” in Proc.
of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR),
2013.
[6] W. Matusik, H. Pfister, M. Brand, and L. McMillan, “A data-driven
reflectance model,” in Proc. of ACM SIGGRAPH, 2003, pp. 27–31.
[7] P. N. Belhumeur, D. J. Kriegman, and A. L. Yuille, “The bas-relief
ambiguity,” Int’l Journal of Computer Vision, vol. 35, no. 1, pp. 33–
44, 1999.
[8] S. Barsky and M. Petrou, “The 4-source photometric stereo technique for three-dimensional surfaces in the presence of highlights
and shadows,” IEEE Trans. on Pattern Analysis and Machine Intelligence, vol. 25, no. 10, pp. 1239–1252, 2003.
[9] D. Miyazaki, K. Hara, and K. Ikeuchi, “Median photometric stereo
as applied to the segonko tumulus and museum objects,” Int’l
Journal of Computer Vision, vol. 86, no. 2, pp. 229–242, 2010.
[10] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma,“Robust photometric stereo via low-rank matrix completion andrecovery,” in Proc. of Asian Conf. on Computer Vision (ACCV), 2010,
pp. 703–717.
[11] S. Ikehata, D. Wipf, Y. Matsushita, and K. Aizawa, “Robust
photometric stereo using sparse regression,” in Proc. of IEEE Conf.
on Computer Vision and Pattern Recognition (CVPR), 2012, pp. 318–
325.
[12] T. Wu, K. Tang, C. Tang, and T. Wong, “Dense photometric stereo:
A markov random field approach,” IEEE Trans. on Pattern Analysis
and Machine Intelligence, vol. 28, no. 11, pp. 1830–1846, 2006.
[13] N. Alldrin and D. Kriegman, “Toward reconstructing surfaces
with arbitrary isotropic reflectance : A stratified photometricstereo approach,” in Proc. of Int’l Conf. on Computer Vision (ICCV),
2007, pp. 1–8.
[14] M. Holroyd, J. Lawrence, G. Humphreys, and T. Zickler, “A photometric approach for estimating normals and tangents,” ACM
Transactions on Graphics, vol. 27, p. 133, 2008.
[15] T. Higo, Y. Matsushita, and K. Ikeuchi, “Consensus photometricstereo,” in Proc. of IEEE Conf. on Computer Vision and Pattern
Recognition (CVPR), 2010, pp. 1157–1164.
[16] B. Shi, P. Tan, Y. Matsushita, and K. Ikeuchi, “Elevation angle fromreflectance monotonicity: Photometric stereo for general isotropic
reflectances,” in Proc. of European Conf. on Computer Vision (ECCV),
2012.

PREPRINT OF A REGULAR PAPER @PAMI2015[17] P. Tan, L. Quan, and T. Zickler, “The geometry of reflectance symmetries,” IEEE Trans. on Pattern Analysis and Machine Intelligence,
vol. 33, no. 12, pp. 2506–2520, 2011.
[18] D. Goldman, B. Curless, A. Hertzmann, and S. Seitz, “Shape andspatially-varying brdfs from photometric stereo,” in Proc. of Int’l
Conf. on Computer Vision (ICCV), 2005, pp. 341–348.
[19] H. Chung and J. Jia, “Efficient photometric stereo on glossy
surfaces with wide specular lobes,” in Proc. of Int’l Conf. on
Computer Vision (ICCV), 2008, pp. 1–8.
[20] S. Nayar, K. Ikeuchi, and T. Kanade, “Determining shape andreflectance of hybrid surfaces by photometric sampling,” IEEE
Transactions on Robotics and Automation, vol. 6, no. 4, pp. 418–431,
1990.
[21] N. Alldrin, T. Zickler, and D. Kriegman, “Photometric stereo with
non-parametric and spatially-varying reflectance,” in Proc. of IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR), 2008, pp.
1–8.
[22] P. Favaro and T. Papadhimitri, “A closed-form solution to uncalibrated photometric stereo via diffuse maxima,” in Proc. of IEEE
Conf. on Computer Vision and Pattern Recognition (CVPR), 2012, pp.
821–828.
[23] O. Drbohlav and R. Šára, “Specularities reduce ambiguity of
uncalibrated photometric stereo,” in Proc. of European Conf. on
Computer Vision (ECCV), 2002, pp. 644–645.
[24] O. Drbohlav and M. Chaniler, “Can two specular pixels calibrate
photometric stereo?” in Proc. of Int’l Conf. on Computer Vision(ICCV), 2005, pp. 1850–1857.
[25] R. Basri, D. Jacobs, and I. Kemelmacher, “Photometric stereo
with general, unknown lighting,” Int’l Journal of Computer Vision,
vol. 72, no. 3, pp. 239–257, 2007.
[26] N. Alldrin, S. Mallick, and D. Kriegman, “Resolving the generalized bas-relief ambiguity by entropy minimization,” in Proc. ofIEEE Conf. on Computer Vision and Pattern Recognition (CVPR),
2007.
[27] M. Chandraker, F. Kahl, and D. Kriegman, “Reflections on thegeneralized bas-relief ambiguity,” in Proc. of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR), 2005, pp. 788–795.
[28] B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan, “Self-calibrating
photometric stereo,” in Proc. of IEEE Conf. on Computer Vision and
Pattern Recognition (CVPR), 2010, pp. 1118–1125.
[29] Z. Wu and P. Tan, “Calibrating photometric stereo by holisticreflectance symmetry analysis,” in Proc. of IEEE Conf. on Computer
Vision and Pattern Recognition (CVPR), 2013.
[30] Z. Zhou and P. Tan, “Ring-light photometric stereo,” in Proc. ofEuropean Conf. on Computer Vision (ECCV), 2010, pp. 265–279.
[31] T. Papadhimitri and P. Favaro, “A new perspective on uncalibrated photometric stereo,” in Proc. of IEEE Conf. on Computer Vision
and Pattern Recognition (CVPR), 2013.
[32] W. Silver, “Determining shape and reflectance using multiple
images,” Master’s thesis, M.I.T., Dept. of Electrical Engineeringand Computer Science, 1980.
[33] P. Ren, J. Wang, J. Snyder, X. Tong, and B. Guo, “Pocket reflectometry,” in Proc. of ACM SIGGRAPH, 2011, pp. 45:1–45:10.
[34] A. Georghiades, “Incorporating the torrance and sparrow modelof reflectance in uncalibrated photometric stereo,” in Proc. of Int’l
Conf. on Computer Vision (ICCV), 2003, pp. 816–823.
[35] M. Chandraker, J. Bai, and R. Ramamoorthi, “A theory of differential photometric stereo for unknown isotropic brdfs,” in Proc.
of IEEE Conf. on Computer Vision and Pattern Recognition (CVPR),
2011, pp. 2505–2512.
[36] T. Okabe, I. Sato, and Y. Sato, “Attached shadow coding: Estimating surface normals from shadows under unknown reflectanceand lighting conditions,” in Proc. of Int’l Conf. on Computer Vision(ICCV), 2009, pp. 1693–1700.
[37] J. Tenenbaum, V. De Silva, and J. Langford, “A global geometricframework for nonlinear dimensionality reduction,” Science, vol.
290, no. 5500, pp. 2319–2323, 2000.
[38] H. Hayakawa, “Photometric stereo under a light source with
arbitrary motion,” JOSA A, vol. 11, no. 11, pp. 3079–3089, 1994.
[39] A. Yuille, D. Snow, R. Epstein, and P. Belhumeur, “Determining
generative models of objects under varying illumination: Shapeand albedo from multiple images using svd and integrability,”Int’l Journal of Computer Vision, vol. 35, no. 3, pp. 203–222, 1999.
[40] J. B. Kruskal and M. Wish, Multidimensional scaling. Sage, 1978,
vol. 11.
[41] “Thestanford
3dscanning
repository,”
http://graphics.stanford.edu/data/3Dscanrep/.

14[42] R. Woodham, “photometric method for determining surface orientation from multiple images,” Optical Engineering, vol. 1, no. 7,
pp. 139–144, 1980.

Feng Lu received the BS and MS degrees inAutomation from Tsinghua University in 2007
and 2010, and the PhD degree in information
science and technology from the University ofTokyo in 2013 respectively. He is currently aproject research associate in the Institute ofIndustrial Science, the University of Tokyo. His
research interests include human gaze analysis,
photometric stereo and reflectance analysis.

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 theeditorial board of IEEE Transactions on PatternAnalysis and Machine Intelligence (TPAMI), International Journal ofComputer Vision (IJCV), IPSJ Journal of Computer Vision and Applications (CVA), The Visual Computer Journal, and Encyclopedia of
Computer Vision. He served/is serving as a Program Co-Chair of PSIVT
2010, 3DIMPVT 2011, ACCV 2012, and ICCV 2017, and a General CoChair for ACCV 2014. He is a senior member of IEEE.
Imari Sato received her B.S. degree in policy
management from Keio University in 1994 and
M.S. and Ph.D. degrees in interdisciplinary information studies from the University of Tokyoin 2002 and 2005. She is an Associate Professor at the National Institute of Informatics. She
serves concurrently as Associate Professor at
the Graduate university for Advanced Studies
and as Visiting Associate Professor at Department of Information Processing, Tokyo Instituteof Technology. Her research interests include
physics-based vision and image-based modeling.
Takahiro Okabe received the B.S. and M.S.
degrees in physics, and the Ph.D. degree in
information science and technology from the
University of Tokyo, Japan, in 1997, 1999 and2011 respectively. After working at the Instituteof Industrial Science, the University of Tokyo, hejoined Kyushu Institute of Technology, Japan, as
an associate professor in 2013. His research interests include computer vision, image processing, pattern recognition, and computer graphics,
in particular their physical and mathematical aspects.
Yoichi Sato received the BS degree from the
University of Tokyo in 1990, and the MS and PhD
degrees in robotics from the School of Computer
Science, Carnegie Mellon University, in 1993
and 1997, respectively. He is a professor at theInstitute of Industrial Science, the University ofTokyo. His research interests include physicsbased vision, reflectance analysis, and gaze andgesture analysis. He served/is serving in several conference organization and journal editorial
roles including IEEE Transactions on PatternAnalysis and Machine Intelligence, International Journal of Computer
Vision, IPSJ Journal of Computer Vision and Applications, ECCV2012
Program Co-Chair and MVA2013 General Chair.