IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

1

Patch-based Uncalibrated Photometric Stereo
under Natural Illumination
Heng Guo† , Zhipeng Mo† , Boxin Shi∗ , Senior Member, IEEE, Feng Lu, Member, IEEE, Sai-Kit Yeung,
Member, IEEE, Ping Tan, Senior Member, IEEE, and Yasuyuki Matsushita, Senior Member, IEEE
Abstract—This paper presents a photometric stereo method that works with unknown natural illumination without any calibration
objects or initial guess of the target shape. To solve this challenging problem, we propose the use of an equivalent directional lighting
model for small surface patches consisting of slowly varying normals, and solve each patch up to an arbitrary orthogonal ambiguity. We
further build the patch connections by extracting consistent surface normal pairs via spatial overlaps among patches and intensity
profiles. Guided by these connections, the local ambiguities are unified to a global orthogonal one through Markov Random Field
optimization and rotation averaging. After applying the integrability constraint, our solution contains only a binary ambiguity, which
could be easily removed. Experiments using both synthetic and real-world datasets show our method provides even comparable
results to calibrated methods.
Index Terms—Uncalibrated photometric stereo, natural lighting, patch-based method, rotation averaging, intensity profile.

F

1

I NTRODUCTION

I ven an image sequence of a Lambertian object illuminatedby three non-coplanar directional lights, surface normalsof the object could be estimated by photometric stereo [56].
The pixel-level details of surface normal estimates are of greatinterest for applications in 3D computer vision such as visualinspection [17] and augmented reality [11].
The classic photometric stereo setup has two assumptions
on lighting – directional and calibrated lighting – restricting theapplicability of conventional photometric stereo. The directional
lighting model assumes a point light source placed far awayfrom the target object, and typically requires the data capture
to be conducted in a dark lab setting. The calibrated lighting
assumption needs an external step for measuring both lighting
intensities and directions, and calibrating lighting itself is also
an ongoing research problem [44]. If the former assumption isrelaxed, the problem becomes calibrated photometric stereo undernatural illumination, while relaxing the latter assumption leads touncalibrated photometric stereo under directional lighting. A fully
calibration-free method under general lighting is desired because
it will push photometric stereo from the laboratory setup to thepractical wild environment and simplify the effort of 3D scanningfor non-experts.

G

•
•
•

•
•
•
•

H. Guo and Y. Matsushita are with the Department of Multimedia Engineering, Graduate School of Information Science and Technology, Osaka
University, Japan. E-mail: {heng.guo, yasumat}@ist.osaka-u.ac.jp.
Z. Mo and P. Tan are with the School of Computing Science, Simon Fraser
University, Canada. E-mail: {zhipeng mo, pingtan}@sfu.ca
B. Shi is with the National Engineering Laboratory for Video Technology,
Department of Computer Science and Technology; Institute for Artificial
Intelligence, Peking University, China and Beijing Academy of Artificial
Intelligence, China. E-mail: shiboxin@pku.edu.cn.
F. Lu is with the State Key Laboratory of Virtual Reality Technology
and Systems, School of Computer Science and Engineering, Beihang
University, China. E-mail: lufeng@buaa.edu.cn.
S.-K. Yeung is with the Division of Integrative Systems and Design and theDepartment of Computer Science and Engineering, Hong Kong Universityof Science and Technology, Hong Kong. E-mail: saikit@ust.hk† H. Guo and Z. Mo contributed equally to this work.
∗ B. Shi is the corresponding author.

Natural lighting

Directional lighting

Fig. 1: Natural lighting vs. directional lighting.
However, generalizing calibrated and directional lighting assumptions at the same time will make the problem rather complicated. As shown in Fig. 1, scene points with unique surface normals are illuminated by different lighting directions and intensities
under natural lighting, whereas in the directional lighting case, thewhole object is lighted by a single lighting direction. Besides, withuncalibrated directional lighting, a 3 × 3 linear ambiguity [51]exists with the estimated surface normal fields after factorizing
the image observations, but the ambiguity in uncalibrated naturallighting case will further be extended to a higher-dimensional
linear ambiguity and cannot be fully removed [7].
To solve photometric stereo with uncalibrated natural illumination, existing methods [3], [28], [41], [47] require a roughshape of the target object. Although the initial shape can beobtained from multiview geometry [47], object shape prior (e.g.,face [28]) or RGBD cameras [19], it either needs extra systemsetup or restricts the application into pre-defined shapes. Recent
work [20] use a balloon-like perspective depth map for the shapeinitialization and estimate the surface shape and reflectance byan end-to-end variational optimization framework. However, therecovered shape accuracy after the optimization is still sensitive todepth initialization. On the other hand, Jung et al. [27] restrict
the natural illumination as the skylight and solve the outdoor
photometric stereo based on the prior of skylight distribution.
Brahimi et al. [10] provide a closed-form solution for photometricstereo under general unknown lighting and perspective camera
projection. However, the environment lighting is approximatedby global first-order spherical harmonics, which has a gap with

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

2

Normal estimation in local patch

Graph-based patch merging method (GPM)

Env. map𝒍ҧ𝑘Graph 𝒢𝑏
Consistent surface
normal set 𝒮𝑛𝑇෡ 1𝑇 𝑶2 𝑵෡ 𝑇2
𝑶1 𝑵...𝑶𝑝 𝑵෡ 𝑝

෡𝑘
𝑵

𝑑𝑖

𝑶𝑖

Surfacenormalclustering

SVD,
Constant
albedo

𝒩𝑖

𝑶𝑖

𝒑
𝑶𝑗

𝒩𝑗

Build
graph

෡ 𝑖,𝑗 = 𝑶𝑗𝑇 𝑶𝑖
𝑶Equivalent directional
lighting

Local ambiguity

෡ 𝑖,𝑗
𝑶

𝑑𝑗

𝑑መ𝑗,𝑘

MRF
Optimization

Global ambiguity

𝑑𝑘

𝑶𝑗
෡ 𝑖,𝑘
𝑶
෡𝑗,𝑘
𝑶

Integrability
constraint

𝑶𝑘

Graph 𝒢𝑜

Relative orthogonaltransformation set 𝒮𝑜

Input images𝑑መ 𝑖,𝑗𝑑መ 𝑖,𝑘ෝ𝒋 𝒑 ><ෝ𝒏𝒊 𝒑 , 𝒏𝒩𝑘
Patchextraction

෡
𝑵

𝑹𝑖

෡ 𝑖,𝑗
𝑹
෡ 𝑖,𝑘
𝑹

𝑹𝑗

෡𝑗,𝑘
𝑹

Rotationaveraging𝑵

𝑹𝑘

Graph 𝒢𝑟
OutputFig. 2: Complete pipeline of patch-based uncalibrated photometric stereo. Variables shown in red represent the unknowns in our method.
real-world natural illumination [8]. In summary, existing methods
for uncalibrated photometric stereo under natural illumination are
still limited due to the requirement of initial shapes and restrictive
lighting approximation models.
In this paper, we propose a photometric stereo method for
uncalibrated natural illumination that relieves the requirements
used in existing methods. We develop a “divide and conquer”
approach to first “divide” the problem into tractable sub-problems
with locally-resolvable ambiguity, and then “conquer” them jointly
by merging all sub-results as a complete solution. Our key observation is that for a small surface patch with slowly varying
normals, the visible hemisphere of an environment map also showssmooth changes. Therefore, the environment lighting for that local
patch could be approximated as equivalent directional lighting bysumming up all samples on the visible hemisphere. Based on theabove observations, we present our method as shown in Fig. 2,
where the “divide” and “conquer” processes are corresponding to
the two modules in the pipeline.
In the “normal estimation in local patch” stage, we assume thesurface normals in an extracted patch have similar directions sothat the patch illumination can be approximated by equivalent directional lights. Following directional lighting-based uncalibrated
photometric stereo techniques [22], we further assume the local
patch has uniform albedo and non-planar shape, then the surface
normals can be recovered up to an orthogonal ambiguity. In the“graph-based patch merging” stage, we resolve the orthogonal ambiguity in each patch and merge local shapes to a complete surface.
Specifically, we first cluster consistent (equal) surface normal pairs
and use them to calculate relative orthogonal transformations,which describe the geometry relationship among patches. Thenan orthogonal ambiguity graph Go is constructed with nodes andedges being set to the unknown patch-wise orthogonal ambiguities
and the known relative orthogonal transformations. As the 3 × 3
orthogonal ambiguity can be decomposed into a binary part and a
3D rotation part, we divide the orthogonal ambiguity graph into abinary ambiguity graph Gb and a rotation ambiguity graph Gr . We
formulate the binary ambiguity estimation on Gb as a per-patchlabeling problem, and solve it by a Markov Random Field (MRF)
optimization [31]. Guided by the relative rotations (edges of therotation ambiguity graph Gr ), we solve the rotation ambiguities
in each patch by introducing rotation averaging [21] algorithms
which has been widely applied in structure from motion framework [55]. After the patch merging stage, the unknown patchwise orthogonal ambiguities can be determined up to a global
orthogonal ambiguity. By further assuming the whole surface to be

integrable, this global orthogonal ambiguity can be finally reduced
to a convex/concave ambiguity.
An earlier version of this work appeared in [35]. Different
from the graph-based patch merging method presented in thispaper, the patch merging process in [35] takes consistent surface
normal pairs as constraints and constructs an angular distancematrix with the element calculated by propagating angular distance along the shortest path between any two surface normal
directions. Then the complete surface normal map is solved up toa global orthogonal ambiguity by conducting matrix factorizationon this angular distance matrix. We refer this method as matrixbased patch merging method (MPM) corresponding to our newly
proposed graph-based patch merging method (GPM). Comparedwith [35], this work improves the surface normal estimation
accuracy by replacing MPM with GPM, and provides analysis ofsurface normal clustering under natural illumination via consistentorthogonality condition. To demonstrate the effectiveness of our
new method, additional experiments on both synthetic and real
data are also presented. To summarize, the main contributions of
our work are as follows:
1)

2)

3)

We explore the equivalent directional lighting model tosolve patch-wise surface normal up to local ambiguities,
bypassing the explicit requirement of global informationof environment maps.
We extend the surface normal clustering via intensity profiles from directional lighting to general lighting case, andpropose a consistent orthogonality condition to extract
consistent surface normal pairs.
We introduce rotation averaging and MRF optimization tosolve patch-wise orthogonal ambiguities and merge local
surface normal solutions to a complete surface normal
map up to a global orthogonal ambiguity.

Our output surface normal map only contains a concave/convex binary ambiguity. As proved in [10], it is an inherent
ambiguity in uncalibrated photometric stereo under orthogonal
camera projection and cannot be solved with image cues only.
However, it can be either manually removed with little effort or
resolved with shape prior. Together with our previous version [35],the proposed equivalent lighting model and the “divide and conquer” framework is the first strategy solving photometric stereo
up to minimum inherent ambiguities under natural illumination
without relying on shape priors.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

2

R ELATED W ORK

There are two major restricting assumptions that need to be relaxed
for photometric stereo [56] to be applied to practical applications – calibrated directional lighting assumption and Lambertianreflectance assumption. Correspondingly, to make photometricstereo work in more realistic scenes, there are two directions togeneralize the conventional approach – generalization of lighting
assumption and generalization of the reflectance model. This
paper focuses on the former problem, thus both calibrated anduncalibrated photometric stereo methods with non-Lambertian
objects (e.g., [5], [37], [32], [49], [13]) are beyond the scope,
and we refer the readers to [50] for a comprehensive review andcomparison of non-Lambertian photometric stereo methods.

3

TABLE 1: Summary of uncalibrated photometric stereo methods
under natural illumination, where f , o, p, and k represent the
numbers of images, spherical harmonic (SH) lighting basis, valid
pixels, and extracted patches, respectively. Our method solves uncalibrated photometric stereo under a moderately flexible lighting
model without requiring a initial shape prior.
MethodInitial shape

Lighting
modelLighting
parameter

Representation
power

[10]
[19], [41][20][47], [46]
[43]
[33]
Ours

None
Depth sensor
Visual hull [38]MVGa
Planarc
Depth sensor
None

Global SH
Global SH
Global SH
Global SH
SV-SHbSV-SH
SV-directional

f ×4
f ×o
f ×o
f ×o
f ×3×pf ×o×kf ×3×k

Weak
Weak
Weak
Weak
StrongStrongModeratea2.1

Calibrated, directional lightingbc

The calibrated Lambertian photometric stereo with directional
lighting assumption is the most classic setup. The first photometric stereo work [56] and its robust extensions rely onthese assumptions. Various robust approaches have been proposed
to eliminate deviations from the classic model by treating thecorrupted measurements as outliers, such as Random Sample
Consensus (RANSAC) [36], [52], median-based approach [34],low-rank matrix factorization (Robust-PCA) [57], and expectation
maximization [58]. Wu et al. [59] formulate the calibrated dense
photometric stereo problem as a Markov network. The per-pixel
surface normal, encoded in the graph node, is optimized byminimizing the surface geometry smoothness (smoothness term)
and the distance between its normal initialization (data term).
Similar to their method, our proposed GPM also formulates thepatch merging problem to a graph structure, but with different
graph content and optimization scheme. As we have no lightcalibration or initial normal map, our optimization target encodedin the node is the patch-wise orthogonal ambiguity constrainedby the relative orthogonal transformations assigned to the graphedges. The optimizations for the binary ambiguity and the rotation
ambiguity are solved with MRF and rotation averaging, separately.
2.2

Calibrated, natural lighting

Natural illumination can be calibrated directly by using a mirror
sphere as a light probe or indirectly by approximating sunlight asa dominant directional source. With mirror sphere measured environment maps, Yu et al. [62] show photometric stereo results bydirectly sampling the captured natural illumination. Ackermann etal. [2] implement photometric stereo for outdoor webcams using atime-lapse video, and Abrams et al. [1] show the necessity of using
images taken over many months (thousands of images) for sufficiently observing illumination variations. Jung et al. [26] develop
parameterized sun and sky lighting models to apply photometricstereo under outdoor illumination. Their latter work [27] refines
the sky model and obtains better normal estimates on cloudy
days. Shen et al. [45] provide an analysis about the limitation ofpoint light source modeling for 1-day outdoor photometric stereo.
Hold-Geoffroy et al. [25] show that outdoor observations recorded
within a few hours could constrain a reliable normal estimation.
2.3

Uncalibrated, directional lighting

Photometric stereo without calibrated lighting as known input
is called uncalibrated photometric stereo. Even if the lighting
assumption is directional lighting, the solutions to both surface

Multiview geometry
Spatially-varying (SV) spherical harmonic (SH) lighting
The method [43] is validated by near-planar objects.

normal and lighting are not unique due to some inherent ambiguities. The shape (or lighting) can be estimated up to a 3 × 3 linear
ambiguity [22]. When the surface is integrable, this ambiguity
further reduces to a 3-parameter Generalized Bas-Relief (GBR)
ambiguity under orthographic projection [9], [63] and vanishes
under perspective camera projection [39]. Existing methods focus
on the estimation of the 3 unknowns in GBR ambiguity to recover
the normal estimates by using priors on albedo [6], [48], detecting
local maximum diffuse points [40], or reflectance symmetry [54],[60]. If multiview inputs are available, the directional lighting
directions could also be indirectly estimated, and photometricconstraints are used to refine the shape [23], [24]. The lighting
can also be semi-calibrated with directions being provided andintensities remaining unknown [14].
2.4Uncalibrated, natural lightingThis is the most challenging category of lighting conditions since
it is general and unknown. Table 1 summarizes existing uncalibrated photometric stereo methods under natural illumination.
Brahimi et al. [10] approximate the shading of the whole surface
with a first-order SH globally, where the number of lighting
parameters to be optimized is 4f . With the integrability constraint,
they show that the uncalibrated natural light photometric stereo
under perspective camera projection is well-posed. However, thefirst-order SH is a simplified natural lighting representation. For
the second-order SH representation (o = 9), there is a 9×3 (= 27
unknowns) linear ambiguity in estimated surface normals [7].
Unfortunately, this high-dimensional ambiguity cannot be completely removed without additional information. Existing methods
require initial shapes from depth sensor [19], [41], or multiview
geometry [4] to make the problem solvable. Recently, Haefner et
al. [20] propose a variational optimization framework to recover
shape, reflectance, and illumination jointly. Although their method
automatically initializes shape from silhouette [38], the embedded
non-convex optimization framework is still sensitive to the initialization of depth, albedo, and lighting vectors. Compared tomodeling the natural illumination for the whole surface with a
global SH lighting, existing methods [43], [33] propose spatiallyvarying spherical harmonic (SV-SH) lighting models. Maier et al.
[33] divide a surface into k patches and model the per-patch
SH lighting independently. Quéau et al. [43] further optimize
the per-pixel SV-lighting direction directly. Although the SV-SH
model can accurately represent the natural light, the uncalibrated

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

photometric stereo becomes highly ill-posed due to numerous unknown lighting parameters to be optimized. To make the problem
solvable, these methods require dedicated shape initialization [33],[43] and non-physical lighting regularization [43].
Similar to Maier et al. [33], our method models the natural illumination with SV-directional lighting in local patches. Comparedto the global SH lighting approximation, our lighting model has astronger representation power for real-world natural illumination.
Even with this flexibility, based on physical assumptions onpatches (uniform albedo, non-planar shape), we can directly obtaina patch shape up to an orthogonal ambiguity without requiring
shape initialization.

3

Env.
Maps

Patch
size

𝑣𝑣 = 1.71°

𝑣𝑣 = 0.25°
𝑟𝑟 = 1

𝑣𝑣 = 5.10°
𝑟𝑟 = 10

𝑟𝑟 = 30Beach

N ORMAL E STIMATION IN L OCAL PATCH

Our method is based on the Lambertian image formation model
under natural light. We ignore the cast shadow (self-occlusion)and assume the camera is radiometrically calibrated or has alinear response, i.e., the pixel brightness equals to the radianceof the scene. Let us consider a photometric stereo image sequence
illuminated by f different environment maps. In default, for eachvalid pixel, we extract a patch Nk (k = {1, 2, . . . , p}, where p is
the total number of pixels) centered at the pixel location.
In the following, we first approximate the illumination at a
local patch with an equivalent directional lighting model, and thenestimate the surface normals within the patch by conventional
uncalibrated photometric stereo algorithms.
3.1

4

Equivalent Directional Lighting Model

Given a scene point with Lambertian albedo ρ and surface normaln = [nx , ny , nz ]> ∈ S2 ⊂ R3 , its pixel brightness is written asZb=
ρL(ω) max((n> ω), 0)dω,(1)Ω
2where ω ∈ S ⊂ R3 is a unit vector in the visible hemisphere Ω,and L(ω) is the environment lighting intensity from direction ω .
For any surface normal vector nk , it uniformly receives illumination from direction ω sampled on the visible hemisphereΩk = {ω | n>k ω ≥ 0} of the environment map. Then for anyω ∈ Ωk we may perform the spherical integration over Ωk toobtain the pixel brightness:Zbk = ρn>
L(ω)ωdω = ρn>(2)kk l̄k ,
Ωkwhere l̄k denotes an equivalent directional lighting as the integralof all samples in Ωk , and the subscript k indicates that for different surface normals, they face different visible hemispheres andtherefore correspond to different equivalent directional lighting.
Note here n is a unit vector, but l̄ is not necessary of length one
since it encodes intensity scaled directional lighting direction.
We assume the surface normals in a small patch have similar
directions. In this way, the natural illumination does not show
abrupt changes for scene points within the patch. Given two
surface normals within patch Nk , we measure their angular difference as hnk,i , nk,j i = arccos(nk,i , nk,j ). To evaluate surface
normal’s variation in aP
patch, we define the mean patch angulardifference by vk = p1k i hnk,i , nk,c i, where c is the patch center
index and pk is the number of scene points in that patch. Then fora surface patch with small mean angular difference vk , all surface
normals should share approximately the same visible hemisphereΩk as well as l̄k , so their brightness could be modeled by a

Uffizi0

Grace

Fig. 3: Illustration of environment lighting approximation. Patches
highlighted by concentric red circles contain varying radii r andmean angular difference of surface normal v . For each patch, wedraw the equivalent directional lighting directions (dot on spheres)and intensities (red means strong while black means weak) for
three different environment maps (figures courtesy of [16]).

single directional light as illustrated in Fig. 2. Similar lighting
representation has been applied in [25], [43].
We illustrate and verify our lighting assumption using a
synthetic experiment. Given a surface normal, we calculate its
equivalent directional lighting by summing up all samplings on its
visible hemisphere of the environment map, and draw the intensityand direction of such a lighting vector on the sphere as shown in
Fig. 3. We use a sphere normal map of 256×256 pixels (the radiusof the sphere is 128 pixels in the image domain) and calculate
the equivalent lighting under three light probes from [16]. By
selecting central patches with the radius of {1, 10, 30} pixels(indicated as red circles), the mean angular difference of surface
normals increases from 0.25◦ to 5.10◦ , leading to more scatteredequivalent directional lighting distributions. For relatively smaller
patches (radius ≤ 10, around 300 pixels) whose surface normals
having smaller variation, the corresponding lighting vectors arehighly concentrated. In such case, it is safe for us to apply
directional lighting assumptions in a patch-wise manner. In thefollowing computation, we neither know the direction and intensity of equivalent directional lighting nor solve them explicitly,
while we develop an uncalibrated photometric stereo method tosolve the surface normal directly.

3.2 Uncalibrated Photometric Stereo based on Equivalent Directional Lighting
Assume a local surface patch Nk is illuminated by f different equivalent directional lighting Lk = [l̄k,1 , l̄k,2 , · · · , l̄k,f ].
Denote the matrix stacking all surface normal vectors n> inpatch Nk in a row-wise manner as Nk , denote the patch albedo

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE5

as ρk = [ρ1 , · · · , ρpk ], then the image brightness of this patch,denoted as Bk , could be written as follows

Bk|pk ×f = diag(ρk )Nk|pk ×3 Lk|3×f ,(4)where ñk,i ∈ R3 is the i-th (i = {1, 2, . . . , pk }) row vectorof Ñk . Without losing generality, we set αk = 1. As Yk =Q>
k Qk is an symmetric matrix, we can solve it if the local patch
contains at least 6 pixels with varying surface normals (Pleaserefer Appendix C for planar patch):>>tri(ñk,1 ñ>tri(Yk ) = 1,(5)k,1 ) · · · tri(ñk,pk ñk,pk ){z} | {z }|Eywhere tri(·) operator extracts the upper triangle matrix elements
as a vector. The residue kEy − 1k22 is recorded as eka to measure
the reliability of uniform albedo assumption. Wepconduct SVD on
Yk such that Yk = ŨΣ̃Ũ> and assign Q̂k as Σ̃Ũ> . Then we
obtain pseudo surface normal map as N̂k = Ñk Q̂k . It has been
proved in [9], [48] that the uniform albedo constraint reduces the3 × 3 linear ambiguity in Ñk to an orthogonal one in N̂k suchthat

Bk = N̂k O>
k Ok L̂k ,(6)where Ok ∈ O(3) is the orthogonal ambiguity that varying frompatch to patch. N̂k and L̂k are the pseudo surface normals andequivalent directional lighting up to an orthogonal ambiguity w.r.t.
their corresponding ground truth, i.e.,>N>k = Ok N̂k ,

Lk = Ok L̂k .

1

0

0

(3)

where diag(·) is a diagonalization operator and pk is the total
number of pixels in patch Nk . This representation is different fromspherical harmonics for natural light, where a high-dimensional
matrix decomposition (9D decomposition for a second order
spherical harmonics) exists with unknown lighting [7], [47].
According to Eq. (3), for each patch the equivalent directional lighting model relieves the problem to be the Lambertian
photometric stereo under unknown directional lighting, which
is a well-studied research area with tractable solutions. So we
perform SVD on Bk , as it was done in classic uncalibrated
photometric stereo methods [22]. The SVD decomposition gives
us Bk = UΣV> , wherein ideal case Σ only contains√ three
non-zerodiagonal
elements.
Wefurtherdenote
Ñ
=
U
Σ and
k
√
L̃k = ΣV> , where Ñk and L̃k are pseudo surface normals
and pseudo equivalent directional lighting for each patch. Here,both the normal and lighting solutions contain an unknown 3 × 3
linear ambiguity, denoted as Qk , since any invertible matrix canbe inserted between Ñk and L̃k to maintain the equality.
As we work on small patches, it is safe to assume a piecewise uniform albedo. Suppose pixels within the patch Nk havethe same albedo αk , i.e. ρk = αk 1, the pseudo surface normals
for this patch should satisfy>kñk,i Qk k22 = ñk,i Q>
k Qk ñk,i = αk ,

20

(7)

Non-uniform albedo across patches. As shown in Fig. 4(a),patches across the boundary of different albedos cannot keep
uniform albedo assumption. Also, the natural illumination within
patches near the feet and the ear part of the B UNNY object cannot
be treated as equivalent directional lighting since the surface
normals at these regions vary significantly. Therefore, patch-wise

(a)

(b)

(c)

Fig. 4: An example of non-uniform albedo causing large errors
across patches. (a) Image observation of the B UNNY object with
non-uniform albedo. (b) Mean angular errors (degree) of thepatch-wise pseudo surface normals w.r.t. the true surface normals.
Each pixel value encodes the mean angular error of the estimatedsurface normals for the patch centered at that pixel location. (c)Confidence map of patch surface normal estimation.

surface normal estimates in these regions are inaccurate, as visualized in Fig. 4(b). Here we first fit an orthogonal matrix to align thepatch-wise pseudo surface normals with the corresponding groundtruth. Then the mean angular error between aligned pseudo surface
normals and the truth surface normals are calculated to measure
the patch-wise surface normal estimation accuracy. Hereafter wedenote this error map as patch surface normal error map.
To reduce the influence of these inaccurate local surface
normal estimates in the following patch merging process, wedefine a confidence metric to measure the reliability of normal
estimation. For a surface patch Nk , we evaluate the equivalentdirectional lighting approximation by defining a normalized patch
re-rendering error ekr = kBk − N̂k L̂k k2F /kBk k2F and test theuniform albedo assumption by the residue eka calculated from
Eq. (5). Based on these two metrics eka and ekr , we define thesurface normal estimation confidence of patch Nk as follows,kkck = e−(βer +γea ) ,(8)where β and γ are the coefficients used to balance ekr and eka , andwe set them as 5 and 0.5 empirically. As shown in Fig. 4(b-c),
the confidence values of all patches are consistent with the patch
surface normal angular error map.

4

G RAPH - BASED PATCH M ERGING M ETHOD

For each patch, now we have estimated pseudo surface normal N̂kup to an orthogonal ambiguity Ok , with surface normal confidence
measured by ck . In this section, we will discuss how to merge
all the patches into an entire surface. Specifically, we first show
consistent surface normal pair extraction from patch overlapping
regions and intensity profiles, followed by the calculation ofrelative orthogonal transformations among patches. Taking relative
orthogonal transformations as constraints, we introduce MRF optimization and rotation averaging to solve the patch-wise orthogonal
ambiguities and merge the whole surface normal up to a global
orthogonal ambiguity. This global ambiguity is finally reduced toa concave/convex ambiguity by addressing integrability.
The MPM proposed in our early work [35] conducts thisstep by taking consistent surface normal pairs as constraints andcreating an angular distance matrix with its element filled bypropagating angular distance along the shortest path between every
surface normal pairs. The whole surface is then obtained by matrix
factorization on this angular distance matrix. Please refer to theoriginal paper in [35] for details. However, The MPM suffers fromIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

error accumulation during the propagation process. Also, the angular distance between two surface normals could be constrained by
all possible paths connecting the corresponding scene points, only
selecting the shortest path to constraint surface normals cannot
guarantee a globally optimized result. As discussed below, thenewly proposed GPM avoids the accumulative error in MPM [35]and optimizes all the patch connections simultaneously.

6

𝒏𝑘1 (𝐰)

𝒏𝑘1 (𝐩)𝒏𝑘2 (𝐩)

𝒩𝑘1 𝐰

𝒏𝑘3 (𝐬)𝒏𝑘2 (𝐪)𝐩𝚯

𝑳𝑘2 (𝐪)𝐬

𝑳𝑘3 (𝐬)

𝐪

𝒩𝑘3

𝒩𝑘2

(a)

4.1

Consistent Surface Normal Clustering

𝑵𝑘1 (Θ)As shown in Fig. 5(a), we provide three surface patches
Nk1 ∼ Nk3 covering scene points {w, p, q, s}. For any scene
point p located at the overlapping region Θ = Nk1 ∪ Nk2
shown in the highlight area, the true surface normals from different
patches at this point are consistent, i.e., hnk1 (p), nk2 (p)i = 0.
In Fig. 5(b), we show the relationship of surface normals in theoverlapping region between two patches. The unknown orthogonal
ambiguities in patch Nk1 and Nk2 are denoted as Ok1 and Ok2 .
Since the true surface normals in the overlapping region Θ areconsistent, i.e., Ok1,k2 = I, the pseudo surface normals of two
patches in this region can be aligned by>Ôk1,k2 = O>
k2 Ok1,k2 Ok1 = Ok2 Ok1 .

(9)

Obviously Ôk1,k2 ∈ O(3) encodes the relationship between
unknown orthogonal ambiguities of the two patches. We nameit relative orthogonal transformation and it can be solved byaligning pseudo surface normals in the overlapping regions, i.e.,>2Ô∗k1,k2 = argmin kÔk1,k2 N̂>
k1 (Θ) − N̂k2 (Θ)kF ,Ôk1,k2(10)
s.t. Ôk1,k2 ∈ O(3).

Equation (10) is an Orthogonal Procrustes problem and we follow
Gower et al. [18] to solve Ôk1,k2 .
Besides finding consistent surface normal pairs via spatial
overlaps, the existing method [30] shows that pixels with strong
correlation in their intensity profiles (an ordered sequence of scene
irradiance at a pixel across images) have the same surface normals.
This observation is proved to be valid under distant directional
lighting. However, given natural illumination, correlated intensity
profiles do not necessarily lead to consistent surface normals. A
counter-example is a constant environment map, i.e., L(ω) = c,under which all surface normals have correlated intensity profiles.
To approximately extend intensity profile constraint to naturallighting, we propose a consistent orthogonality condition. For
two disconnected scene points q and s as shown in Fig. 5(a),if their surface normals and equivalent directional lighting can betransformed by an orthogonal matrix simultaneously, i.e., Ok2,k3 nk2 (q) Lk2 (q) = nk3 (s) Lk3 (s) ,(11)where Ok2,k3 ∈ O(3) is the orthogonal transformation, thenO = I and both surface normals and equivalent directional lighting for q and s should be consistent. Please refer to Appendix A
for a detailed analysis.
However, what we know from Sec. 3.2 are pseudo equivalent
lighting and surface normals of scene points, with unknown
orthogonal ambiguities to the corresponding ground truth. So we
extend the consistent orthogonality condition to the pseudo normal and lighting case. Assume the pseudo equivalent directional
lighting and surface normals of q and s can be aligned by anorthogonal matrix Ôk2,k3 , i.e.,i hihÔk2,k3 n̂k2 (q) L̂k2 (q) = n̂k3 (s) L̂k3 (s) .
(12)

𝑶𝑘1
෡ 𝑘1 (Θ)
𝑵

𝑶𝑘1,𝑘2𝑵𝑘2 (Θ)෡ 𝑘1,𝑘2
𝑶𝑶𝑘2
෡ 𝑘2 (Θ)
𝑵

𝒏𝑘2 (𝐪), 𝑳𝑘2 (𝐪)𝑶𝑘2
ෝ 𝑘2 (𝐪), 𝑳෠ 𝑘2 (𝐪)𝒏(b)𝑶𝑘2,𝑘3
෡ 𝑘2,𝑘3𝑶

𝒏𝑘3 (𝐬), 𝑳𝑘3 (𝐬)𝑶𝑘3ෝ 𝑘3 (𝐬), 𝑳෠ 𝑘3 (𝐬)𝒏(c)Fig. 5: Illustration of consistent surface normal clustering. (a)Consistent surface normal pairs are extracted from overlapping
patch region Θ and scene point pair (q, s) satisfies consistentorthogonality condition. (b) The relative orthogonal transformation between patch Nk1 and Nk2 is calculated following the relationship between surface normals in the overlapping patch region
Θ (Eqs. (9) and (10)). (c) The relative orthogonal transformation
between patch Nk2 and Nk3 is extracted based on the consistentorthogonality condition between scene points q and s (Eqs. (13)and (14)).

Following the relationship shown in Fig. 5(c), the truth surface
normals and equivalent lighting between q and s can be simultaneously aligned by

Ok2,k3 = Ok3 Ôk2,k3 O>
k2 ∈ O(3),(13)

where Ok2 and Ok3 are the orthogonal ambiguities of the surface
patches Nk2 and Nk3 covering point q and s. Since Eq. (13)makes the consistent orthogonality condition true (Ok2,k3 ∈O(3)), the truth surface normals at scene points q and s are consistent, i.e., Ok2,k3 = I, hnk2 (q), nk3 (s)i = 0. Therefore, Eq. (12)
is an extended consistent orthogonality condition to cluster consistent surface normals from pseudo surface normals and pseudoequivalent directional lighting.
Based on Eq. (3), if the surface normals and equivalent lightingof two scene points fit to the consistent orthogonality condition,
their intensity profiles are correlated. Therefore, to cluster consistent normals on the whole surface, we first filter scene point
pairs with correlated intensity profiles, and then check whether
the pseudo surface normals and equivalent directional lighting ofeach filtered point pair satisfy Eq. (12).
Similar to Ôk1,k2 , Ôk2,k3 = O>
k3 Ok2 also encodes the relationship of orthogonal ambiguities between two surface patches.
To calculate relative orthogonal transformation between patch
Nk2 and Nk3 , we minimize the following energy function as an
Orthogonal Procrustes problem [18].

Ô∗k2,k3 = argmin kÔk2,k3 Dk2 (q) − Dk3 (s)k2F ,Ôk2,k3(14)

s.t. Ôk2,k3 ∈ O(3),where Dk2 (q) = [n̂k2 (q) L̂k2 (q)] and Dk3 (s) follows thesame definition.
To summarize, we collect consistent surface normal pairs fromoverlapping patch regions and scene points satisfying consistentorthogonality conditions. Based on these consistent surface normal

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

7

pairs, we extract the relative orthogonal transformations which
describe the relationship between unknown patch-wise orthogonal
ambiguities. All relative orthogonal transformations form a set
So , which will be applied as edges in the following orthogonal
ambiguity graph building process.

where r is the index of the root node. Following Eq. (16), we
define the smoothness term as
∞di dj 6= di,jE2 (di , dj ) =,(19)
1 − ci,jdi dj = di,j4.2

Constructing Orthogonal Ambiguity Graph

We create an orthogonal ambiguity graph Go = {V, E} (where V
is the set of all nodes and E is the set of all edges connecting
nodes) to build the connections among patches. As shown in
Fig. 2, the nodes of the orthogonal ambiguity graph are filled withthe unknown orthogonal ambiguities Ok from all patches. The
relationship between orthogonal ambiguities can be represented byrelative orthogonal transformations, as shown in Eq. (9). Therefore
we apply all the elements in relative orthogonal transformation set
So to build the edges E of Go . Given surface normal estimation
confidence ci and cj of patch Ni and Nj calculated from Eq. (8),we further define the edge confidence as ci,j = ci cj . Intuitively,
if the normal estimations of two patches are reliable, we tend totrust the relative orthogonal transformation between them.
Based on the orthogonal ambiguity graph, we optimize thepatch-wise orthogonal ambiguities via the following minimization:
XO∗1 , · · · , O∗p = argminµ(O>j Oi − Oi,j ),O1 ,··· ,Op i,j∈E(15)s.t. Oi ∈ O(3).
where µ(·) is a distance measure between two orthogonal matricesin O(3). Directly solving Eq. (15) is non-trivial, so we decompose
the orthogonal ambiguity O into two parts: binary ambiguity
d = |O| ∈ {+1, −1} and rotation ambiguity R ∈ SO(3). Correspondingly, the orthogonal ambiguity graph can also be divided
into binary ambiguity graph Gb and rotation ambiguity graph Gr asshown in Fig. 2. Based on these two graphs, we recover the patchwise orthogonal ambiguities by solving their binary ambiguity partand rotation ambiguity part one after another.
4.3

Optimizing Binary Ambiguity Graph

In binary ambiguity graph Gb , the node value di and the edge
value di,j are calculated from the determinate of the orthogonal
ambiguity Oi and the relative orthogonal transformation Oi,j ,respectively. Following Eq. (15), the binary ambiguities existingin nodes should satisfyX{d∗1 , · · · , d∗p } = argmin(di dj − di,j )2 ,d1 ,··· ,dp i,j∈E(16)s.t. di ∈ {−1, 1}.
Equation (16) can be interpreted as assigning each node of theundirected graph Gb a label defined on {−1, 1}. Therefore we
formulate the problem as maximum a posteriori estimation ofbinary MRF [53], with the energy function defined as
XX
E(d) =
E1 (di ) + η
E2 (di , dj ),
i∈V
i,j∈E
(17)
s.t. di ∈ {−1, 1},
where coefficient η is used to balance the data term E1 and thesmoothness term E2 , i and j represent the node index. We define
the node with maximum degrees in Gb as the root node and set its
binary ambiguity value as 1, then our data term is defined as∞di = −1, i = r
E1 (di ) =,(18)0otherswhere ci,j is the confidence of the edge connecting i-th and j -th
node. Given the definition of the data term and the smoothness
term, we minimize the energy function Eq. (17) with TRW-S
algorithm [29]. Note that, since the binary ambiguity in the root
node could be either −1 or 1, the solved binary ambiguities in all
nodes can only be optimized up to a global binary ambiguity.

4.4

Optimizing Rotation Ambiguity Graph

With binary ambiguity solved, the orthogonal ambiguity in eachnode is reduced to rotation ambiguity. Guided by Eq. (15), we
solve the rotation ambiguity via the following optimization:{R∗1 , · · · , R∗p } = argmin

P

R1 ,··· ,Rp i,j∈E

χ(µ(R>
j Ri , Ri,j )),s.t. Ri ∈ SO(3),(20)

where µ(·) is a distance measure between two rotations in SO(3)and χ(·) is a loss function defined over this distance measure.
This optimization belongs to the rotation averaging problem [21].
Similar to Sec. 4.3, we fix the rotation ambiguity of the root
node as identity, and follow Chatterjee et al. [12] to optimize
the rotation ambiguity in each node. During the rotation averaging
optimization, we apply geodesic distance measurement for µ(·)
and choose Cauchy loss function rather than `2 loss function forχ(·) to improves the robustness when outliers exist in relativerotation transformation Ri,j . Since the true rotation ambiguity of
the root node is unknown, we can only solve per-patch rotation
ambiguities up to a global rotation ambiguity.
After solving rotation ambiguities, we rotate all the patchwise pseudo surface normals and average the normals in theoverlapping regions to get a complete pseudo surface normal map
N̂. Compared to the ground truth, pseudo surface normal map N̂
has two ambiguities left: a global binary ambiguity and a global
rotation ambiguity. We combine the two ambiguities as a global
orthogonal ambiguity Og .

4.5

Resolving Global Ambiguity

So far, estimated N̂ contains only a global ambiguity w.r.t. thetrue surface normal map. This ambiguity can be reduced to aconvex/concave ambiguity by forcing integrability constraint assuggested in [32]. The corresponding proof and the detailed steps
for estimating the global ambiguity can be found in the appendix.
The remained binary convex/concave ambiguity in our surface
normal estimation result could be easily removed manually.

5

E XPERIMENTAL R ESULTS

We first use synthetic data to verify the quantitative accuracyof our method, followed by a comparison between the newly
proposed graph-based patch merging method (GPM) and ourprevious matrix-based patch merging method (MPM) [35]. Finally,
we show the comparison with existing methods on real-world data.

8

B UNNYB EAR

S PHERE

Env. maps

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCEGT Surface normal

Rendered images

Fig. 6: Synthetic dataset. Environment maps (visualized as light probes) from sIBL Archive are shown in the top row. Below we show
ground truth normals for three objects in the first column and examples of rendered images in other columns corresponding to theenvironment maps above.

GT. shading

SV-SH

Ours

Global SH

Env. Map

0.03

0.10

1.89

0.5

0.5

5

0

0

0

Fig. 7: Comparisons between different environment lighting approximation model shown in Table 1. The top row shows theshading maps and the bottom row provides the absolute error maps
and the mean absolute error value of approximated shadings.

5.1

Synthetic Data Setup

We collect 31 real-world environment maps from the sIBL
Archive1 as natural illumination sources, which include diverse
natural illumination from both indoor and outdoor scenarios. We
use Blender [15] as rendering engine and choose three objects –S PHERE, B EAR (from [50]) and B UNNY (with increasing geometric complexity) – to render Lambertian reflectance with white
albedo under natural illumination. The image resolution of thethree objects is fixed to 160 × 160. Ground truth surface normals
and sample images in our synthetic dataset are shown in Fig. 6.
5.2

example environment map and its corresponding shading under
Lambertian reflectance.
Taking the ground-truth shading and surface normal as input,
we extract 3 × 3 patches and calculate our equivalent lighting
direction for each patch. Then we assign it as our approximatedlighting direction at the patch center. To compare with SH-based
lighting models, we render the shading map with our approximatedlighting directions as shown in the third column of Fig. 7. At
the same time, we calculate the second-order global SH lighting
coefficients to approximate the shading given the ground truth
surface normal. We also divide the image into patches and estimate
SV-SH lighting to approximate the patch shading in a similarmanner to our SV-directional model. The absolute error maps w.r.t.
the ground-truth shading are shown in the second row of Fig. 7,
revealing that the SV-SH model and our equivalent lighting model
have close lighting approximation accuracy, and both models are
more accurate than the global SH lighting.

Representation Power of Lighting ModelAs shown in Table 1, the lighting models used in existing methods
include global SH [20] and SV-SH [33]. We compare the representation power of these two models with our SV-directional
equivalent lighting model. Figure 7 shows the comparison on an
1. http://www.hdrlabs.com/sibl/archive.html

5.3

Lighting Model Verification

The local surface normal estimation in our method requires theillumination on a patch to be directional light. Theoretically, if thesurface normals within a patch have the same direction, its natural
illumination is equivalent to a single directional light. However,when surface patches contain diverse normal directions, it isunclear whether a single equivalent lighting direction represents
the patch illumination accurately. In Sec. 3.1, we have defined
the mean angular difference of surface normals (denoted as vn )to evaluate the normal variations in a local patch. Similarly,
we can also define the mean angular difference of equivalent
lighting directions (denoted as vl̄ ) corresponding to the patch
surface normals. This metric can be seen as the error using asingle equivalent lighting direction to approximate the natural
illumination within the patch. As shown in Fig. 8, we provide
a statistic analysis for vl̄ w.r.t. that of surface normals vn on local
patches. The green bar indicates the median value, the top andbottom bounds of the black box indicate the first and third quartile

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

9

vl̄

15

20

Ours (MPM)

10

20

Ours (GPM)#Light

0

        Mean angular difference of surface normals: vn

Normal estimation error




ObjectS PHERE
B EAR



B UNNY


        Mean angular difference of surface normals: vn

Fig. 8: Evaluation of equivalent lighting model. Top row provides
the mean angular difference of equivalent lighting directions vl̄
within a 3 × 3 patch w.r.t. that of surface normals vn in thecorresponding range shown in x-axis. The bottom row provides
the mean angular error of patch normal estimation w.r.t. vn .

values, and the top and bottom ends of the vertical blue line
indicate the minimum and maximum mean angular difference ofequivalent lighting direction. Generally, if a local patch has a larger
normal variation, its illumination is less accurately approximatedby a single equivalent lighting direction.
We also investigate the influence on local surface normal
estimation accuracy when we treat the patch illumination as
an equivalent directional light. Given a surface patch, we first
approximate the patch illumination with the equivalent lighting
direction of the patch center, then estimate the patch surface
normals with this approximated lighting and calculate the mean
angular error w.r.t. the ground truth. The second row of Fig. 8
shows the statistic summary of patch surface normal estimation
error w.r.t. the mean angular difference of patch surface normals.
Although larger surface normal variation will make our equivalent
lighting model approximation more difficult, the patch surface
normal estimation errors from the approximated lighting direction
remain at a low level ((mean angular error < 1.5◦ ).
5.4

Performance under Varying Lighting Conditions

We provide the evaluation of Ours (MPM) [35] and Ours (GPM)on synthetic data under varying numbers of environment lights.
As shown in Fig. 9, we select 10 and 15 subsets of environment
maps out of 20 in our dataset to test how the normal estimation
accuracy varies with lighting conditions. From the angular error
distributions of the S PHERE object, by increasing the image observations under varying natural lights, the surface normal estimation
errors become smaller. The table shown in Fig. 9 further provides
the evaluation of MPM [35] and GPM on all three synthetic
objects. The error values become larger for the B EAR and B UNNY

MethodOurs (MPM)
Ours (GPM)
Ours (MPM)
Ours (GPM)
Ours (MPM)
Ours (GPM)

f = 10
8.87
4.47
17.13
8.93
18.09
10.44

f = 15
6.42
3.29
15.20
8.39
16.78
9.01

f = 20
3.97
2.75
14.78
7.89
15.21
8.69

Fig. 9: Comparison between different patch merging methods (MPM & GPM) under varying numbers of lights (10, 15,
and 20). The top two rows show the angular error distributions
from Ours (MPM) and Ours (GPM) of the S PHERE object. The
table below provides the mean angular errors (in degree) w.r.t.
to the ground truth shown in Fig. 6. Our newly proposed GPM
outperforms our previous MPM [35] on all the three objects.

compared to the smooth shape of S PHERE, which is caused by
the difficulty in approximating equivalent directional lighting on
shape patches with rapid normal variation. The mean angularerrors shown in the table tell that generally a larger number of
input images and more diverse lighting distributions lead to moreaccurate surface normal recoveries. We have also tried further
increasing the number of environment maps up to 31, but the
improvement is rather unobvious, so we fix the number of input
images as 20 for the experiments on synthetic data hereafter.
We also observe that under varying lighting conditions andobject shapes, Ours (GPM) has a smaller mean angular error
compared to Ours (MPM) [35]. It verifies that compared to
the local shortest path searching strategy used in MPM, GPM’s
global optimization on all connections among patches via MRF
optimization and rotation averaging can achieve more accurate
surface normal estimation results.
5.5

Ablation StudyAs shown in Fig. 2, our method mainly contains three stages:
local surface normal estimation, patch merging including MRF
optimization and rotation averaging, and global ambiguity determination. Taking the B EAR as an example, we analyze the errorof estimated surface normal maps from each stage.
In the first stage, the patch-wise surface normals are solved up
to local ambiguities. We first resolve these orthogonal ambiguitiesby aligning the estimated surface normal to the ground truthin each patch and then merge aligned patch normals to build
a complete surface normal map. The angular error map of thissurface normal map shown in Fig. 10(a) verifies that normal
estimation error brought by the first stage is 1.63◦ . We can seethat inaccurate local surface normal estimates mainly occur atregions with large normal variations, such as the neck and leg

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE❶

Patch normal
estimation❷ Patch merging

1.63

6.95

❸

7.32

10

Solving global
ambiguity

OT12 [38]HP19 [19]

2.95

4.79

HY19 [20] Ours (MPM) Ours (GPM)

7.89

S PHERE

10

(a)

(b)
GT binary map(c)Est. binary map

(d)

6.24

4.69

2.75
40

0

6.2%

(f)

(g)

Fig. 10: Ablation study on the B EAR case. The top row showsthe error map of normal estimates at each stage of Ours (GPM).
The bottom row gives the binary ambiguity estimation, where theblue and green pixels correspond to binary ambiguity {1, −1},respectively. The difference between (d) and (e) is shown in (h).
About 6.2% patches have wrong binary ambiguity estimation.

parts of the B EAR, where the equivalent lighting approximations
are inaccurate. These inaccurate local surface normal estimatesfurther influence the rotation ambiguity and binary ambiguity
estimation in the second stage. Figure 10(b) and (c) show theangular error maps of GPM’s surface normal estimation (up toa global ambiguity) with estimated and the ground truth binary
ambiguities shown in Fig. 10(e) and (f). After resolving the patchwise orthogonal ambiguities with MRF optimization and rotation
averaging, surface normal estimation error has increased to 7.32◦ ,in which MRF optimization contributes 0.37◦ and rotation averaging contributes 5.32◦ . As shown in Fig. 10(g), 6.2% of the patches
have wrong binary ambiguity estimation and they are mainly
distributed at regions with inaccurate local normal estimates. Our
complete solution achieves 7.89◦ , where the error brought bysolving the global ambiguity with integrability is 0.57◦ .

5.6

Comparison with Existing Methods

Based on the synthetic dataset shown in Fig. 6, we compare
our method with existing methods [20] (denoted as “HY19”)
and [19] (denoted as “HP19”) on surface normal estimation.
HY19 [20] is the state-of-the-art method for uncalibrated photometric stereo under natural lighting. HP19 [19] also recovers
detailed shape from a rough depth map taking image observations
under unknown natural light as reference. As both HY19 [20]and HP19 [19] require depth initialization, we apply OT12 [38]to generate rough shapes from input images and use them toinitialize HY19 [20] and HP19 [19]. As shown in Fig. 11, thesurface normal estimates from HP19 [19] are influenced by its
initial shape generated by visual hull [38]. HY19 [20] returns
better results from the shape initialization. However, its normal
estimates’ accuracy is still limited by the global SH lighting
approximation. In comparison, Ours (GPM) adopts a more flexible
SV-directional equivalent lighting model and merges local patchesby optimizing all patch connections, therefore achieves the mostaccurate surface normal recoveries among all the methods.

28.68

15.99

11.55

14.78

7.89
40

0

B UNNY

(e)

B EAR

0

28.38

16.69

12.28

15.21

8.69
40

0

Fig. 11: Comparison with existing methods on synthetic data
shown in Fig. 6. Even rows show the estimated surface normal
maps, and odd rows show the corresponding angular error maps
and the mean angular error values in degree.

5.7

Influence of SV-albedos and Shadows

We have tested our method on objects with uniform albedo.
Our method can also be applied to objects with spatially-varying
albedos as long as abrupt albedo changes are not observed within
the patch. In Fig. 12, the image observations of the first twocolumns are rendered with non-uniform albedo maps. Comparedto the normal estimates for uniform albedo shown in Fig. 11,
the piecewise constant albedo distribution only increases theestimation error from 8.69◦ to 9.31◦ . These errors are mainly
brought by the patches across the albedo variation edges. Our
method fails to output accurate surface normal estimates (error
becomes 34.70◦ ) for general spatially-varying albedo distribution
as uniform albedo assumption is not valid for all patches.
We also evaluate the influence of cast shadows. Comparing
the error distribution with and without the cast shadow, we
observe less accurate surface normal estimates (mean angularerror increases to 11.28◦ ) around the B UNNY’s neck and footregions, where cast shadows bring errors to the local surface
normal estimation. On the other hand, when the environment maps
include abrupt changes such as high-frequency light sources, thelighting directions for a surface patch have more variations as thevisibility hemispheres of two surface normals may include/exclude
a high-frequency (the extreme case is a single point light source)
light source. As shown in the last column of Fig. 12, we add 100

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

Cast
shadow

Env. map w/
abrupt changes

GT

Image observations

AlbedoEnv. maps

Piecewise
Generaluniform albedo SV-albedo

11

MPM

GPM

20.83

17.60

Shading

Shading

OT12 [38] HP19 [19] HY19 [20] HJ15 [25]

9.31

34.70

11.28

14.59

26.62

31.00

23.88

27.62

Err. maps

Err. maps

40

40

0
0

5.8

Real-world Experiment

We evaluate and compare our method using real-world datafrom [62] (denoted as “YY13”), [25] (denoted as “HJ15”) andHP19 [19]. The data from YY13 [62] are captured by fixing therelative position between the target objects and the camera whilemoving the whole setup to different places with different natural
illumination, and the data from HJ15 [25] are captured within
one day in an outdoor environment; both datasets and methods
have a mirror sphere to calibrate the environment maps, but suchinformation is not used in our method. The dataset from HP19 [19]includes 8 challenging scenes with complex object shapes andalbedos. Each data includes 20 high resolution (1280 × 720)
images captured under uncalibrated daylight and moving point
light sources.
Surface normal estimation results using OWL (66 images) object from [25] is shown in Fig. 13. We use the ground truth normal
provided by the authors and make a quantitative comparison withexisting methods. Surface normal estimates from HP19 [19] andHY19 [20] have large mean angular errors due to the inaccurate
initial shape from OT12 [38]. The result obtained from MPM [35]generally looks noisier, but it is quantitatively better than thecalibrated result from HJ15 [25], especially in local regions near
the OWL’s eyes where HJ15 [25] shows large errors. Compared to
MPM [35], GPM further improves the surface normal estimation
accuracy at the body and contour region of the OWL, since all the
local patch connections are globally optimized during the MRF
optimization and rotation averaging.
We show the shape estimation results using H ORSEHEAD (7images), C HEF (multi-albedo, 7 images), and M OTHER & BABY
(10 images) objects from YY13 [62] in Fig. 14. Since we do nothave the ground truth for these data, we can only qualitatively
compare our results with them by integrating estimated normal
fields to the depth map with [42]. Referring to the geometry

YY13

Ours (MPM)

Ours (GPM)

HORSEHEAD

Input images

CHEF

small synthetic point light sources to all of the 20 environment
maps. As the illumination for a local patch cannot be treated as
equivalent directional lighting, the normal estimation error has
increased from 8.69◦ to 14.59◦ .

Fig. 13: Quantitative comparison with real data from HJ15 [25].
The numbers on the top of error maps are mean angular errors.

MOTHER&BABY

Fig. 12: The accuracy of Ours (GPM) is influenced by nonuniform albedos, cast shadows and environment maps with abruptchanges.

Fig. 14: Qualitative comparison with real data from YY13 [62].

recovering results, the shape of the H ORSEHEAD and the C HEF
recovered from YY13 [62] are near flat, and the estimated headpart of the M OTHER & BABY by YY13 [62] is distorted. Compared
to YY13 [62], both MPM and GPM produce more visuallyplausible shape estimates without knowing anything about thelighting conditions. Besides, on the neck part of the H ORSEHEAD,
the body part of the C HEF and the arm and baby part of the
M OTHER & BABY, GPM further outperforms MPM [35].
Figure 15 shows the shape recovery results of 8 challenging
scenes from HP19 [19]. PH17 [41], HP19 [19] and HY19 [20]directly produce depth output with given initial object shape, andachieve visually plausible shape recoveries on all the eight scenes.
Compared with these three methods, MPM and GPM are freeof depth initialization. Following Quéau et al. [42], we integrate
estimated surface normal from GPM and MPM to the depth.
The results show that both methods obtain reasonable results on
FACE 1, FACE 2 and BACKPACK. Compared to the state of the art
HY19, our GPM produces comparable shape estimation in thecase of S HIRT and even better result on the VASE by providingIEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

PH17 [41] HP19 [19]HY19 [20]

MPM

GPM

FACE 1BACKPACK

Input

through a graph-based patch merging method (GPM), which
introduces MRF optimization and rotation averaging to solve
the patch-wise ambiguities up to a global orthogonal ambiguity
for the whole surface. Finally, we resolve the global ambiguity
to become a concave/convex ambiguity, which could be easily
removed manually. We believe such a method has great potential
to bring photometric 3D modeling techniques from lab setup with
controlled lighting to wild and large datasets on the Internet.

RUCKSACKTABLETCASEOVENMITT

VASE

SHIRT

FACE 2

6.1

Fig. 15: Qualitative comparison on real data from HP19 [19].
GPM achieves comparable results with PH17 [41], HY19 [19] andHP19 [20] on scenes above the dotted line. The three scenes below
the dotted line are the failure cases for both GPM and MPM. Notethat PH17 [41], HY19 [19] and HP19 [20] require a shape prioras initialization, while GPM and MPM avoid that requirement anddirectly estimate surface normal from images.

richer geometry details. The TABLETCASE and OVENMITT are
two flat objects with SV-albedos. Surface normals of these two
objects have nearly the same directions within a local patch,which is a degenerate case for solving uncalibrated photometricstereo (Sec. 3.2). Also, the complex shape of RUCKSACK brings
large surface normal variations in local patches, which violates theequivalent directional lighting assumption. Therefore, both GPM
and MPM fail on these three scenes.

6

12

C ONCLUSION

We propose a uncalibrated photometric stereo method under
unknown natural illumination. Our method simplifies the natural
illumination using the equivalent directional lighting model thatis valid for local patches. We then solve each patch up to anarbitrary orthogonal ambiguity. The patches are further unified

Limitation and Future Work

Limitation. The limitation of our method mainly lies in the local
surface normal estimation process, in which three assumptions
need to be satisfied: uniform albedo, patch surface normals having
small angular difference, and non-planar surface with 6+ distinct
surface normals. Our method cannot handle complex albedo
maps (e.g., general SV-albedo shown in Fig. 12 and TABLETCASE
in Fig. 15) when the uniform albedo assumption becomes invalid
for most patches. Also, if the surface normals vary significantly
or cast shadow exists in a local patch, the equivalent directional
lighting approximation is less accurate, which leads to wrong
shape estimation results such as the RUCKSACK as shown in
Fig. 15. When given near planar surface such as the case of
TABLETCASE and OVENMITT, our method is also not reliable
due to the degeneration in the local surface normal estimation.
Future work. We hope to develop a more robust local surface
normal estimation method that can handle more complex illumination, texture and cast shadows. To deal with the degenerate casesof near-planar local patches, we can detect them by calculating
the numerical rank as discussed in Appendix C. As this detection
process is not stable due to the shadow and noise in the realcaptured images, we left it as one of our future works. Besides,
to resolve the linear ambiguity in local patches as discussed inSec. 3, an alternative way is applying integrability rather thanuniform albedo constraint to solve the pseudo normals up to a
GBR ambiguity [9], [63]. But how to solve the patch-wise GBR
ambiguities up to a global GBR ambiguity is beyond the scope ofcurrent method. Finally, our patch-wise processing shares similar
spirits with local shading analysis in [61], where they find local
shapes have simple parametric approximation under directional
lighting. In contrast, we explore how local shapes simplify natural
illumination representation. It might be interested to combine local
shape constraints in [61] to further narrow the solution space.

7

ACKNOWLEDGMENT

This work is supported by JSPS KAKENHI Grant Number
JP19H01123, National Natural Science Foundation of China underGrant No. 62136001, 62088102, 61872012. Sai-Kit Yeung ispartially supported by an internal grant from HKUST (R9429).

R EFERENCES
[1]
[2]
[3]
[4]

A. Abrams, C. Hawley, and R. Pless. Heliometric stereo. In Proc. ofEuropean Conference on Computer Vision (ECCV), 2012.
J. Ackermann, F. Langguth, S. Fuhrmann, and M. Goesele. Photometricstereo for outdoor webcams. In Proc. of IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2012.
J. Ackermann, M. Ritz, A. Stork, and M. Goesele. Removing theexample from example-based photometric stereo. In Proc. of European
Conference on Computer Vision (ECCV), 2010.
J. Ackermann, M. Ritz, A. Stork, and M. Goesele. Removing the example
from example-based photometric stereo. In Proc. of ECCV Workshop
RMLE, 2010.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE[5]

N. Alldrin and D. Kriegman. Toward reconstructing surfaces with
arbitrary isotropic reflectance: A stratified photometric stereo approach.
In Proc. of International Conference on Computer Vision (ICCV), 2007.
[6] N. Alldrin, S. Mallick, and D. Kriegman. Resolving the generalized basrelief ambiguity by entropy minimization. In Proc. of IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2007.
[7] R. Basri, D. Jacobs, and I. Kemelmacher. Photometric stereo with general
unknown lighting. International Journal of Computer Vision, 72(3):239–
257, 2007.
[8] R. Basri and D. W. Jacobs. Lambertian reflectance and linear subspaces. IEEE Transactions on Pattern Analysis and Machine Intelligence,25(2):218–233, 2003.
[9] P. Belhumeur, D. J. Kriegman, and A. L. Yuille. The bas-relief ambiguity.
International Journal of Computer Vision, 35(1):33–44, 1999.
[10] M. Brahimi, Y. Quéau, B. Haefner, and D. Cremers. On the wellposedness of uncalibrated photometric stereo under general lighting. In
Advances in Photometric 3D-Reconstruction, pages 147–176. Springer,
2020.
[11] J. Carmigniani, B. Furht, M. Anisetti, P. Ceravolo, E. Damiani, andM. Ivkovic. Augmented reality technologies, systems and applications.
Multimedia tools and applications, 51(1):341–377, 2011.
[12] A. Chatterjee and V. M. Govindu. Robust relative rotation averaging. IEEE Transactions on Pattern Analysis and Machine Intelligence,40(4):958–972, 2017.
[13] L. Chen, Y. Zheng, B. Shi, A. Subpa-Asa, and I. Sato. A microfacet-based
refectance model for photometric stereo with highly specular surfaces. In
Proc. of International Conference on Computer Vision (ICCV), 2017.
[14] D. Cho, Y. Matsushita, Y.-W. Tai, and I. S. Kweon. Photometric stereo
under non-uniform light intensities and exposures. In Proc. of European
Conference on Computer Vision (ECCV), 2016.
[15] B. O. Community. Blender - a 3D modelling and rendering package.
Blender Foundation, Stichting Blender Foundation, Amsterdam, 2018.
[16] P. E. Debevec. Rendering synthetic objects into real scenes: Bridging
traditional and image-based graphics with global illumination and highdynamic range photography. In Proc. of SIGGRAPH, 1998.
[17] A. R. Farooq, M. L. Smith, L. N. Smith, and S. Midha. Dynamic
photometric stereo for on line quality control of ceramic tiles. Computers
in industry, 56(8-9):918–934, 2005.
[18] J. C. Gower. Generalized procrustes analysis. Psychometrika, 40(1):33–
51, 1975.
[19] B. Haefner, S. Peng, A. Verma, Y. Quéau, and D. Cremers. Photometric
depth super-resolution. IEEE Transactions on Pattern Analysis andMachine Intelligence, 2019.
[20] B. Haefner, Z. Ye, M. Gao, T. Wu, Y. Quéau, and D. Cremers. Variational
uncalibrated photometric stereo under general lighting. In Proc. ofInternational Conference on Computer Vision (ICCV), 2019.
[21] R. Hartley, J. Trumpf, Y. Dai, and H. Li. Rotation averaging. International Journal of Computer Vision, 103(3):267–305, 2013.
[22] H. Hayakawa. Photometric stereo under a light source with arbitrary
motion. Journal of the Optical Society of America, 11:3079–3089, 1994.
[23] C. Hernandez, G. Vogiatzis, and R. Cipolla. Multiview photometricstereo. IEEE Transactions on Pattern Analysis and Machine Intelligence,30(3):548–554, 2008.
[24] T. Higo, Y. Matsushita, N. Joshi, and K. Ikeuchi. A hand-held photometric stereo camera for 3-D modeling. In Proc. of International Conference
on Computer Vision (ICCV), 2009.
[25] Y. Hold-Geoffroy, J. Zhang, P. F. U. Gotardo, and J.-F. Lalondey. x-hourphotometric stereo. In Proc. of International Conference on 3D Vision(3DV), 2015.
[26] J. Jung, J.-Y. Lee, and I. S. Kweon. One-day outdoor photometric stereo
via skylight estimation. In Proc. of IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2015.
[27] J. Jung, J.-Y. Lee, and I. S. Kweon. One-day outdoor photometric stereo
using skylight estimation. International Journal of Computer Vision,
127(8):1126–1142, 2019.
[28] I. Kemelmacher-Shlizerman and R. Basri. 3d face reconstruction from asingle image using a single reference face shape. IEEE Transactions on
Pattern Analysis and Machine Intelligence, 33(2):394–405, 2010.
[29] V. Kolmogorov. Convergent tree-reweighted message passing for energy
minimization. IEEE Transactions on Pattern Analysis and MachineIntelligence, 28(10):1568–1583, 2006.
[30] S. J. Koppal and S. G. Narasimhan. Clustering appearance for scene
analysis. In Proc. of IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2006.
[31] S. Z. Li. Markov random field models in computer vision. In Proc. ofEuropean Conference on Computer Vision (ECCV), 1994.
[32] F. Lu, Y. Matsushita, I. Sato, T. Okabe, and Y. Sato. Uncalibrated
photometric stereo for unknown isotropic reflectances. In Proc. of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
[33] R. Maier, K. Kim, D. Cremers, J. Kautz, and M. Nießner. Intrinsic3d: High-quality 3d reconstruction by joint appearance and geometry

13

optimization with spatially-varying lighting. In Proc. of International
Conference on Computer Vision (ICCV), 2017.
[34] D. Miyazaki, K. Hara, and K. Ikeuchi. Median photometric stereo as
applied to the segonko tumulus and museum objects. International
Journal of Computer Vision, 86(2):229–242, 2010.
[35] Z. Mo, B. Shi, F. Lu, S.-K. Yeung, and Y. Matsushita. Uncalibrated photometric stereo under natural illumination. In Proc. of IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2018.
[36] Y. Mukaigawa, Y. Ishii, and T. Shakunaga. Analysis of photometricfactors based on photometric linearization. Journal of the Optical Societyof America, 24(10):3326–3334, 2007.
[37] T. Okabe, I. Sato, and Y. Sato. Attached shadow coding: estimating
surface normals from shadows under unknown reflectance and lighting
conditions. In Proc. of International Conference on Computer Vision(ICCV), 2009.
[38] M. R. Oswald, E. Töppe, and D. Cremers. Fast and globally optimal single view reconstruction of curved objects. In Proc. of IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2012.
[39] T. Papadhimitri and P. Favaro. A new perspective on uncalibrated
photometric stereo. In Proc. of IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2013.
[40] T. Papadhimitri and P. Favaro. A closed-form, consistent and robust
solution to uncalibrated photometric stereo via local diffuse reflectancemaxima. International Journal of Computer Vision, 107(2):139–154,
2014.
[41] S. Peng, B. Haefner, Y. Quéau, and D. Cremers. Depth super-resolution
meets uncalibrated photometric stereo. In Proc. of International Conference on Computer Vision Workshop (ICCVW), 2017.
[42] Y. Quéau, J.-D. Durou, and J.-F. Aujol. Variational methods for normal
integration. Journal of Mathematical Imaging and Vision, 60(4):609–
632, 2018.
[43] Y. Quéau, F. Lauze, and J.-D. Durou. A l1-tv algorithm for robust
perspective photometric stereo with spatially-varying lightings. In
International Conference on Scale Space and Variational Methods in
Computer Vision, pages 498–510, 2015.
[44] H. Santo, M. Waechter, M. Samejima, Y. Sugano, and Y. Matsushita.
Light structure from pin motion: Simple and accurate point light calibration for physics-based modeling. In Proc. of European Conference on
Computer Vision (ECCV), 2018.
[45] F. Shen, K. Sunkavalli, N. Bonneel, S. Rusinkiewicz, H. Pfister, and
X. Tong. Time-lapse photometric stereo and applications. 33(7):359–
367, 2014.
[46] L. Shen and P. Tan. Photometric stereo and weather estimation using
internet images. In Proc. of IEEE Conference on Computer Vision and
Pattern Recognition (CVPR), 2009.
[47] B. Shi, K. Inose, Y. Matsushita, P. Tan, S.-K. Yeung, and K. Ikeuchi.
Photometric stereo using internet images. In International Conference
on 3D Vision (3DV), 2014.
[48] B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan. Self-calibrating
photometric stereo. In Proc. of IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2010.
[49] B. Shi, P. Tan, Y. Matsushita, and K. Ikeuchi. Bi-polynomial modeling of
low-frequency reflectances. IEEE Transactions on Pattern Analysis andMachine Intelligence, 36(6):1078–1091, 2014.
[50] B. Shi, Z. Wu, Z. Mo, D. Duan, S.-K. Yeung, and P. Tan. A benchmark
dataset and evaluation for non-lambertian and uncalibrated photometricstereo. IEEE Transactions on Pattern Analysis and Machine Intelligence,
2018.
[51] W. M. Silver. Determining shape and reflectance using multiple images.
Master’s thesis, MIT, 1980.
[52] K. Sunkavalli, T. Zickler, and H. Pfister. Visibility subspaces: uncalibrated photometric stereo with shadows. Proc. of European Conference
on Computer Vision (ECCV), 2010.
[53] R. Szeliski, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov, A. Agarwala, M. Tappen, and C. Rother. A comparative study of energy
minimization methods for markov random fields with smoothness-basedpriors. IEEE Transactions on Pattern Analysis and Machine Intelligence,30(6):1068–1080, 2008.
[54] P. Tan, S. P. Mallick, L. Quan, D. Kriegman, and T. Zickler. Isotropy,
reciprocity and the generalized bas-relief ambiguity,. In Proc. of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2007.
[55] S. Ullman. The interpretation of structure from motion. Proceedings of the Royal Society of London. Series B. Biological Sciences,203(1153):405–426, 1979.
[56] R. Woodham. Photometric method for determining surface orientationfrom multiple images. Optical Engineering, 19(1):139–144, 1980.
[57] 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.
of Asian Conference on Computer Vision (ACCV), 2010.
[58] T. Wu and C. Tang. Photometric stereo via expectation maximiza-

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

14tion. IEEE Transactions on Pattern Analysis and Machine Intelligence,32(3):546–560, 2010.
[59] T.-P. Wu, K.-L. Tang, C.-K. Tang, and T.-T. Wong. Dense photometricstereo: A markov random field approach. IEEE Transactions on PatternAnalysis and Machine Intelligence, 2006.
[60] Z. Wu and P. Tan. Calibrating photometric stereo by holistic reflectancesymmetry analysis. In Proc. of IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2013.
[61] Y. Xiong, A. Chakrabarti, R. Basri, S. J. Gortler, D. W. Jacobs, andT. Zickler. From shading to local shape. IEEE Transactions on PatternAnalysis and Machine Intelligence, 37(1):67–79, 2015.
[62] L.-F. Yu, S.-K. Yeung, Y.-W. Tai, D. Terzopoulos, and T. F. Chan.
Outdoor photometric stereo. In IEEE International Conference on
Computational Photography (ICCP), 2013.
[63] A. Yuille, D. Snow, R. Epstein, and P. Belhumeur. Determining generative models of objects under varying illumination: Shape and albedo
from multiple images using SVD and integrability. International Journalof Computer Vision, 35(3):203–222, 1999.

Feng Lu received the B.S. and M.S. degrees
in automation from Tsinghua University, in 2007
and 2010, respectively, and the Ph.D. degree
in information science and technology from The
University of Tokyo, in 2013. He is currently aProfessor with the State Key Laboratory of Virtual Reality Technology and Systems, Schoolof Computer Science and Engineering, Beihang
University. His research interests including 3D
shape recovery, reflectance analysis, and human gaze analysis.

Heng Guo received his B.E. and M.S. degrees
in signal and information processing from University of Electronic Science and Technology
of China in 2015 and 2018, respectively. He iscurrently working toward the Ph.D. degree in
Osaka university. His research interests include
physics-based vision and machine learning.

Zhipeng Mo received his B.E. degree from Tianjin University in 2014 and Ph.D. degree fromSingapore University of Technology and Designin 2019, under the supervision of Prof. Sai-Kit
Yeung and Prof. Boxin Shi. He joined Simon
Fraser University as a postdoctoral fellow. His
research interest is photometric methods in computer vision.

Boxin Shi Boxin Shi received the BE degreefrom the Beijing University of Posts and Telecommunications, the ME degree from Peking University, and the PhD degree from the University ofTokyo, in 2007, 2010, and 2013. He is currently
a Boya Young Fellow Assistant Professor andResearch Professor at Peking University, where
he leads the Camera Intelligence Group. Before
joining PKU, he did postdoctoral research withMIT Media Lab, Singapore University of Technology and Design, Nanyang Technological University from 2013 to 2016, and worked as a researcher in the National
Institute of Advanced Industrial Science and Technology from 2016
to 2017. He won the Best Paper Runner-up award at International
Conference on Computational Photography 2015. He has served as an
editorial board member of IJCV and an area chair of CVPR/ICCV. He isa senior member of IEEE.

Sai-Kit Yeung is currently an Associate Professor at the Division of Integrative Systems andDesign (ISD) at the Hong Kong University of Science and Technology (HKUST). Before joining
HKUST, he was an Assistant Professor at theSingapore University of Technology and Design(SUTD) and founded the Vision, Graphics andComputational Design (VGD) Group. During histime at SUTD, he was also a Visiting AssistantProfessor at Stanford University and MIT. Prior
to that, he had been a Postdoctoral Scholar in
the Department of Mathematics, University of California, Los Angeles(UCLA). He was a visiting student at the Image Processing ResearchGroup at UCLA in 2008 and the Image Sciences Institute, University
Medical Center Utrecht, the Netherlands in 2007. He received his
PhD degree in Electronic and Computer Engineering, MPhil degree in
Bioengineering, and BEng degree (First Class Honors) from HKUST,
respectively in 2009, 2005 and 2003.

Ping Tan is an associate professor with the
School of Computing Science at Simon Fraser
University (SFU). Before that, he was an associate professor at National University of Singapore (NUS). He obtained his PhD degree fromthe Hong Kong University of Science and Technology (HKUST) in 2007, and his Master andBachelor degrees from Shanghai Jiao Tong University (SJTU), China, in 2003 and 2000 respectively. His research interests include computer
vision, computer graphics, and robotics. He hasserved as an editorial board member of the IEEE Transactions on
Pattern Analysis and Machine Intelligence (PAMI), International Journalof Computer Vision (IJCV), Computer Graphics Forum (CGF), and theMachine Vision and Applications (MVA), and served as area chairs for
CVPR/ICCV, SIGGRAPH, SIGGRAPH Asia, and IROS. He is a senior
member of IEEE.

Yasuyuki Matsushita received his B.S., M.S.
and Ph.D. degrees in EECS from the Universityof Tokyo in 1998, 2000, and 2003, respectively.
From April 2003 to March 2015, he was with
Visual Computing group at Microsoft Research
Asia. In April 2015, he joined Osaka University
as a professor. His research area includes computer vision, machine learning and optimization.
He is/was an Editor-in-Chief for International
Journal of Computer Vision and on editorial
board of IEEE Transactions on Pattern Analysis
and Machine Intelligence (TPAMI), The Visual Computer journal, IPSJ
Transactions on Computer Vision Applications (CVA), and Encyclopediaof Computer Vision. He served/is serving as a Program Co-Chair of
PSIVT 2010, 3DIMPVT 2011, ACCV 2012, ICCV 2017, and a General
Co-Chair for ACCV 2014 and ICCV 2021. He is a senior member ofIEEE.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

1

Patch-based Uncalibrated Photometric Stereo
under Natural Illumination
Heng Guo† , Zhipeng Mo† , Boxin Shi∗ Senior Member, IEEE, Feng Lu Member, IEEE, Sai-Kit
Yeung Member, IEEE, Ping Tan Senior Member, IEEE, and Yasuyuki Matsushita Senior Member, IEEEFA PPENDIX AD ISCUSSION OF INTENSITY PROFILE BASED SUR FACE NORMAL CLUSTERING

Assume surface normals at two scene points q and s are n(q) andn(s). The equivalent lighting of these two scene points under fenvironment maps areZl̄t (q) =Lt (ω)ωdω,
Ω(q)
Z
(1)
t
l̄ (s) =Lt (ω)ωdω,
Ω(s)
t

3

where L (ω) : R → R represents light intensity of t-th environment map under a unit direction ω ∈ S2 ⊂ R3 , Ω(q) and Ω(s) are
the visible hemispheres corresponding to the two surface normals,Ω(q) = {ω | n> (q)ω ≥ 0},Ω(s) = {ω | n> (s)ω ≥ 0}.

(2)

As shown in Fig. 14(d-e), Since n(q) and n(s) are two unit
direction, there must exists an orthogonal matrix O ∈ O(3) suchthat On(q) = n(s). As a result, the two corresponding visible
hemispheres can also be aligned by O,Ω(s) = {ω | n> (s)ω ≥ 0}
= {ω | n> (q)O> ω ≥ 0}

(3)>

= {Oω | n (q)ω ≥ 0}.
In other word, for any ω ∈ Ω(q), Oω ∈ Ω(s). If the consistentorthogonality condition is satisfied on these two scene points, both
•

•
•
•
•
•
•

H. Guo and Y. Matsushita are with the department of Multimedia Engineering, Graduate School of Information Science and Technology, Osaka University, Japan. E-mail: {heng.guo, yasumat}heng.guo@ist.osakau.ac.jp.
Z. Mo and P. Tan are with the School of Computing Science, Simon Fraser
University, Canada. E-mail: {zhipeng mo, pingtan}@sfu.ca
B. Shi is with National Engineering Laboratory for Video Technology,
Department of Computer Science and Technology, Peking University,
China. E-mail: shiboxin@pku.edu.cn.
F. Lu is with State Key Laboratory of Virtual Reality Technology and Systems, School of Computer Science and Engineering, Beihang University,
China. E-mail: lufeng@buaa.edu.cn.
S.-K. Yeung is with the Division of Integrative Systems and Design and theDepartment of Computer Science and Engineering, Hong Kong Universityof Science and Technology, Hong Kong. E-mail: saikit@ust.hk† H. Guo and Z. Mo contributed equally to this work.
∗ B. Shi is the corresponding author.

surface normal and equivalent lighting can be aligned by the same
orthogonal matrix, i.e.,On(q) = n(s)Ol̄t (q) = l̄t (s) ∀t ∈ (1, f ).
(4)

 >
O O=I
With consistent surface normals and equivalent lighting, the intensity profiles of scene point q and s are correlated. Therefore,
correlated intensity profiles between scene points are the necessary
condition of the consistent orthogonality condition.
Combining Eq. (1) with Eq. (4) we haveOl̄t (q) − l̄t (s)ZZ=O
Lt (ω)ωdω −Lt (ω)ωdω
ZΩ(q)
ZΩ(s)
t
=O
L (ω)ωdω −Lt (Oω)(Oω)dω
Ω(q)
Ω(q)
Z
=O
[Lt (ω) − Lt (Oω)]ωdω

(5)

Ω(q)

= Oδ l̄t = 0.
We denote δ l̄t as equivalent differential lighting, which represents
the spherical integral of the differential environment lighting intensity over the visible hemisphere Ω(p), as shown in Fig. 14. Since
O is an invertible matrix, the consistent orthogonality condition
leads to zero equivalent differential lighting.
Obviously, when surface normals at q and s are consistent,
zero equivalent differential lighting can be satisfied as they havethe same equivalent lighting. Under the case of n(q) 6= n(s), theconsistent orthogonality condition requires all the f environment
maps illuminating these two scene points satisfying the following
condition:Z[Lt (ω) − Lt (Oω)]ωdω = 0, ∀t ∈ (1, f ).
(6)
Ω(q)

It is difficult to analytically prove that unequal surface normal
pairs cannot satisfy the consistent orthogonality condition, sinceEq. (6) is related to the light intensity Lt (ω) from f unknown
environment maps. As shown in Fig. 14, real-world environment
maps are natural HDR images without following any regular
distribution. Also, as we increase the environment lighting amountf , Eq. (6) becomes more difficult to achieve. Therefore, we
provide a statistical analysis on real-world environment maps toverify whether inconsistent surface normal paris may satisfy theconsistent orthogonality condition.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

𝑳2 (𝝎)

𝑳1 (𝝎)

𝑳1 (𝑶𝝎)2

𝑳2 (𝑶𝝎)𝑳𝑓 (𝝎)

𝑳𝑓 (𝑶𝝎)

⋯

⋯

(a)
𝒏(𝒑)

⋯

(b)
Ω𝒑

𝑳𝑡 (𝝎) − 𝑳𝑡 (𝑶𝝎)

𝒏(𝒔)

𝑶

(d)

(c)

(e)

Ω𝒔

𝛿 𝒍1ҧ 𝒏(𝒑) 𝛿 𝒍ҧ𝛿 𝒍ҧ 𝑓
⋯

2

Ω𝒑

(f)

Mean angular error (degree)

Fig. 14: (a) and (b) visualize the f environment maps of two visible hemispheres under distinct surface normals shown in (d) and (e).
(f) shows the equivalent differential lighting δ l̄1 ∼ δ l̄f defined by the spherical integral of environment lighting intensity difference
illustrated in (c).

0.08
0.06
0.04
0.02
0.00

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Number of environment maps

Fig. 15: Surface normal clustering error w.r.t. different environment lighting numbers f , with the consistent orthogonality condition satisfied. The statistics of angular errors are displayed using
the box-and-whisker plot: The green bar indicates the medianvalue, the top and bottom bounds of the black box indicate thefirst and third quartile values, and the top and bottom ends of thevertical blue line indicate the minimum and maximum errors.

To begin with, we first define a “consistent orthogonal error” d(q, s) to evaluate whether the light and surface normal of scene points q and s fit to the consistentorthogonal condition. It can be defined as the mean angular error between vector set [n(s), l̄1 (s), l̄2 (s), · · · , l̄f (s)] and[On(q), Ol̄1 (s), Ol̄2 (s), · · · , Ol̄f (s)], where O is an orthogonal matrix that aligning the equivalent lights and surface normals
between the two scene points. We consider the two scene points
satisfy the consistent orthogonal condition if their consistentorthogonal error is less than a setting threshold.
We collect 31 real-world environment maps from sIBL
Archive, and uniformly sample 151256 distinct surface normal
directions from a sphere and pre-compute the equivalent lightingsfor each normal direction of all the 31 environment maps. We
set the threshold for consistent orthogonal error as 0.01◦ andsummarize the mean angular error of surface normals satisfying
the consistent orthogonal condition w.r.t. different numbers ofenvironment lights in Fig. 15. As an example, in the case of 15
natural lightings, we randomly sample 15 out of 31 environment
maps for 20 times to obtain 20 different environment map groups.
For each group, we first extract scene point pairs whose surface
normal and equivalent lighting directions satisfy the consistent orthogonality condition. Then we record the surface normal angularerror of matched scene point pairs. The mean angular error on
all the 20 groups are only 0.027 degrees, which are quite small
numbers.
From the statistic in Fig. 15, the surface normal clustering
errors with the consistent orthogonality condition on different
number of environment maps are near zero. Therefore, From apractical point of view, it is sufficiently safe to say that underreal-world natural lighting, if surface normals and the equivalentdistant lighting directions of two scene points satisfy the consistentorthogonality condition, the two surface normals and the equivalent lightings have the same directions.

A PPENDIX B
S OLVING G LOBAL A MBIGUITY WITH I NTEGRABILITY
After local surface normal estimation and patch merging process,we can obtain a complete surface normal map up to a global
orthogonal ambiguity. As discussed in [3], [1], [6], this global
orthogonal ambiguity can be reduced to a convex/concave binary
ambiguity by addressing the surface integrability constraint. In thefollowing, we first give the proof and then present the steps tosolve the global orthogonal ambiguity with integrability.

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

B.1

B.2

Proof of Resolving Orthogonal Ambiguity

After merging patch surface normals by optimizing the binary
ambiguity graph and the rotation ambiguity graph, we obtain acomplete surface normal map N̂ up a global orthogonal ambiguity
Og such that

N = Og N̂.

(7)

Following Belhumeur et al. [1], if the surface normal n satisfies
the integrability constraint,then

 
n2
∂ n1
=
,
n3
∂y n3
∂n2
∂n3
∂n1
∂n3
n3
− n2
= n3
− n1
.
∂x
∂x
∂y
∂y

∂
∂x



(8)

>
n̂3 n̂2y − n̂2 n̂3yn̂1 n̂3y − n̂3 n̂1y  n̂2 n̂1y − n̂1 n̂2y 
o1 × o3= 0,n̂2 n̂3x − n̂3 n̂2x 
o2 × o3n̂3 n̂1x − n̂1 n̂3x n̂1 n̂2x − n̂2 n̂1x
(9)

o1 × o3 = kc13
o2 × o3 = kc23 .

(10)

On the other hand, as the rows of orthogonal ambiguity Og , o1 ∼o3 are unit vectors. Therefore, we can solve the scale k up to asign ambiguity. As the determinant of Og can be either 1 or −1,
there are 4 candidates of Og satisfying both orthogonal constraint
and integrability constraint, denoted as Og1 ∼ Og4 below:




c>
−c>
23
23
1
1
 , Og2 = 
,
Og1 =
−c>
c>
13
k −c ×k −c 13
c
×
c
23
13
23
13




>
−c>
c
23
23
1
1
 , Og4 =  −c>  .
Og3 =
c>
13
k c ×
k c ×13cc
23

Ideally, the estimated c13 and c23 from the homogeneous linear
system derived from Eq. (9) should comply with the constraint thatc>13 c23 = 0. Due to the error introduced by the finite difference
and the inaccurate pseudo surface normals included in Eq. (9), theabove constraint cannot be satisfied. Therefore, we formulate an
optimization to address both orthogonal constraint of Og and theintegrability constraint:argmin kAxk22 ,xs.t. x> C>
1 C1 x = 1,

(11)

13

Obviously, we have Og1 = −Og3 , Og2 = −Og4 . As the groundtruth surface normal should have the same direction of camera
view, two of the four orthogonal ambiguity candidates are chosen
to guarantee the recovered shape to be in front of the camera.
Therefore, the global orthogonal ambiguity can be resolved up toa binary choice of the remained two candidates. In the geometryside, this binary ambiguity corresponds to the classical concaveconvex ambiguity that occurs in shape from shading [5]. The
conclusion is also consistent with previous methods [3], [1], [6].

(12)x > C>
1 C2 x = 0,
where C1 = [I, 0], C2 = [0, I], and I ∈ R3 is an identity matrix.
However, the above minimization is non-convex and hard to
optimize. As the O(3) group is compact, we solve the orthogonal
ambiguity in a discrete Hypothesis-and-Test manner.
As the integrability constraint is invariant from the concave/convex ambiguity, we can first solve the rotation ambiguity
from the global orthogonal ambiguity and then choose the correct
one from Eq. (11). Without loss of generality, we decompose therotation ambiguity matrix R into three sub-rotations

R = Rz Ry Rx ,where × denotes the cross product operator, the subscript x, y
represent the partial derivatives in two directions. As the pseudo
surface normal for the complete surface is known, we stack Eq. (9)for all scene points and obtain a homogeneous linear system
Ax = 0. The non-trivial solution of x is then obtained via SVD
on A, result in the cross product estimates c13 and c23 up to ascale k , i.e.

13

How to Solve the Orthogonal Ambiguityx > C>
2 C2 x = 1,

Denoting the three rows of the orthogonal ambiguity Og as o1 , o2 ,and o3 , Substituting Eq. (7) to Eq. (8), we obtain the following
equality:

23

3

(13)

where Rx , Ry , Rz are rotation matrices along x, y, z axes, corresponding to the rotation angle θx , θy , θz , respectively. Suppose theground-truth rotations of Rx , Ry are known, we can get pseudo
surface normal Ñ up to an rotation along the z -axis such thatcos θ − sin θ 0
cos θ 0 Ñ. (14)N = Rz Ry Rx N̂ = Rz Ñ =  sin θ001
Combine Eqs. 9 with Eqs. 14, we have> ñ2 ñ3y − ñ3 ñ2y − ñ3 ñ1x + ñ1 ñ3x
sin θz= 0.
ñ2 ñ3x − ñ3 ñ2x − ñ1 ñ3y + ñ3 ñ1ycos θz(15)

Stacking Eq. (15) for all pixels, we obtain a constrained optimization system about θz as follows

argmin kByk22 ,ys.t. kyk22 = 1,(16)

where y = [sin θz , cos θz ]> . The above optimization can beformulated as a generalized Eigenvalue problem which has aunique solution [2].
Therefore, we sample different pairs of rotations Rx and Ry ,and then solve Rz with the integrability constraint. For each groupof Rx , Ry , Rz , we record integrability cost as kByk22 from theoptimization in Eq. (16). The rotation ambiguity is corresponding
to the group with the smallest integrability cost. As shown in
Algorithm 1, we summarize how to resolve the rotation ambiguity
based on the integrability constraint.
With estimated rotation ambiguity R, we can now build thefour possible orthogonal ambiguity candidates shown in Eq. (11).
Given the prior that the ground-truth surface normals have positivez elements in the viewer-oriented coordinate system, we remove
two of the four candidates, and the remaining orthogonal ambiguity candidates correspond to the convex/concave ambiguity.

4

Est. light

GT. light

Env. maps

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE

0.32

0.29

0.60

0.28

0.36

0.39

0.35

0.25

0.30

10

Err. maps

0.31

0Fig. 16: Accuracy of our equivalent lighting model on 10 environment maps. The lighting directions are visualized in the same way of
surface normal map. The values on the top of error maps are the mean angular errors in degree.

Algorithm 1: Solve rotation ambiguity with integrability
Input : Max rotation angles θxm , θym along x and y axes
Output: Rotation ambiguity R
Initialization: Initial best cost cs
m m
1 for θx ∈ (−θx , θx ) do2for θy ∈ (−θym , θym ) do
3for Ṅ ∈ {N̂, −N̂} do
4Calculate Rx , Ry from θx , θy ;5
Rotate pseudo surface normal toÑ = Ry Rx Ṅ;6
Calculate Rz and record the integrability costc from Eq. (16);7
if c < cs then8cs = c;9
R = Rz Ry Rx ;
10
end
11
end
12
end
13 endA PPENDIX C
N ORMAL E STIMATION FOR P LANAR PATCH
In Sec. 3.2 of the main paper, we have shown that surface normal
can be solved up to an orthogonal ambiguity if there are at least 6
scene points within the patch sharing the same albedo but distinct
surface normals. In this way, there exists a unique solution for the
linear system shown below,>>tri(ñk,1 ñ>tri(Q>k,1 ) · · · tri(ñk,pk ñk,pk )k Qk ) = 1, (17)
{z
}
|
{z
}|
E

y

where ñk is the pseudo surface normal up to a linear ambiguity
Qk in k -th surface patch. When there are no more than 6 diverse
surface normals on the patch (e.g. planar surface), E becomesrank deficient, which reveals the degeneration in our local surface
normal estimation. Such case implies the corresponding surface

GT. Normal

w/o Planar Approx.

w/ Planar Approx.

Surface shape

12.88

4.44

20

20

0

0

Fig. 17: Normal estimation for surface with planar patches asshown in the left column. The middle and right columns show
the surface normal estimates, mean angular errors, and the error
distributions w/o and w/ “planar approximation” strategy.

patch is near flat especially for large patch size (e.g. 5 × 5, 7 × 7).
Therefore, we use a “planar approximation” strategy to force thepseudo surface normals of the patch sharing the same direction
rather than following the estimation in Sec. 3.2 for degenerate
patches. As the ground truth patch surface normal map is alsonear-planar, the pseudo surface normal map as a plane can beapproximately aligned to its ground truth with an orthogonal
matrix. Therefore, we can still apply our GPM to solve the perpatch orthogonal ambiguity including degenerate patches.
As shown in Fig. 17, we shown an object with planar local
patches. For surface normal estimates labeled by “w/ Planar
Approx.”, we first detect degenerate patches by checking whether
E in Eq. (17) is rank-deficient. After that we use “planar approximation” strategy to generate the pseudo surface normals
for degenerate patches and apply Ours (GPM) to obtain thecomplete surface normal estimates. From the error maps shownin the bottom row, the “planar approximation” strategy enables

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCETABLE 2: Comparison on Memory and Time Usage.

Memory Usage(MB/pix)

Time Usage(ms/pix)

MethodS PHEREB UNNY

B EAR

AverageOurs (GPM)
Ours (MPM)

0.029
0.157

0.041
0.189

0.043
0.175

0.038
0.174

Ours (GPM)
Ours (MPM)

11.64
52.35

11.14
27.20

11.16
28.07

11.31
35.87

our method handle the degenerate cases and obtain more accurate
surface normal estimates at the near-planar surface regions.

A PPENDIX D
V ISUALIZATION OF SV-L IGHTINGAs shown in Fig. 16, we visualize the ground-truth per-pixel equivalent lighting directions on a sphere and the approximated lighting
directions from our lighting model. To model the spatially-varying
lighting directions under the natural illumination, our method
assumes the lighting in local patches can be approximated as asingle directional light. As shown in the third row of Fig. 16, eachpixel position encodes the approximated lighting direction for a
local patch centered at that pixel, which is calculated from theshading and surface normals in that patch. From the angular error
maps shown in the last row, the estimated equivalent lighting fromour model is close to the ground-truth lighting directions, whichshows that our method is flexible to model the spatially-varying
environment light.

A PPENDIX E
T IME AND M EMORY C ONSUMING
Table 2 shows the memory and time usage of our method withtwo different patch merging strategies on the three synthetic data
shown in Fig. 6. Given p nodes and q edges of the orthogonal
ambiguity graph, the memory usage in the GPM is 9(p + q)
because the 3D orthogonal matrix has 9 elements. Since thenodes are locally connected in our orthogonal ambiguity graph,it’s reasonable to assume q = kp, k  p, therefore the memory
complexity in GPM is belong to O(9p + 9kp) = O(p). On theother hand, given p pixels, the angular distance propagation matrix
has a dimension of p × p, which leads the memory complexity in
MPM [4] to be O(p2 ). Therefore, our GPM is more memory
efficient compared to MPM [4]. Besides, the angular propagation
matrix in MPM is built from the shortest path searching between
every surface normal pairs, which leads to the time complexity
Cp2 O(p log(p)) = O(p3 log(p)). The time complexity in GPM isrelated to the iterations of rotation averaging and MRF optimization. Therefore it is hard to represent the time cost in a theoreticalway. From the experiments on the three objects shown in Table 2,
GPM runs about 3 times faster than MPM [4] in average.

R EFERENCES[1] P. Belhumeur, D. J. Kriegman, and A. L. Yuille. The bas-relief ambiguity.
International Journal of Computer Vision, 35(1):33–44, 1999.
[2] B. Ghojogh, F. Karray, and M. Crowley. Eigenvalue and generalized
eigenvalue problems: Tutorial. arXiv preprint arXiv:1903.11240, 2019.
[3] F. Lu, Y. Matsushita, I. Sato, T. Okabe, and Y. Sato. Uncalibrated
photometric stereo for unknown isotropic reflectances. In Proc. of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
[4] Z. Mo, B. Shi, F. Lu, S.-K. Yeung, and Y. Matsushita. Uncalibrated photometric stereo under natural illumination. In Proc. of IEEE Conference
on Computer Vision and Pattern Recognition (CVPR), 2018.
[5] J. Oliensis. Uniqueness in shape from shading. International Journal ofComputer Vision, 6(2):75–104, 1991.
[6] T. Papadhimitri and P. Favaro. A closed-form, consistent and robust
solution to uncalibrated photometric stereo via local diffuse reflectancemaxima. International Journal of Computer Vision, 107(2):139–154,2014.

5