Photometric Stereo via Discrete Hypothesis-and-Test Search
Kenji Enomoto1

Michael Waechter1 Kiriakos N. Kutulakos2 Yasuyuki Matsushita1
1
Osaka University 2 University of Toronto

Abstract

BRDFs & known lights…

…

In this paper, we consider the problem of estimating
surface normals of a scene with spatially varying, general
BRDFs observed by a static camera under varying, known,
distant illumination. Unlike previous approaches that are
mostly based on continuous local optimization, we cast the
problem as a discrete hypothesis-and-test search problem
over the discretized space of surface normals. While a naı̈ve
search requires a significant amount of time, we show thatthe expensive computation block can be precomputed in ascene-independent manner, resulting in accelerated inference for new scenes. It allows us to perform a full searchover the finely discretized space of surface normals to determine the globally optimal surface normal for each scene
point. We show that our method can accurately estimate
surface normals of scenes with spatially varying different
reflectances in a reasonable amount of time.

Surface normal
candidates

Hypothesize - and -Test

Target scene

Figure 1: An overview of our approach. We hypothesize
a surface normal and test whether it can explain the target
measurements. By conducting the hypothesis-and-test for
all possible surface normals, our method is able to find a
globally optimal surface normal.
results in a relatively small number of surface normal candidates. For example, even if we discretize the angles in
one-degree intervals, it results in 32, 400 = 360×90 normal
candidates. Our method uses each surface normal candidate
as hypothesis and tests its suitability to the image formation
model; thus, we call it a hypothesis-and-test search strategy. In this manner, our method searches for the globally
optimal surface normal from all (discretized) possible ones.
To alleviate the issue of computing cost in our discrete
search, we developed a precomputation method that performs expensive computations in a scene-independent manner prior to the inference for a new scene. To deal witha diverse set of reflectances, we use a non-parametric, discrete table of bidirectional reflectance distribution functions(BRDFs), whose axes are the space of surface normals, light
directions, and materials, for a fixed viewing direction. The
table of BRDFs, which we call a BRDF tensor, can containan arbitrary number of materials, and importantly, the number of reference materials considered in the BRDF tensor
does not influence the computation time during inference.
Our method is motivated by the success of examplebased [9] and virtual exemplar-based [11] approaches.
However, unlike example-based methods, we do not require

1. Introduction
Photometric stereo recovers fine surface details in theform of surface normals from images taken by a static camera under varying lightings. While traditional photometricstereo methods [25] assume Lambertian reflectance or simplified parametric reflectance models, it is understood that
their deviation from real-world reflectances introduces errors in surface normal estimates. In the past, other studies
used more sophisticated reflectance models [1, 15, 21, 13, 4]
for more accurate surface normal recovery; however, theygenerally encounter an issue of non-convex optimization in
determining the surface normals. The problem is rooted in
the fact that these methods frame the estimation problem asa continuous optimization problem.
In this paper, we cast surface normal estimation as a discrete hypothesis-and-test search problem. Instead of treating surface normals to be estimated as a continuous quantity, our method finely discretizes the space of surface normals and finds the best surface normal by a hypothesis-andtest search. Since a surface normal vector has only two degrees of freedom (a unit 3D vector) represented by its azimuth and elevation angles in a hemisphere, discretization1

placing a reference object in the scene. Also, unlike the virtual exemplar-based method that performs a continuous local search using a non-convex objective function, we treat
the problem as a discrete search problem and perform an
exhaustive search over the discretized space. It enables us
to find a globally optimal surface normal within the boundsof our objective function.
The chief contributions of this paper are twofold. First,
we propose a discrete hypothesis-and-test search strategy
for photometric stereo. By finely discretizing the space of
surface normals, our method finds the globally optimal surface normal through exhaustive search. Second, we showthat expensive computation can be performed prior to thesurface normal estimation, allowing the global hypothesisand-test search to work in a reasonable amount of time. We
assess the accuracy of the proposed method using both synthetic and real-world data and show its favorable performance in determining surface normals of a scene. In particular, the proposed method achieves a stable estimate, i.e.,
superior average/variance of mean angular error over a diverse set of materials.

2. Related work
Calibrated photometric stereo methods for diverse materials can be roughly divided into three categories; modelbased, learning-based, and example-based approaches. In
the following, we discuss the corresponding related works.

Model-based photometric stereo Model-based approaches use parametric expressions for BRDFs andthe model parameters including the surface normal are
estimated, typically, via optimization. Key for thesemodel-based methods is the choice of a parametric BRDF
model. Woodham’s original work [25] assumed Lambertian reflectance, which allows using convex least-squares
optimization to determine surface normals and albedos.
Parametric modeling of non-Lambertian BRDFs is actively
studied, particularly in the graphics community. For
example, the Blinn-Phong model [23], the TorranceSparrow model [7], the Ward model [6], the specular spikemodel [5, 26], and a microfacet BRDF with ellipsoidal
normal distributions [4] have been developed. However,each of these models is limited to a class of materials, andsuch models are highly nonlinear, resulting in non-convex
photometric stereo problems. Thus, some recent methods
use a bivariate function instead. For representing lowfrequency reflectances, Shi et al. [21] use a bi-polynomial
function and Ikehata and Aizawa [13] use a sum of lobes
with unknown center directions. Although these modelbased approaches can be used in a relatively wide range ofmaterials, there are always problematic materials.

Learning-based photometric stereo Recently, deeplearning-based photometric stereo methods have been proposed. They learn a mapping from measured intensities
under known lightings to surface normals using a neural network [19, 3, 13]. These methods showed promising results on various scenes owing to the network being
trained with diverse shapes and materials. Santo et al. [19]and Chen et al. [3] created a training dataset by rendering the Blobby [14] and Sculpture [24] shape datasets with
100 BRDFs from the MERL dataset [16]. Ikehata [12]also introduced a training dataset for calibrated photometricstereo, called CyclesPS dataset, containing several objects
with a diverse set of materials rendered using Disney’s principled BSDF [2]. Learning-based methods require a large
number of shapes & BRDFs for training. When they areapplied to new scenes with shapes & BRDFs that are very
different from the training data, it is necessary to do costly
re-training of the network.
Example-based photometric stereo Example-based
photometric stereo relies on the concept of orientationconsistency [9], i.e., two surfaces with the same surface
normal and BRDF will have the same appearance under
the same illumination. An early work along this direction
is found in Horn and Ikeuchi [10]. In example-based
approaches, a reference object with known surface normals
is placed in a target scene. Further, the BRDF of thereference object is assumed to be the same as that of thetarget object. Then, a surface normal is recovered for eachpoint of the target object by searching the corresponding
pixel intensity of the reference object that best matches thetarget’s appearance. To relax the assumption of identical
BRDF between reference and target, Hertzmann andSeitz [9] introduced two reference objects, a diffuse and a
specular sphere, placed in the target scene, and approximate
the target BRDF by a non-negative linear combinationof the reference BRDFs. Although this approach makes
example-based photometric stereo applicable to morediverse materials, it is still inaccurate to approximate adiverse set of materials by a linear combination of two
BRDFs. In addition, in many practical applications it isundesirable to place reference objects in a target scene.
Hui and Sankaranarayanan [11] introduced virtual
exemplar-based photometric stereo that performs examplebased photometric stereo without actually introducing reference objects into a target scene. They render virtualreference spheres under the target scene illumination with
MERL BRDFs [16] and assume that the target BRDF lies in
the non-negative span of the MERL BRDFs. In the virtual
exemplar-based approach, however, there are many timeconsuming processes such as rendering virtual spheres,
an iterative optimization for solving a non-negative leastsquares problem, and searching over all possible surface

normals. To reduce the computation cost, they proposed
an efficient search algorithm which however eliminates theguarantee of finding the optimal solution.
Our method shares the assumption that the target BRDF
can be represented by a combination of several reference
BRDFs. However, we cast the problem as a discrete
hypothesis-and-test search problem, which gives a guarantee of reaching the globally optimal solution within thebound of the objective function. Additionally, our method
enables search all surface normal candidates in reasonable
time owing to an efficient precomputation.

3. Image formation and problem statement
Suppose a surface point with a unit surface normaln ∈ S 2 ⊂ R3 is illuminated by an incoming directional
light l ∈ S 2 , without ambient lighting or global illumination effects such as cast shadows or inter-reflections. When
this surface point is observed by a camera with linear response, the measured intensity m ∈ R+ can be written asm ∝ ρ(n, l) max(n> l, 0),(1)where ρ(n, l) : S 2 × S 2 → R+ is a general isotropic bidirectional reflectance distribution function (BRDF).
In calibrated photometric stereo, a static camera records
multiple, say L0 , measurements {m1 , . . . mL0 } for each surface point under various light directions {l1 , . . . lL0 }. Then,
Eq. (1) can be written in matrix form as 
max(n> l1 , 0)m1
 .  ..
 .. ∝.


mL0| {z } |0max(n> lL0 , 0)

0

}|

{z

m



E



ρ(n, l1 )

..
,
.

ρ(n, lL0 ){z}ρwhere m is a measurement vector, E is a diagonal irradiance matrix, and ρ is a reflectance vector. We model thereflectance ρ by a linear combination of BRDF basis vectors in a similar manner to Hertzmann et al. [9], and Huiand Sankaranarayanan [11]. By stacking M known BRDF
basis vectors in a BRDF basis matrix B, ρ can be written asρ1 (n, l1 ) . . . ρM (n, l1 )


..
..
..
ρ=
 c,
.
.
.
|

ρ1 (n, lL0 ) . . . ρM (n, lL0 ){z}Bwhere c = [c1 , . . . , cM ]> is a BRDF coefficient vector.
With this, the image formation model can be simplified tom = EBc.

(2)

Problem statement Our goal is to find the optimal surface normal n and BRDF coefficients c for each surface

surface normal 𝒏hypothesized normal 𝒏𝑖

light 𝒍material
light 𝒍

known
light
directions…

material
BRDF tensor 𝒯 ∈ ℝ𝑁×𝐿×𝑀+

sampled BRDF′
matrix 𝑩𝑖 ∈ ℝ𝐿 ×𝑀

Figure 2: Starting from the BRDF tensor T that represents
reflectances for a comprehensive set of light directions, surface normals and materials (BRDFs), we slice out a sampled
BRDF matrix Bi for a set of known light directions and ahypothesized surface normal ni . The column space of Bi is
the space of reflectances over all possible materials for the
hypothesized normal under the known light directions.
point, given observations m and associated light directions{l1 , . . . lL0 } based on the model of Eq. (2). The irradiancematrix E and BRDF matrix B are functions of the surface
normal n.

4. Proposed method
Our method casts the photometric stereo problem as
a discrete search where the space of surface normals is
discretized. We hypothesize a surface normal n and test
whether it (approximately) satisfies the image formation
model of Eq. (2). By conducting this hypothesis-and-test
for all possible surface normals, our method is able to find a
globally optimal surface normal n and the associated BRDF
coefficients c that best satisfy Eq. (2).

4.1. Hypothesis-and-test strategy
Let N = {ni | i = 1, . . . , N } be the discretized spaceof surface normals, which we call a set of surface normal
candidates. We prepare a tensor representation for diverse
BRDFs whose axes are (1) surface normals, (2) light directions, and (3) materials. Suppose the spaces of surface
normals and light directions are discretized into N and L
bins, respectively, and there are M distinct BRDFs. Then,×L×M
the BRDF tensor T can be defined as T ∈ RN
(see+
the left of Fig. 2).
For simplicity, let us assume that the BRDF tensor contains the actual light directions of the observed scene. If wehypothesize a certain surface normal ni ∈ N for a scene
point, using L0 ≤ L known light directions, we can slice a0
×M
sampled BRDF matrix Bi ∈ RL
from the BRDF tensor+
T along the hypothesized surface normal ni and a set of L0
known light directions as illustrated in Fig. 2. We can also-form an irradiance matrix Ei for the hypothesized surface
normal ni . Using Bi and Ei instead of B and E, Eq. (2)becomesran(𝐃#)𝐙$𝐦
𝐙#

ran(𝐃$)𝐦Ω ⊂ ℝ,𝐦def
m ' Ei Bi c =
Di c,

𝐧$

𝐧#

0

where Di (= Ei Bi ) ∈ RL ×M . For the overdeterminedcase L0 > M , the least-squares solution for the BRDF coefficients c that best explains the measurements, is−1 >c∗ = D > D
D m = D†i m,where D†i is Di ’s pseudo-inverse. The estimated BRDF coefficients c∗ are least-squares optimal for the hypothesized
normal ni and the space of sampled BRDFs Bi . We cantest the validity of the hypothesized ni by evaluating the `2
reconstruction error as2

ei = km − Di c∗ k2 .

(3)

Therefore, the optimal surface normal n∗ can be found as
the minimizer of the following objectiven∗ = ni∗ ,i∗ = argmin ei .

(4)i∈{1,...,N }A naı̈ve implementation may require a significant computation effort for solving this problem. We therefore introduce
an efficient precomputation strategy now.

Figure 3: Geometric interpretation of the reconstruction er2
ror. The reconstruction error of measurements kZi mk2 canbe seen as distance between the measurement vector m andthe subspace spanned by Di in the L0 -dimensional space Ω.
As illustrated in Fig. 3, a measurement vector m exists in an L0 -dimensional space Ω. The column vectorsof Di span a rank(Di )-dimensional subspace in Ω, andthe measurement reconstructions Di c∗ = Di D†i m residein this subspace. Thus, geometrically, the reconstruction
2
error kZi mk2 can be considered as the distance between
the measurement vector m and the subspace spanned by
Di . From this perspective, if rank(Di ) = L0 , the columnsof Di span the entire Ω; therefore, the reconstruction error
becomes always zero regardless of the correctness of the
hypothesized surface normal n.
To avoid such a situation, we reduce the column dimension of Di from M to M 0 < L0 if L0 < M . Specifically,
we apply singular value decomposition (SVD) to Di as
Di = USV> ,4.2. Precomputation
The reconstruction error ei for the hypothesized surface
normal ni in Eq. (3) can be further simplified as2
ei = km − Di c∗ k2 =

=





I − Di D†i m
0m − Di D†i m
2

2

def
kZi mk2 ,=2and reduce the dimensionality of Di using the first M 0 singular vectors / values as>
Di ← UM 0 SM 0 VM
0,

2
2

(5)

0

where Zi (= I − Di D†i ) ∈ RL ×L is uniquely determined
given a hypothesized normal ni . We can precompute a setof {Zi } for all normal candidates in N and at test time we
simply need to assess the magnitude of Zi m for all i.
This precomputation happens only once and the result
can be used for any new scene as long as the light configuration is unchanged.

4.3. Dimensionality reduction of Di
In principle, minimizing Eq. (4) gives us a correct solution for the surface normal n. In practice, however, weneed to pay attention to the dimension and range of matrix0
Di ∈ RL ×M . Specifically, when L0 < M or m ∈ ran(Di )(the range of Di ) for all Di , there exists one or more BRDF
coefficient vectors c∗ that makes all reconstruction errors{ei } zero.

>
where UM 0 , SM 0 , and VM
0 are the truncated singular vectors / values.
We empirically found that the proper value of M 0 is related to the noise level in the observations. In Sec. 5.4,
we examine the accuracy of surface normal estimation withvarying M 0 and discuss the choices for M 0 .

5. Experiments
This section describes the results of experiments withsynthetic and real-world data. We further discuss the computation time, the effect of dimensionality reduction and thediscretization of the space of light directions. We begin withdescribing the construction of the BRDF tensor and the synthetic and real-world datasets that we use for evaluation.
BRDF tensor: The BRDF tensor is constructed fromthree components; materials, surface normal candidates,
and light directions. As materials, we used the MERLBRDF dataset [16] which consists of 100 distinct BRDFs

Model-based
CNN-based
Ours

violet-rubber
light-red-paintblue-rubber
dark-red-paintorange-paint
pure-rubber
yellow-paintgold-paintyellow-plastic
pink-felt
pink-fabricneoprene-rubber
teflonred-plastic
pvc
pink-plastic
delringray-plastic
violet-acrylic
specular-blue-phenolicsilicon-nitrade
silver-paintwhite-diffuse-bballpolyethylene
white-paintgreen-plastic
pearl-paintcherry-235
white-acrylic
specular-white-phenoliccolonial-maple-223
alumina-oxidegreen-metallic-paintred-specular-plastic
fruitwood-241
specular-maroon-phenolic
gold-metallic-paint3
purple-paintspecular-green-phenolicpink-jasper
green-acrylic
yellow-matte-plastic
red-fabric2
specular-violet-phenolicnylon
red-phenolicdark-specular-fabricipswich-pine-221
two-layer-gold
polyurethane-foam
beige-fabricpickled-oak-260
specular-red-phenolicmaroon-plastic
natural-209
aventurnine
white-fabricwhite-marble
blue-acrylic
specular-orange-phenolic
two-layer-silver
red-metallic-paintgold-metallic-paint2
specular-yellow-phenolic
blue-fabricgreen-fabrichematitealum-bronze
black-oxidized-steel
black-soft-plastic
dark-blue-paintyellow-phenolic
green-metallic-paint2
black-obsidian
red-fabric
blue-metallic-paint2
specular-black-phenoliccolor-changing-paint2
white-fabric2
nickelcolor-changing-paint3
silver-metallic-paintgold-metallic-paintpink-fabric2
light-brown-fabric
blue-metallic-paintgreen-latex
silver-metallic-paint2
chromeblack-phenolic
ss440
special-walnut-224
color-changing-paint1
brass
black-fabricchrome-steel
tungsten-carbide
steel
aluminium
grease-covered-steel

Mean angularerror (degrees)

25
20
15
10
5
0

Materials

Figure 4: Evaluation on the MERL sphere dataset. We compared our method against a model-based method [13] and a
CNN-based method [12]. The average/variance of the mean angular errors over all materials are as follows: Ours: 2.23/3.24,
model-based method [13]: 4.31/29.43, CNN-based method [12]: 3.15/9.19.
including diffuse, specular, and metallic materials. For
surface normal candidate sampling we followed Hui’s
method [11] and obtained 20001 candidates by 0.5◦ equiangular sampling over the hemisphere [8]. In all experimentsof this paper, we assume that the BRDF tensor contains the
known light directions.
MERL sphere: We created a MERL sphere dataset,
which is synthetic data rendered with 100 isotropic BRDFs
from the MERL BRDF database [16]. We generated images of a sphere object illuminated from 20, 40, 60, 80, 100
different known directions1 , respectively. We also created anoisy MERL sphere dataset by adding signal-independentand signal-dependent noise [17] to the MERL
sphere√dataset. The noise model is m̃ = m + (µ + λ m)X wherem̃ and m are image signals with and without noise, µand λ are weighting factors for signal-independent andsignal-dependent noise, respectively, and X is a N (0, 1)distributed random variable.
MERL bunny: We rendered the Stanford bunny withspatially varying materials using 11 distinct MERL BRDFs
as test BRDFs and 80 different known light directions.
Real-world benchmark: We took an existing real-world
dataset, the DiLiGenT dataset [20], which contains 10 real
objects of general reflectance illuminated from 96 different known directions. This dataset provides ground truth
surface normal maps for all objects measured by highprecision laser scanning, enabling quantitative evaluation.

5.1. Evaluation on synthetic data
We performed experiments on synthetic data to confirm
that our method works with diverse materials. Since there
is no global illumination in the MERL sphere & bunnydataset, we can evaluate only the ability of our method toadapt to diverse materials. To support that our method isapplicable to a more diverse set of materials than existing methods, we compared it against a model-based [13], a
CNN-based [12], and a dictionary-based method [11]. Additionally, we investigated the accuracy of our method withvarying surface normal candidates and lights, and compared
it to the dictionary-based method [11]. All following experiments were performed in the 80 lights setting unless otherwise specified.

MERL sphere: We performed a comparative experiment
between our method, the model-based method [13], and the
CNN-based method [12], and verified that our method ismore accurate and stable on diverse material data. For thematerials we used a leave-one-out scheme, testing it on one
MERL BRDF while constructing the BRDF tensor from theremaining 99 BRDFs.2
The results are shown in Fig. 4. The average/variance ofthe mean angular errors over all materials are 2.23/3.24 forour method, 4.31/29.43 for the model-based method [13],and 3.15/9.19 for the CNN-based method [12]. This result
suggests that our method is more accurate and stable thanexisting methods on most materials in the MERL database.
2 Due

1 We used the generalized spiral points algorithm for uniformly dis-tributing lights on the hemisphere [18]. See the supplementary material
for details of the light distributions.

to the high computational costs of the dictionary-based
method [11], a leave-one-out scheme on it with the author’s implementation cannot be performed in a reasonable amount of time. We do a comparison with it in the following experiments.

20°+
10°
0°
Ground truth &observation

1.36/2.24

2.16/4.07

4.98/36.13

1.37/2.39

Ours

CNN-based [12]Model-based [13]Dictionary-based [11]Figure 5: Visual comparison for the Stanford bunny with spatially varying materials. We show estimated normals and angularerror maps of our and 3 comparison methods. Values under the error maps indicate the mean/variance of angular errors.
MERL bunny: We evaluated our method, the modelbased [13], the CNN-based [12], and the dictionary-based
method [11] on the Stanford bunny with spatially varying materials. Our BRDF tensor and the BRDF dictionary of the dictionary-based method are constructed from80 MERL BRDFs that do not include the 11 test BRDFs
mentioned above.
In Fig. 5 we show each method’s estimated surface normals and angular error map. The values under the error
maps indicate the mean/variance of their angular errors.
This result suggests that our method can be applied to materials that are difficult for the CNN-based and the modelbased method. Additionally, from the low variance of angular errors we can see that our method is more stable thanthe existing methods on diverse materials.
Varying number of surface normal candidates andlights: We investigate the angular errors of the estimatedsurface normals for varying numbers of surface normal
candidates and lights, and compare our method to the existing search-based method [11]. In the experiment witha varying number of surface normal candidates, we use
20001, 10001, 1501, 251 surface normal candidates yieldedby 0.5◦ , 1◦ , 3◦ , 5◦ equiangular sampling, respectively. In
the experiment with a varying number of lights, we use the
MERL sphere dataset with 20, 40, 60, 80, 100 lights. Dueto the high computational costs of the existing dictionarybased method [11], we perform a 5-fold cross-validation,using 80 BRDFs for computing our method’s BRDF tensor
and the dictionary method’s material dictionary [11], andtesting on the MERL sphere dataset with the remaining 20
BRDFs.

Table 1: Mean angular errors for estimated surface normals
in degrees for varying numbers of normal candidates.
number of normal candidates N
20001 10001 1501
251
Ours
HS17 [11]

2.25
2.34

2.31
2.41

2.86
3.00

4.62
4.68

Table 2: Mean angular errors for estimated surface normals
in degrees for varying numbers of lights.

Ours
HS17 [11]100number of lights L0
80
60
40

20

2.22
2.30

2.25
2.34

2.69
2.81

2.28
2.39

2.35
2.49

Tables 1 and 2 show the experimental results for varying numbers of surface normal candidates and lights, respectively. Both results suggest that our method consistently estimates surface normals more accurately than theexisting method. One of the reasons is a difference in the
search strategies. Hui and Sankaranarayanan [11] proposeda search strategy that avoids searching parts of the surface
normal space since their method is computationally expensive due to iterative optimization. In contrast, our method
searches over all surface normal candidates in reasonable
computation time facilitated by the precomputation.

5.2. Evaluation on real-world data
We show quantitative comparisons on the DiLiGenT
dataset in Tab. 3 where we compare our method with existing methods in terms of mean angular error. Here, ourTable 3: Comparisons on the DiLiGenT dataset. We apply a least-squares method as baseline. The values represent mean
angular error (MAE).
Ball

Bear

BuddhaCat

Cow

GobletHarvest

Pot1

Pot2

Reading

Avg.

Ours

1.58

6.38

13.69

6.30

7.80

11.42

18.74

6.67

7.26

15.49

9.53

SI18 [12]HS18 [11]CH18 [3]TM18 [22]ST14 [21]IA14 [13]Baseline

2.20
1.33
2.80
1.47
1.74
3.34
4.10

4.10
5.58
7.60
5.79
6.12
7.11
8.39

7.90
8.48
7.90
10.36
10.60
10.50
14.90

4.60
4.88
6.20
5.44
6.21
6.74
8.41

7.90
8.23
7.30
6.32
13.90
13.10
25.60

7.30
7.57
8.60
11.47
10.10
9.71
18.50

13.90
15.80
15.90
22.59
25.40
26.00
30.60

5.40
5.16
7.10
6.09
6.51
6.64
8.89

6.00
6.41
7.30
7.76
8.78
8.77
14.70

12.60
12.08
13.30
11.03
13.60
14.20
19.80

7.20
7.55
8.40
8.83
10.30
10.60
15.40

Cat

20°+
10°

Reading0°

Observation

Ground truth

CNN-based [12]Ours

Figure 6: Visual comparison for the non-convex objects Cat and Reading. We show estimated normals and angular error
maps of our method and a state-of-the-art CNN-based method [12].
method does not achieve the best results. This is considered
to be due to factors not modeled in our method, namely
cast shadows or inter-reflections. Figure 6 shows visualcomparisons between our method and Ikehata’s CNN-based
method [12] for two non-convex objects with cast shadows and inter-reflections. Our method causes a large angular error in pixels where cast shadows or inter-reflections
are likely to occur. However, in convex parts our method
outperforms the CNN-based state-of-the-art method and estimates the surface normals well for the Reading object’s
complex specular material. We show more visual comparisons in the supplementary material.

5.3. Computation cost
For inference, our method evaluates the reconstruction
2
error kZi mk2 in Eq. (5) for each surface normal candidate
ni ∈ N . All matrices Zi are precomputed; therefore, atinference time we only need to evaluate the reconstruction
error of each ni and find the minimizer. The dimension of00matrix Zi ∈ RL ×L only depends on the number of lights
0
L , but not the number of materials. Table 4 summarizes

Table 4: Computation time in seconds for surface normal
estimation of one pixel with 20001 surface normal candidates, measured on one core of an Intel Xeon CPU E5-2680
v4 @ 2.40GHz.

Time [sec.]

100

number of lights L0
80
60
40

20

0.31

0.21

0.024

0.15

0.062

the computation time per pixel for 20001 surface normal
candidates and a varying number of lights L0 . The computation is highly parallelizable, e.g., by pixel-wise or normal
candidate-wise parallelization.

5.4. Choice of dimension M 0 for noisy data
It is difficult to avoid imaging noise in real-world dataand it must therefore always be considered. We empiricallyobserved that M 0 is related to our method’s robustness towards noise. Thus, we determine an optimal M 0 by a validation on the noisy MERL sphere dataset.

Table 5: Surface normal estimation with varying M 0 on the noisy MERL sphere dataset under three sets of lights, (a) 60 lights,(b) 80 lights, and (c) DiLiGenT’s 96 lights. Values indicate mean angular errors in degrees.
(a)

(b)

(c)

M0

0/0

5/30

µ/λ
30/5

30/30

M0

0/0

5/30

µ/λ
30/5

30/30

M0

0/0

5/30

µ/λ
30/5

30/30

3
5
10
20

2.21
1.83
1.71
1.75

3.99
4.10
4.29
4.95

3.65
3.52
3.53
3.84

4.37
4.50
4.71
5.33

3
5
10
20

2.18
1.76
1.62
1.62

3.99
4.12
4.30
4.68

3.83
3.71
3.72
3.88

4.43
4.57
4.77
5.12

3
5
10
20

2.59
2.38
2.33
2.37

5.60
6.19
7.08
8.07

5.15
5.20
5.39
5.70

6.22
6.75
7.62
8.65

Setup: We apply 5-fold cross-validation on the MERL
BRDFs; the 100 MERL BRDFs are divided into 80
BRDFs for the BRDF tensor and 20 BRDFs for testing. We test varying M 0 = {3, 5, 10, 20} and varying noiseµ/λ = {0/0, 5/30, 30/5, 30/30} under three sets of lights:
60 lights, 80 lights, and DiLiGenT’s 96 lights.
Results: Table 5 shows mean angular errors of estimatedsurface normals in degrees for varying M 0 . In the noiselesscases (µ = λ = 0), a larger M 0 performs better. Contrarily,
in most of the noisy cases, M 0 = 3 produces lower angularerrors, indicating that M 0 = 3 is most robust to noise. For
this reason, we apply M 0 = 3 in all following experiments.

5.5. Discretized light directions
In all experiments so far, we assumed that the BRDF
tensor contains the known light directions. In practice, the
BRDF tensor rarely contains all of the known light directions and we should use pre-defined light directions closest
to the known light directions instead. Here, we examine
how the surface normal estimation accuracy is affected by
the discretization of light directions.
Setup: As pre-defined light directions in the BRDF tensor
we used 20001 discretized directions from a 0.5◦ equiangular sampling over the hemisphere. When a set of known
light directions is given, as illustrated in Fig. 2, we can sliceout a sampled BRDF matrix for a hypothesized surface normal and the set of light directions that are closest to the
known light direction in terms of cosine distance. We canthen follow the same estimation process used so far. We
performed this experiments with 20, 40, 60, 80, 100 known
lights. We used the MERL sphere dataset and performed
5-fold cross-validation over the BRDFs.
Results: Table 6 shows surface normal estimation resultsfor exactly known and for discretized light directions in the
BRDF tensor. We observe that there are only slight differences between both cases (< 0.1◦ in the mean angularerrors), which suggests that it is sufficient to precompute a
BRDF tensor for sufficiently finely discretized light directions and there is no need to recalculate the BRDF tensor

Table 6: The effect of discretizing the light directions andpicking the directions closest to a set of known directions.
The values indicate mean angular error in degrees.

Using known lights
Using discretized lights

100

number of lights L0
80
60
40

20

2.15
2.16

2.18
2.19

2.69
2.73

2.21
2.22

2.31
2.33

with known light directions each time one works with a newlight setup.

6. Discussion
In this paper, we have presented a photometric stereo
method based on discrete hypothesis-and-test search. The
proposed method can work with a diverse set of BRDFs thatare represented in a BRDF tensor and can determine surface
normals of a scene with spatially varying general BRDFs.
By putting most of the computation into a precomputation
step, we enabled a full search over all surface normal candidates, leading to a solution guaranteed to be optimal within
the bounds of the objective function and the discretization.
The approach is also supported by the fact that with the continuing increase of computation power, memory size, andthe availability of many-core processors, the applicabilityof the full search strategy is expanding. We are interested inseeing more applications along the direction.
One interesting question to investigate in the futurewould be how many BRDFs are needed to fully represent
real-world materials. Our BRDF tensor can hold an arbitrary number of BRDFs; therefore, we are interested in increasing the number of BRDFs in the tensor to study thisquestion, while it may only be empirically understood.

Acknowledgements
This work was supported by JSPS KAKENHI Grant
Number JP19H01123. Michael Waechter was supported
through a postdoctoral fellowship by the Japan Society for
the Promotion of Science (JP17F17350). Kiriakos N. Kutulakos was supported by the Natural Sciences and Engineering Research Council of Canada under the RGPIN program.

References[1] Neil Alldrin, Todd Zickler, and David Kriegman. Photometric stereo with non-parametric and spatially-varying
reflectance. In Computer Vision and Pattern Recognition(CVPR), 2008. 1[2] Brent Burley. Physically-based shading at Disney. In SIGGRAPH, 2012. 2[3] Guanying Chen, Kai Han, and Kwan-Yee Kenneth Wong.
PS-FCN: A flexible learning framework for photometricstereo.
In European Conference on Computer Vision(ECCV), 2018. 2, 7
[4] Lixiong Chen, Yinqiang Zheng, Boxin Shi, Art Subpa-Asa,and Imari Sato. A microfacet-based reflectance model forphotometric stereo with highly specular surfaces. In Computer Vision and Pattern Recognition (CVPR), 2017. 1, 2
[5] Tongbo Chen, Michael Goesele, and Hans-Peter Seidel.
Mesostructure from specularity. In Computer Vision and Pattern Recognition (CVPR), 2006. 2[6] Hin-Shun Chung and Jiaya Jia. Efficient photometric stereo
on glossy surfaces with wide specular lobes. In Computer
Vision and Pattern Recognition (CVPR), 2008. 2[7] Athinodoros S. Georghiades. Incorporating the Torrance
and Sparrow model of reflectance in uncalibrated photometric stereo. In International Conference on Computer Vision(ICCV), 2003. 2[8] Radoslav Harman and Vladimı́r Lacko. On decompositionalalgorithms for uniform sampling from n-spheres and n-balls.
Journal of Multivariate Analysis, 101(10):2297–2304, 2010.
5[9] Aaron Hertzmann and Steven M. Seitz. Example-based photometric stereo: Shape reconstruction with general, varying
BRDFs. Transactions on Pattern Analysis and Machine Intelligence (PAMI), 27(8):1254–1264, 2005. 1, 2, 3
[10] Berthold KP Horn and Katsushi Ikeuchi. The mechanicalmanipulation of randomly oriented parts. Scientific American, 251(2):100–113, 1984. 2[11] Zhuo Hui and Aswin C. Sankaranarayanan. Shape andspatially-varying reflectance estimation from virtual exemplars. Transactions on Pattern Analysis and Machine Intelligence (PAMI), 39(10):2060–2073, 2017. 1, 2, 3, 5, 6, 7
[12] Satoshi Ikehata. CNN-PS: CNN-based photometric stereo
for general non-convex surfaces. In European Conference
on Computer Vision (ECCV), 2018. 2, 5, 6, 7
[13] Satoshi Ikehata and Kiyoharu Aizawa. Photometric stereo
using constrained bivariate regression for general isotropic
surfaces. In Computer Vision and Pattern Recognition(CVPR), 2014. 1, 2, 5, 6, 7
[14] Micah K. Johnson and Edward H. Adelson. Shape estimation in natural illumination. In Computer Vision and Pattern
Recognition (CVPR), 2011. 2[15] Feng Lu, Xiaowu Chen, Imari Sato, and Yoichi Sato.
SymPS: BRDF symmetry guided photometric stereo forshape and light source estimation. Transactions on PatternAnalysis and Machine Intelligence (PAMI), 40(1):221–234,
2018. 1[16] Wojciech Matusik, Hanspeter Pfister, Matt Brand, andLeonard McMillan. A data-driven reflectance model. Transactions on Graphics (TOG), 22(3):759–769, 2003. 2, 4, 5
[17] Netanel Ratner and Yoav Y. Schechner. Illumination multiplexing within fundamental limits. In Computer Vision and
Pattern Recognition (CVPR), 2007. 5[18] Edward B. Saff and Arno Kuijlaars. Distributing many points
on a sphere. The Mathematical Intelligencer, 19(1):5–11,
1997. 5[19] Hiroaki Santo, Masaki Samejima, Yusuke Sugano, Boxin
Shi, and Yasuyuki Matsushita. Deep photometric stereo network. In Workshops of the International Conference on Computer Vision (ICCVW), 2017. 2[20] Boxin Shi, Zhipeng Mo, Zhe Wu, Dinglong Duan, Sai-Kit
Yeung, and Ping Tan. A benchmark dataset and evaluation
for non-Lambertian and uncalibrated photometric stereo.
Transactions on Pattern Analysis and Machine Intelligence(PAMI), 41(2):271–284, 2019. 5[21] Boxin Shi, Ping Tan, Yasuyuki Matsushita, and Katsushi
Ikeuchi. Bi-polynomial modeling of low-frequency reflectances. Transactions on Pattern Analysis and MachineIntelligence (PAMI), 36(6):1078–1091, 2014. 1, 2, 7
[22] Tatsunori Taniai and Takanori Maehara. Neural inverse rendering for general reflectance photometric stereo. In International Conference on Machine Learning (ICML), 2018. 7[23] Silvia Tozza, Roberto Mecca, Marti Duocastella, and Alessio
Del Bue. Direct differential photometric stereo shape recovery of diffuse and specular surfaces. Journal of Mathematical Imaging and Vision, 56(1):57–76, 2016. 2[24] Olivia Wiles and Andrew Zisserman. SilNet: Singleand multi-view reconstruction by learning from silhouettes.
British Machine Vision Conference (BMVC), 2017. 2[25] Robert J. Woodham. Photometric method for determining
surface orientation from multiple images. Optical Engineering, 19(1):139–144, 1980. 1, 2
[26] Sai-Kit Yeung, Tai-Pang Wu, Chi-Keung Tang, Tony F.
Chan, and Stanley J. Osher. Normal estimation of a transparent object using a video. Transactions on Pattern Analysis
and Machine Intelligence (PAMI), 37(4):890–897, 2014. 2