Photometric Stereo using Constrained Bivariate Regression for General
Isotropic Surfaces
Satoshi Ikehata∗
The University of Tokyo, Japan

Kiyoharu Aizawa
The University of Tokyo, Japanikehata@hal.t.u-tokyo.ac.jp

aizawa@hal.t.u-tokyo.ac.jp

Abstract
This paper presents a photometric stereo method thatis purely pixelwise and handles general isotropic surfaces
in a stable manner. Following the recently proposed sumof-lobes representation of the isotropic reflectance function,
we constructed a constrained bivariate regression problem
where the regression function is approximated by smooth,bivariate Bernstein polynomials. The unknown normal vector was separated from the unknown reflectance function
by considering the inverse representation of the image formation process, and then we could accurately computethe unknown surface normals by solving a simple and efficient quadratic programming problem. Extensive evaluations that showed the state-of-the-art performance using
both synthetic and real-world images were performed.

1. Introduction
Photometric stereo estimates the surface normals of an
object from appearance variations under different lighting
conditions. Since Woodham [20] first introduced the photometric stereo for Lambertian scenes, the extension to realworld objects which exhibit diverse appearances beyond
Lambertian model has drawn significant interest.
Traditionally, certain parametric reflectance models are
assumed to inversely solve the photometric stereo problem.
One of the most popular classes assumes a basic Lambertianmodel but augmented with outlier detection for handling all
non-Lambertian regions of the scene [21, 11]. This strategy
is numerically stable and relatively robust to shadows andimage noises, but complex reflections such as rough specularities can be highly disruptive. In contrast to the first class,a second class of methods treats non-Lambertian reflections
as inliers using the nonlinear bidirectional reflectance distribution function (BRDF) [9, 16]. While these methods are
more capable of handling a wide variety of objects includ∗ This work was supported by the Grants-in-Aid for JSPS Fellows (No.
248615).

ing rough surfaces, they may suffer from numerical instabilities derived from the complex nonlinear optimization.
Instead of explicitly modeling the parametric form ofreflectance, the monotonicity property of reflectance has
been recently integrated in the photometric stereo problem [2, 10, 17]. Chandraker and Ramamoorthi [6] showthat an isotropic BRDF consists of a sum of lobes whose
contribution to the reflected intensity decreases monotonically as the surface normal deviates away from the direction
where the reflectance lobe is concentrated (i.e., referred toas a preferred direction). Following this remark, the surface
normal has been recovered utilizing the monotonicity of thereflectance function under the assumption that the numberof lobes is one and its preferred directions are known (e.g.,
the lighting direction in [10] and half vector in [2, 17]).
While effective, these methods are highly disruptive when
the assumption on the preferred direction is incorrect or
the reflectance function is composed of two or more lobes.
Furthermore, to our knowledge, simultaneous estimation ofboth azimuth and elevation angles has never been achievedby enforcing the monotonicity of a reflectance function thathas a preferred direction that is different from the lighting
vector (note that [2, 17] assume that the azimuth angle ofthe surface normal is known).
This paper presents a photometric stereo algorithm foraccurate estimation of surface normals of a general isotropic
scene by enforcing the monotonicity of a reflectance function with an unknown lobe number and preferred directions.
For this purpose, the bivariate reflectance model is developed in Section 2 where pixelwise appearances are wellapproximated by a bivariate monotonic (and therefore invertible) smooth function of the dot products between thesurface normal and the lighting direction, and between thelighting and viewing directions. We may then consider
the inverse representation of the image formation process,where the unknown normal vector is separated from the unknown monotonic inverse reflectance function. By parameterizing the latter using the Bernstein polynomials [13], we
obtain a set of constrained linear equations in both the surface normals and reflectance parameters, leading to a sim-ple, quadratic programming problem.
The proposed framework benefits from an efficient pixelwise optimization that is easily amenable to parallel processing and does not require typical smoothness constraints
for both object structure and reflectance, which could disrupt the recovery of fine details.

2. Photometric stereo using constrained bivariate regression
In this section, we formulate the photometric stereo as
a constrained bivariate regression problem. Henceforth werely on the following assumptions:(1) The relative position between the camera and the object
is fixed across all images.
(2) The object is illuminated by a point light source at infinity from varying and known directions.
(3) The camera view is orthographic, and the radiometricresponse function is linear.

2.1. Problem Statement
Diverse appearances of real-world objects can be encoded by a BRDF (ρ) that relates the observed intensity I at
a given point on the object to the associated surface normaln ∈ R3 , the incoming lighting direction l ∈ R3 , and theoutgoing viewing direction v ∈ R3 via
I = ρ(n, l, v) max (nT l, 0),(1)

T

where max (n l, 0) accounts for attached shadows. Thereis a problem for photometric stereo in recovering the surface normal n of a scene by inversely solving Eq. (1) from acollection of m observations under different lighting conditions. Note that except for uncalibrated photometric stereo
problems such as [8], l and v are usually known.
Recently Chandraker and Ramamoorthi [6] have presented a semiparametric model of an isotropic BRDF thatis represented as a sum of K different functions giving
ρ=KXρk (nT αk ).

(2)

k=1

Here ρk are (unknown) nonlinear functions, and αk (i.e.,||αk || = 1) are called preferred directions, along whichρk are concentrated. It is known that physically valid reflectance functions satisfy the following requirements.
(L1) Monotonicity: ρk 0 > 0.
(L2) Nonnegativity: ρk ≥ 0.
(L3) Passing thorough the origin: ρk (0) = 0.
Chandraker and Ramamoorthi [6] have shown that inversely solving Eq. (2) under known surface normals gives
good estimates of a wide variety of isotropic BRDFs without suffering from the curse of dimensionality. Unfortunately, however, solving Eq. (2) directly is prohibitively dif-ficult in the context of the photometric stereo problem because of the numerous unknown parameters, some of which
are coincident in the same term (i.e., n, αk , ρk ). Therefore, most of the conventional photometric stereo algorithms have assumed that the dominant preferred directionof the reflectance function is unique and known [10, 2, 17].
Instead of this approach, we only assume that the preferred direction (αk ) of each function (ρk ) is lying on theplane spanned by the lighting and viewing directions as
αk =pk l + qk v,||pk l + qk v||(3)

where pk and qk are nonnegative unknown values (i.e.,pk ≥ 0 and qk ≥ 0). The degree of freedom of αk isactually 1 because ||αk || = 1. 1 Then, this assumption
provides us following important result.
Theorem: Suppose there is no shadow at a surface
point (i.e., ∀i nT li , liT v ≥ 0 and Ii ≥ 0, where i is theindex of the light) and ρ(n, l, v) in Eq. (1) has the formof Eq. (2), whose parameters satisfy the requirements ofa physically valid BRDF (L1)–(L3) and Eq. (3). Then,
it is guaranteed that there exists at least one continuous
bivariate function f (x, y) ∀x, y ∈ [0, 1], which satisfies
f ≥ 0, ∂f /∂x > 0, ∂f /∂y ≤ 0 and ∀i Ii = f (nT li , liT v).
Proof : From Eq. (3), nT αk is transformed into
pk nT l + qk nT v.
nT αk = p 2
pk + qk2 + 2pk qk lT v

(4)

Here we used ||l|| = ||v|| = 1. Eq. (4) illustrates thatnT αk is non-decreasing for nT l with fixed lT v andnon-increasing for lT v with fixed nT l since p, q are
non-negative constant values and nT v is constant overdifferent lightings. From (L1), it is guaranteed that eachρk (nT αk ) is also non-decreasing/non-increasing fornT l and lT v when either of them is fixed. Integratingthese results into Eq. (1) and Eq. (2), it is proved that
I is monotonic increasing for nT l with fixed lT v andnon-increasing for lT v with fixed nT l, which implies
we can always define continuous functions f (x, y)
which satisfy f ≥ 0, ∂f /∂x > 0,∂f /∂y ≤ 0
and ∀i Ii = f (nT li , liT v) since Ii ≥ 0 and∀i 0 ≤ nT li , liT v ≤ 1.
We illustrate this theorem in Fig. 1-(a). We note thatassuming that f (x, y) is always passing through the y-axis(i.e., f (0, y) = 0) does not limit any kind of isotropic
BRDF represented by Eq. (1) because I = nT lρ.
Following this theorem, we formulate the photometricstereo as a constrained bivariate regression problem whose
1 Note that this assumption does not violate most existing BRDF models
such as [7, 19, 15, 12, 3].

𝑥𝑖 , 𝑦𝑖 , 𝑧𝑖 = {𝒏𝑇 𝒍𝑖 , 𝒍𝑇𝑖 𝒗, 𝐼𝑖 }𝑧 = 𝑓(𝑥, 𝑦)

𝑥 = 𝑔(𝑦, 𝑧)

𝑧

1.0

𝑥

𝑦
(a)

0.0

𝑥

1.0

1.0

1.0

𝑧
(b)

0.0

𝑦

Figure 1. (a) Collections of {x, y, z} = {nT l, lT v, I} are lying
on a continuous function of z = f (x, y) which satisfies ∂f /∂x >0, ∂f /∂y ≤ 0 and f (0, y) = 0. (b) 3-d points lying on f are also
lying on a inverse function (x = g(y, z)).

goal is to recover an unknown surface normal n and a continuous bivariate function f from a collection of lighting
directions li and associated appearances Ii (i = 1, . . . , m),which satisfy the following equations and constraints:
Ii = f (nT li , liT v)

i = 1, . . . , m.

(5)(L4) Monotonicity (x): ∂f /∂x > 0.
(L5) Monotonicity (y): ∂f /∂y ≤ 0.
(L6) Nonnegativity: f ≥ 0.
(L7) Passing through y-axis: f (0, y) = 0.
We call Eq. (5) the forward bivariate reflectance model.
The major benefit of this formulation is that we do not needto explicitly approximate the number of lobes K and theirpreferred directions αk . However, there is a critical issue,
namely the coincidence of unknown parameters n and f
in the same term. We overcome this difficulty by a convenient, inverse representation of the imaging model applied
through a constrained bivariate regression framework.

2.2. Inverse Bivariate Reflectance Model
Strict monotonicity of f (x, y) (L4) guarantees theunique existence of the function giving x = g(y, f (x, y)) =g(y, z), which obeys the following requirements.
(L8) Monotonicity (y): ∂g/∂y ≥ 0.
(L9) Monotonicity (z): ∂g/∂z > 0.
(L10) Nonnegativity: g ≥ 0.
(L11) Passing through the y-axis: g(y, 0) = 0.
The proof, which has been omitted for brevity, is obvious
by seeing Fig. 1-(b). From the definition, each 3-D pointof {x, y, z} = {nT li , liT v, Ii }(i = 1, . . . , m) lying on f isalso lying on g as follows.
nT li = g(liT v, Ii )i = 1, . . . , m.

(6)

In contradiction to Eq. (5), we call Eq. (6) the inverse bivariate reflectance model. Our goal is now updated to recover
the surface normal n and a continuous bivariate function gwith some shape restrictions (L8)–(L11). The fundamental
advantage of Eq. (6) is that unknown variables of n and g
are separated, which contributes to simplifying the problem.
Constraints on g limit the solution space of Eq. (6), butthere are still multiple feasible solutions of a pair of n and g
because {liT v, Ii } are sparsely distributed on the valid rangeof {y, z}. To reduce the inherent ambiguity of the problem, we further assume a parametric model of the inverse
bivariate reflectance function g(y, z). Given that the lefthand side of Eq. (6) is linear in the unknown normal vectorn, for computational simplicity we would like to impose asimilar linearity on the right-hand side in our parameterized
representation of g(lT v, I) (we have omitted subscripts for
simplicity). For this purpose, we then choose to expressg(lT v, I) as a summation over p known nonlinear basis
functions gk (lT v, I) weighted by an unknown coefficient
vector β , [β1 , . . . , βp ]T , leading to the representationTg(l v, I) =p
X

βk gk (lT v, I).

(7)

k=1

While nonlinear in lT v and I, g(lT v, I) is clearly linear
in β. We need to choose gk carefully because estimating
a multivariate regression function subject to shape restrictions with compact support is challenging and usually very
time-consuming [18]. Here we adopt the bivariate Bernstein polynomials [13], where the shape-restricted regression function estimation is shown to be the solution of aquadratic programming problem [5, 18], making it computationally attractive. Furthermore, the Bernstein polynomials approximation naturally selects smooth functions withlittle computational effort, unlike other nonparametric regression functions (e.g., smoothing spline [4]), which implicitly enforces the smoothness of BRDF as in [2].
Bivariate Bernstein polynomials [13] are composed ofmultiple basis functions of the form
bk1 ,k2 (x1 , x2 , N1 , N2 ) = bk1 (x1 , N1 )bk2 (x2 , N2 ), 
Nibki (xi , Ni ) =xki (1 − xi )Ni −ki (i = 1, 2),ki(8)where 0 ≤ xi ≤ 1 and Ni is the order of the polynomial as
for xi , which will be chosen as a function of the sample sizem (e.g., Ni = o(mγi ) with γi > 0 suitably chosen via the
popular V-fold cross-validation method as shown in [18]).
We transform Ii (i = 1, . . . , m) to lie in the unit [0, 1] via
a simple linear equation as Ii ← Ii / max (I). Note thatliT v naturally lies in [0, 1] as we only consider the casev = [0, 0, 1]T and lz > 0. Then, the bivariate Bernsteinpolynomials approximation of g is represented as
x = g(y, z) = β T bNy ,Nz (y, z)
=

Ny Nz
X
X
ky =0 kz =0

βky ,kz bky ,kz (y, z, Ny , Nz ),(9)where bNy ,Nz , [b0,0 , . . . , bNy ,Nz ]T ∈ R(Ny +1)(Nz +1)×1
and β , [β0,0 , . . . , βNy ,Nz ]T ∈ R(Ny +1)(Nz +1)×1 .
Unlike the B-splines procedure (which may require
quadratic constraints on the coefficients) [4], shape restrictions (e.g., monotonicity, nonnegativity) on Eq. (9)are easily encoded via linear constraints, that is Aβ ≥ 0
and Cβ = 0, where A, C are shape restriction matrices.
Following [18], the shape restriction matrices required forour problem are defined as follows.
(1) Monotonicity: ∂g/∂y ≥ 0 and ∂g/∂z ≥ 0. (L8), (L9)The first-order partial derivatives of g with respect to yin Eq. (9) can be represented as∂g(y, z)/∂y
= Ny

Nz NX
y −1
X

(10)

where the coefficients of Bernstein polynomials (β) are restricted via the following equation:


Amono
Aβ =
β ≥ 0, Cβ = 0.
(12)
AnonnegCollecting variations of observations at the same pixel
under different lighting directions, Eq. (11) can be merged
into following linear problem:
LT n = B T β.

T
Here, B , [bNy ,Nz (l1T v, I1 ), . . . , bNy ,Nz (lmv, Im )] andL , [l1 , . . . , lm ]. By merging unknown variables (n, β),this problem is transformed into

(βky+1 ,kz − βky ,kz )bky ,kz (y, z, Ny − 1, Nz ).

P x = [L − B]T x = 0,kz =0 ky =0

From the definition in Eq. (8), it is easy to show thatall Bernstein basis polynomials are nonnegative withrespect to 0 ≤ y, z ≤ 1. Hence, the non-decreasing
constraint (i.e., ∂g/∂y ≥ 0) is simply achieved byenforcing βky+1 ,kz ≥ βky ,kz .
The non-decreasing
constraint with respect to z (i.e., ∂g/∂z ≥ 0) is
also achieved in the same manner.
The restriction
matrix for linear constraint Amono β ≥ 0 is represented as Amono = [ATy ATz ]T , which is composedof submatrices Ay ∈ RNy (Nz +1)×(Ny +1)(Nz +1) andAz ∈ RNz (Ny +1)×(Ny +1)(Nz +1) , where Ar ensures themonotonicity of the function with respect to r. 2 Note thatthe strict monotonicity constraint ∂g/∂z > 0 was eased to∂g/∂z ≥ 0 for computational simplicity.
(2) Nonnegativity: g ≥ 0. (L10)
The nonnegativity of g is guaranteed when ∀i βi ≥ 0.
Hence, the restriction matrix for a linear constraint
Anonneg β ≥ 0 is as Anonneg , diag([1, . . . , 1]) ∈R(Ny +1)(Nz +1)×(Ny +1)(Nz +1) .
(3) Passing through the y-axis: g(y, 0) = 0. (L11)
From the definition in Eq. (8), bky ,kz (y, 0) = 0 for all kz 6=
PNy
0. Therefore g(y, 0) =ky =0 βky ,0 bky ,0 (y, 0, Ny , Nz )becomes zero for all y when ∀ky βky ,0 = 0.
This constraint is encoded via a linear constraint
Cβ = 0, where C ∈ R(Ny +1)(Nz +1)×(Ny +1)(Nz +1) ,diag([1, 0, . . . , 1, 0, . . . , 1, . . .]) with Nz of 0 between 1.

2.3. Solution Method
By substituting Eq. (9) into the inverse bivariate reflectance model, Eq. (6) becomesnT li = β T bNy ,Nz (liT v, Ii )i = 1, . . . , m,

(11)

2 The concrete form of the shape restriction matrix is included in the
supplementary material.

(13)

(14)

where x , [nx , ny , nz , β0,0 , . . . , βNy ,Nz ]T and nx , ny , nz
are the three elements of the surface normal. Without lossof generality, we may avoid thePdegenerate x = 0 solution to Eq. (14) by constraining i xi = 1, which impliescT x = 1 where c = [1, . . . , 1]T .
Given the appearance variations (I1 , I2 , . . . , Im ) under
different known lighting conditions (l1 , l2 , . . . , lm ), the optimal surface normal (n) and model parameters (β) are recovered by solving the constrained linear problem:min ||P x||22 ,xs.t. Ãx ≥ 0 and C̃x = 0,

(15)


cT
. Note that Eq. (15) can0C
be effectively solved by the general quadratic programming.
where Ã , [0 A] and C̃ ,3. Handling Retroreflective Materials
If the reflectance of a target object obeys Eq. (2), ourmethod reasonably recovers the surface normal of the object by solving Eq. (15). However, one limitation of our
method is that this assumption is not satisfied in the presence of retroreflections that are often observed on roughsurfaces because Eq. (2) does not have the ability to represent this kind of reflection. In the presence of retroreflections, our surface normal estimation fails because of theviolation of (L5) by the behavior of retroreflections in thatthe power of reflections increases as lT v increases. While
it may limit available materials, fortunately we have foundthat our method practically handles retroreflective materials by simply reversing the direction of the monotonicity
constraint on lT v in Eq. (15) (i.e., use (L8)’ ∂g/∂y ≤ 0
instead of (L8)) because the retroreflections do not affectthe monotonicity for nT l while they unify the direction of
the monotonicity for lT v over the BRDF space. The problem of course is that we do not know whether the material is retroreflective or not. To overcome this difficulty,=

arg min
a

E

=

m
X

m
X

kn̂T li − aIi k22 ,

(16)

i=1

kn̂T li − ãIi k22 .

(17)

i=1

Here, n̂ is a recovered surface normal by solving Eq. (15)under the monotonicity constraint for lT v in either of two
directions. We simply choose the direction for which E issmaller. This strategy is very simple but very efficient as we
will show in Section 4.

4. Experimental Results
In this section, we evaluate our method on synthetic andreal image data. All experiments were performed on an Intel Core i7-2640M (2.80GHz, single thread) machine with
8GB RAM and were implemented in MATLAB. For thequantitative evaluation, we generated 32-bit HDR images ofa sphere (256 × 256) with foreground masks under different
BRDF settings: (A) common physical or phenomenological BRDF (Cook–Torrance [7], Ward [19] , Lafortune [12],Oren–Nayar [15] and Ashikhmin–Shirley [3]) and (B) the
measured MERL BRDF database [14]. Lighting directions
were randomly selected from a hemisphere with the object placed at the center. Additionally, for the third dataset,
denoted (C), we used real images for qualitatively evaluating our method in practical situations. For each dataset,shadows were removed via simple thresholding as in otherworks such as [16] (Thresholds for shadow removal were
fixed over algorithms but varied over objects manually). Because ground truth surface normals are provided in datasets(A) and (B), we quantitatively evaluated our method by the
angular error between the recovered normal map and theground truth when using these datasets.

4.1. Evaluation with Synthesized BRDF
We evaluated the performance using the synthesized images in dataset (A) generated under 100 different lightings using five common BRDFs3 . Here we compared3 Details of BRDF models are described in the supplementary material.

Mean Angular Errors (degree)

Cook-Torrance
Ward
Lafortune
Oren-Nayer
Ashikhmin-Shirley
Average

020 40 60 80 100 120 140 160 180 200 220 240

Ave. no. of nonshadowed pixels
Figure 3. Experimental results of dataset (A) with a varying number of images. For a fair evaluation, we display the average number of nonshadowed pixels that join our algorithm on the x-axis
instead of the number of images.
Mean Angular Errors (degree)

ã

10
9
8
7
6
5
4
3
2
1
0

8
Ours
Shi. et al [17] (0)Shi. et al [17] (0.05)
Shi. et al [17] (0.1)
Shi. et al [17] (0.15)
Shi. et al [17] (0.20)
Shi. et al [17] (0.25)
Shi. et al [17] (0.30)

7
6

5
4

3
2
1

0
CookTorrance

Ward

Lafortune

OrenNayer

AshikhminShirley

BRDF

Figure 4. Comparison with [17]. The values in parentheses are
the radian values added to the true azimuth angle map. The noisy
azimuth maps were used in [17]. The sign was chosen randomly.
0.04

Computation time/pixel (sec)

we present a practical approach for handling both non- andretroreflective materials. The important observation is thatwhen we incorrectly constrain the problem, the regression
usually fails. Therefore, after we have the regression outputs under both constraints, we can judge which constraint
was optimal by examining regression errors.
However, we have empirically found that comparing regression errors in Eq. (15) does not work since the flexible
Bernstein polynomials are generally well fitted to observations even though the constraint was not correct; instead, we
compute the following linear regression error E for choosing the optimal solution:

Ny=1 (m=100)
Ny=2 (m=100)
Ny=3 (m=100)
Ny=1 (m=300)
Ny=2 (m=300)
Ny=3 (m=300)

0.035
0.03
0.025
0.02
0.015
0.01
0.005
0
0

5

10
15
𝑁𝑦 + 1 × (𝑁𝑧 + 1)

20

25

Figure 5. Evaluation of computational time. We illustrate the perpixel computational time for each combination of the number ofBernstein basis functions and lightings.

our method with the standard Lambertian least-squaresregression-based approach [20] (LS) and a recent Lambertian sparse-regression-based approach implemented withthe sparse Bayesian learning [11] (SBL) (λ is fixed by

3
2.5
2

1

BQ
Ours

0.5
0
10

20

30

40

50

60

70

80

90 10010

Lafortune

7
6
5
4
3

2
1
0
10

20

30

40

50

60

20

30

40

50

60

70

80

Ashikhmin-Shirley

9
8
7
6
5
4
3
2
1
0
10

90 10020

30

40

Tlow (%)
Mean Angular Errors (degree)

Mean Angular Errors (degree)

Tlow (%)
8

Ward

70

80

90

100

4
3.5
3
2.5
2
1.5
1
0.5
0

Oren-Nayar

10

20

30

Tlow (%)

40

50

60

Tlow (%)

50

60

70

80

90 100

Tlow (%)

70

Mean Angular Errors (degree)

LS
SBL

1.5

4.5
4
3.5
3
2.5
2
1.5
1
0.5
0

Mean Angular Errors (degree)

Mean Angular Errors (degree)

Mean Angular Errors (degree)

Cook-Torrance

3.5

80

90 1005

Average over 5 BRDF

4
3
2
1
0
10

20

30

40

50

60

70

80

90 100

Tlow (%)

Figure 2. Experimental results of dataset (A) with five different BRDFs with varying frequencies in observations.

10−6 ). Our method was also compared with a recent parametric non-Lambertian photometric stereo method with thebiquadratic reflectance model [16] (BQ). In this experiment,
we fixed (Ny , Nz ) in Eq. (11) by (1, 5) to examine the robustness of our method against these parameters.
Evaluation with varying frequencies
We first evaluated our algorithm using observations withvarying frequency to present the flexibility of our modelin comparison with models used in previous works (e.g.,
BQ[16]) that only work for low-frequency reflectance values. Here high-frequency specularities were discarded byusing the nonshadowed pixels with intensities that wereranked below the Tlow % (10 ≤ Tlow ≤ 100).
The result is illustrated in Fig. 2. Overall, we observedthat our method outperformed other algorithms for almost
all frequencies for the Lafortune, Ashikhmin–Shirley, andOren–Nayar models, and performed competitively with
SBL for the Cook–Torrance and Ward models. Interestingly, because of our inverse reflectance model and the simple lT v constraint selection strategy, our method works wellfor the Oren–Nayar model which exhibits strong retroreflective reflections, and it also works well when all the frequencies are included unlike other methods4 . Unfortunately,
however, our method did not work when the number of images was very small, as will be discussed below.
Valid number of input images
We also evaluated our algorithm using a varying number ofimages to find the valid number required for effective recovery. The results are displayed in Fig. 3. We observed thatthe minimum number of images required to make the algo4 The result of BQ is similar as LS since the optimization procedure of
BQ is nonlinear and strongly affected by the initial estimation by LS.

rithm work was around 20 and more than 60 were required
for stable reconstruction because our method could sufferfrom the over-fitting when the number of the input images
was very small.
Comparison with another monotonicity-based approach
In this section, we compared our method with the recent
elevation-angle estimation algorithm [17] assuming that thedominant reflectance lobe was pointing at the half-vector
direction. We generated a ground truth azimuth angle
map and azimuth angle maps with some radian errors because [17] requires an azimuth angle map as input whereas
our method simultaneously recovers all elements in the normal. The radian errors varied from 0.05 to 0.3 were added
to the ground truth azimuth angle values, where the sign ofthe error was chosen randomly.
The results are illustrated in Fig. 4. We observed thatwhen true azimuth angles were given, [17] outperformed
our method in the Cook–Torrance and Ward datasets. However, as the amount of errors increased, the differences became smaller and finally our method outperformed [17].
As we expected, [17] did not work for the Lafortune andAshikhmin–Shirley models because these models violate
their assumption.
Evaluation of computational time
Here, we examined the computational time required for our
computation. We tried various combinations of Ny andNz and m in Eq. (9) and solved our optimization problem Eq. (11) using the lsqlin function in MATLAB.
The evaluation results are illustrated in Fig. 5. Here we
present a per-pixel computational time to solve one optimization problem. Therefore, the actual computational time
was twice that in the figure as we applied our algorithm

Mean Angular Errors (degree)

40

LS
BQ (Tlow=100)
Ours w/o retro-reflection detection

SBL
BQ (Tlow=25)
Ours w/ retro-reflection detection

35
30
25

20
15
10
0

yellow-matte-plastic
white-acrylic
white-diffuse-bball
violet-rubber
specular-orange-phenolic
white-marbledark-specular-fabricspecular-maroon-phenolic
specular-red-phenolic
ipswich-pine-221
red-specular-plastic
maroon-plastic
specular-yellow-phenolicgray-plastic
silver-paintpink-plastic
white-paintblue-rubber
pearl-paintred-fabricsilver-metallic-paint2
specular-violet-phenolic
alumina-oxide
gold-paint
purple-paintpink-jasper
yellow-plastic
yellow-phenolic
green-plastic
red-phenolicorange-paintyellow-paintteflonspecular-white-phenolic
red-plastic
light-red-paintgreen-acrylic
pink-fabricaventurninedelrinwhite-fabricdark-red-paintpvc
pink-felt
nylon
specular-green-phenolic
gold-metallic-paint2
polyethylene
violet-acrylic
pure-rubber
two-layer-gold
red-fabric2
gold-metallic-paint3
silver-metallic-paintneoprene-rubber
pink-fabric2
green-metallic-paintblue-acrylic
color-changing-paint1
alum-bronze
gold-metallic-paintchromecolor-changing-paint3
special-walnut-224
black-phenolic
two-layer-silver
specular-black-phenolicsilicon-nitrade
specular-blue-phenolic
red-metallic-painthematitenickelblue-fabricsteelcolor-changing-paint2
green-latex
beige-fabricblack-oxidized-steel
tungsten-carbide
chrome-steel
green-metallic-paint2
fruitwood-241
blue-metallic-paintblue-metallic-paint2
green-fabricdark-blue-paintgrease-covered-steel
aluminium
black-obsidian
polyurethane-foam
black-fabricblack-soft-plasticcolonial-maple-223
natural-209
ss440
cherry-235
brass
white-fabric2
pickled-oak-260
light-brown-fabric5

Figure 6. Comparison between different methods with dataset (B). Results are aligned in ascending order of the mean angular error of our
method with the retroreflection detection algorithm.

twice to distinguish retroreflective materials. We observedthat the computational complexity depends on the numberof basis functions rather than the number of lightings.

4.2. Evaluation with Measured BRDF
Here we evaluated the performance of our method tothe dataset (B). We generated images under 300 different
lightings for 100 different materials from the MERL BRDF
database [14]. In this experiment, our method was also
compared with LS [20], SBL [11], and BQ [16]. LS,
SBL, and our method used all the frequencies in observations (i.e., Tlow = 100) whereas only BQ was performed
with both Tlow = 25 and Tlow = 100 as that model was
originally designed to represent the low-frequency observations. In this experiment, we fixed Ny = 3 and Nz = 5
and performed our method with/without the retroreflection
detection algorithm described in Section 3 to verify the effectiveness of this process.
The results are illustrated in Fig. 6. We observed that
our method with our efficient retroreflection detection outperformed other algorithms for most of the materials. BQ(Tlow = 25) is more effective for some materials, but weemphasize that our method was capable of handling all frequencies in observations because of our flexible reflectancemodel, while BQ (Tlow = 100) does not work for mostmaterials. We also observed that the angular errors of
our method without retroreflection detection were relatively
larger than that of our method with retroreflection detection
for materials exhibiting strong retroreflections (e.g., MERL
fabrics), which indicated that our retroreflection detection
algorithm worked quite well for those materials. Finally, theaverage angular errors over 100 materials were 12.5 (LS),
6.2 (SBL), 13.1 (BQ, Tlow = 100), 1.7 (BQ, Tlow = 25),
2.4 (our method without retroreflection detection) and 1.2(our method with retroreflection detection), respectively.

4.3. Qualitative Evaluation with Real Images
We also evaluated our algorithm using real images: (1)100 images of two-face (this dataset is from [21]), (2) 100
images of doraemon, and (3) 44 images of fatguy (these twodatasets are from [11]). In this experiment, we compared
our algorithm with LS [20], SBL [11] (λ was fixed by 10−1 )and BQ [16] (for Tlow = 25 and 100). Note that our method
used a retro-reflective detection scheme and we fixed (Ny ,Nz ) by (1, 5). The threshold of the shadow removal for eachdataset was chosen manually but fixed over algorithms.
The experimental results are illustrated in Fig. 7. Here
we show both recovered surface normal maps and surface
meshes reconstructed by a poisson solver [1] with a fixedscale. We observe that our method succeeded to estimate
smoother and more reasonable normal maps and surface
meshes. We also observe that BQ (Tlow = 25) worked
poor for those datasets since shadows could not be completely removed by a simple thresholding therefore the lowfrequency component in the observation was not reliable.
In contrast to that, our method performed well since ourmethod could account all observations without discarding
the informative high-frequency component.

5. Conclusion
In this paper, we have proposed the constrained bivariate regression-based photometric stereo, which worked for

(a)

(b)

Input

(c)
LS
(𝑇𝑙𝑜𝑤 = 100)

(d)

(e)

(f)

SBL
(𝑇𝑙𝑜𝑤 = 100)(g)
BQ
(𝑇𝑙𝑜𝑤 = 100)

(h)

(i)
BQ
(𝑇𝑙𝑜𝑤 = 25)

(j)

(k)

Ours
(w/ retro-reflection detection)

Figure 7. Experimental results using real data (two-face, doraemon and fatguy). We illustrate (a) example images of target objects, andnormal maps recovered by (b) LS, (d) SBL, (f) BQ (Tlow = 100), (h) BQ (Tlow = 25) and (j) Ours (with a retro-reflective detection). We
also show surface meshes generated from normal maps in (c), (e), (g), (i) and (k).

various kinds of isotropic surfaces by exploiting various
conditions shared among physically valid BRDFs. Our detailed experimental results have shown the state-of-the-art
performance of our method for both synthetic and real data.
The current limitation is that the proposed method enforces
a global monotonic constraint with regard to lT v, whichmight not be true for surfaces with both the retro and specular reflections though they are rarely observed in the realworld. Another limitation is that we assumed that shadows
were discarded from images in advance, and this might be
impractical in real scenes. To ease this condition, we areinterested in incorporating the data cleansing scheme in asimilar manner as [11].

References[1] A. Agrawal, R. Raskar, and R. Chellappa. What is the rangeof surface reconstructions from a gradient field ? In Proc.
ECCV, 2006. 7[2] N. Alldrin, T. Zickler, and D. Kriegman. Photometric stereo
with non-parametric and spatially-varying reflectance. In
Proc. CVPR, 2008. 1, 2, 3
[3] M. Ashikhmin and P. Shirley. An an isotropic phong brdf
model. Journal on Graphics Tools, 5(2):25–32, 2000. 2, 5
[4] K. Bollaerts, P. H. Eilers, and I. Mechelen. Simple and multiple p-splines regression with shape constraints. British J.
Math. Statist. Psych., 59(2):451–469, 2006. 3, 4
[5] P. Chak. Semi-nonparametric estimation with bernstein polynomials. Economics Letters, 89(2):153–156, 2005. 3[6] M. Chandraker and R. Ramamoorthi. What an image reveals
about material reflectance. In Proc. ICCV, 2011. 1, 2
[7] R. Cook and K. Torrance. A reflectance model for computer
graphics. ACM Trans. on Graph., 15(4):307–316, 1981. 2, 5
[8] P. Favaro and T. Papadhimitri. A closed-form solution touncalibrated photometric stereo via diffuse maxima. In Proc.
CVPR, 2012. 2[9] D. Goldman, B. Curless, A. Hertzmann, and S. Seitz. Shapeand spatially-varying brdfs from photometric stereo. In Proc.
ICCV, October 2005. 1[10] T. Higo, Y. Matsushita, and K. Ikeuchi. Consensus photometric stereo. In Proc. CVPR, 2010. 1, 2
[11] S. Ikehata, D. Wipf, Y. Matsushita, and K. Aizawa. Robust
photometric stereo using sparse regression. In Proc. CVPR,
2012. 1, 5, 7, 8
[12] E. Lafortune, S.-C. Foo, K. Torrance, and D. Greenberg.
Non-linear approximation of reflectance functions. In Proc.
ACM SIGGRAPH, 1997. 2, 5
[13] G. G. Lorentz. Bernstein Polynomials. Chelsea Publishing
Company, New York, 1986. 1, 3
[14] W. Matusik, H. Pfister, M. Brand, and L. McMillan. A datadriven reflectance model. ACM Trans. on Graph., 22(3):759–
769, 2003. 5, 7
[15] M. Oren and S. Nayar. Generalization of lambert’s reflectance model. In In Proc. of the 21st annual conference
on Computer graphics and interactive tecniques, 1994. 2, 5
[16] B. Shi, P. Tan, Y. Matsushita, and K. Ikeuchi. A biquadratic
reflectance model for radiometric image analysis. In Proc.
CVPR, 2012. 1, 5, 6, 7
[17] B. Shi, P. Tan, Y. Matsushita, and K. Ikeuchi. Elevation angle
from reflectance monotonicity. In Proc. ECCV, 2012. 1, 2,
5, 6[18] J. Wang. Shape restricted nonparametric regression with
Bernstein polynomials. PhD thesis, North Carolina State
University, 2012. 3, 4[19] G. Ward. Measuring and modeling anisotropic reflection.
Computer Graphics, 26(2):265–272, 1992. 2, 5[20] P. Woodham. Photometric method for determining surface
orientation from multiple images. Opt. Engg, 19(1):139–
144, 1980. 1, 5, 7
[21] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and
Y. Ma. Robust photometric stereo via low-rank matrix completion and recovery. In Proc. ACCV, 2010. 1, 7