Photometric Stereo using Internet Images
Boxin Shi1,2 Kenji Inose3

Yasuyuki Matsushita4 Ping Tan5 Sai-Kit Yeung1 Katsushi Ikeuchi3

1 Singapore University of Technology and Design
3 The University of Tokyo

4 Microsoft Research Asia

Abstract—Photometric stereo using unorganized Internet
images is very challenging, because the input images are captured under unknown general illuminations, with uncontrolledcameras. We propose to solve this difficult problem by a simple
yet effective approach that makes use of a coarse shape prior.
The shape prior is obtained from multi-view stereo and will
be useful in twofold: resolving the shape-light ambiguity inuncalibrated photometric stereo and guiding the estimatednormals to produce the high quality 3D surface. By assuming
the surface albedo is not highly contrasted, we also propose anovel linear approximation of the nonlinear camera responses
with our normal estimation algorithm. We evaluate our method
using synthetic data and demonstrate the surface improvement
on real data over multi-view stereo results.
Keywords-photometric stereo; Internet images; 3D reconstruction; shape prior

I. I NTRODUCTION
Shape recovery is a fundamental problem in computer
vision. Over the past decade, both the capturing devices and
3D reconstruction algorithms have been improved drasticallythat brings surface reconstruction from small scale desktop
objects to large scale outdoor sculptures. Given multiple
images of the same large scale object, for which Internet
is an important image resource, recent progress in structurefrom motion (SfM) and multi-view stereo (MVS) allow
reconstruction even up to city scale. There are existing works
that recover sparse 3D points [1] and depth map [2] for large
scale objects using Internet images. These works focus moreon acquiring the rough depth using geometric constraints
rather than high quality surface.
Photometric stereo, on the other hand, can recover highly
detailed surface geometry at pixel-level accuracy in the formof surface normal map, by using scene radiances observed
under varying lightings [3]. Recently, photometric stereo inan outdoor setting is possible by using a mirror sphere to
calibrate the natural illumination [4]. However, for Internet
images the natural illumination is completely unknown.
When the lighting conditions are unknown, the problem
becomes uncalibrated photometric stereo, whose solution
can only be derived up to some ambiguity, such as thegeneralized-bas-relief (GBR) ambiguity [5] for unknown
directional lightings, or a high-dimensional linear ambiguity [6] for unknown general lightings. Besides the unknown
illumination, uncontrolled sensor is another difficult issue

2 MIT Media Lab
5 Simon Fraser University

since automatic gain control and nonlinear radiometric response deteriorate the resultant shape. Therefore in most of
the previous photometric stereo approaches, sensor gains
and responses are either pre-calibrated or assumed to be
known; however, sensor parameters are usually inaccessible
in Internet images.
In this paper, we focus on Internet images of large outdoor
sculptures which are difficult to capture by Lidar or flying
drone due to their size or security reason. We propose a
unified approach by using the shape prior – coarse shapeinformation obtained by SfM and MVS – to show its important roles in various steps throughout the whole pipeline:from preparing the organized input images for photometricstereo, resolving the ambiguity in uncalibrated photometricstereo to guiding the final surface reconstruction. We also
show that the effect of nonlinear sensor responses can beapproximated by a high-dimensional linear transformation
applied over the illumination component except for highly
contrasted albedos. This can be viewed as pseudo multiplexing of natural lightings which allows highly accurate shapeestimation without the influence of nonlinear responses of
sensors. The key contribution of this paper is to extend
photometric stereo method to work with a wild setting of
unknown general illumination and uncontrolled sensors on
unorganized Internet images.
II. R ELATED W ORK
The Internet images contain comprehensive contents forthe same place of interest from various viewpoints andilluminations. Considering the geometric constraint frommulti-view images, recent progress in SfM and MVS showthat sparse yet reliable 3D points can be recovered fromInternet images [7], [1], [8]. In more recent works, Shanet al. [9] propose a large scale system combining Internet
photos from many resources (ground-level, aerial images,
and street view) for realistic reconstruction, and Zheng et al.
[2] propose a multi-view depth estimation method with consideration of pixel-level view selection. By integrating thephotometric cues, Shen and Tan [10] obtain sparse normals
that are useful for weather estimation. Ackermann et al. [11]apply MVS using Internet images to compute sparse surface
normals and transfer them to images under varying lightingsfor estimating dense normal. By using Internet face photos
as a shape prior and combining shading constraints from

Internet images

Sparse 3D points…

Shape prior
2

1
SfM

Estimated shapeMVS
6

3 3D warping

Depthnormalfusion4 Depth to normal

Registered images

Normal prior…

+

5
Uncalibrated
photometric stereo

Estimated normal

Figure 1. Pipeline of our method, which contains 6 main steps. We take unorganized Internet images as input to first generate a shape (normal) prior,which is then used to produce high quality surface reconstruction with great details.

photometric stereo, the 3D face models can be reconstructedin the wild [12].
The Internet images are highly unorganized and captured in an uncontrolled setup. To apply photometric stereo
on Internet images, the uncalibrated lighting needs to be
estimated. For uncalibrated photometric stereo, it is wellknown that there exists a 3 × 3 linear ambiguity [13] in therecovered surface normals for general surfaces, and a threeparameter GBR ambiguity for integrable surfaces [5]. Recent
works mainly focus on estimating the three unknowns to
obtain final normal estimates (e.g., [14]). Under a general
unknown lighting, there is a 9 × 3 (= 27 unknowns) linear
ambiguity in surface normals under illuminations modeledby second order spherical harmonics. Unfortunately, thishigh-dimensional ambiguity cannot be completely removed
without additional information [6]. We propose to resolve
this ambiguity by using the shape prior.
Another challenging issue for Internet images is thatthe sensors are uncontrollable and their parameters areinaccessible, so methods that require controlled exposure
time (e.g., [15]) are unsuitable to calibrate the radiometricresponse for linearizing Internet images. Self-calibration toradiometric response can be applied for directional lighting [14] or directional plus ambient lighting [16], but for
natural lighting and uncontrolled sensors it is still an open
problem. Instead of explicitly estimating the response function, we disregard it in a self-contained pipeline.

The complete pipeline of our method is shown in Fig. 1.
We collect Internet images of an outdoor sculpture and apply
SfM [7] (Step 1) and PMVS [8] (Step 2) to obtain sparse
point clouds. Then a Poisson reconstruction method [17] isused for creating a water-tight coarse depth prior. Based on
the depth prior we align multiview images to the reference
view via 3D warping (Step 3) to prepare the input forphotometric stereo. A mask is added manually to exclude thesky. Note that in Fig. 1, the Internet images are unorganizedpictures from multiple viewpoints, but the registered images
contain the object from exactly the same viewpoint andvarying natural illumination.
From the shape prior which is a rough depth map, we first
convert it to a rough normal map that gives us the normal
prior (Depth to normal, Step 4). The normal prior and registered images are then combined to solve the uncalibrated
photometric stereo (Step 5) problem. Finally, the estimatednormal and shape prior are integrated to produce the final
3D surface (Depth Normal Fusion, Step 6) which has moredetails due to the accurate normal information.
In the next section, we will explain the details of solving
the uncalibrated photometric stereo problem with the shapeprior, i.e., Step 4-6, especially Step 5; after that, we will
explain how the shape prior helps to relieve the uncontrolledcamera issue in Sec. V.
IV. N ORMAL E STIMATION FROM I NTERNET I MAGES
A. Depth to normal

III. P ROPOSED M ETHOD : OVERVIEW
The main challenges of solving photometric stereo using
Internet images are uncontrolled illuminations and sensors.
The first problem can be formulated as an uncalibrated
photometric stereo with general unknown lightings, and thelatter one is to deal with unknown exposures and radiometricresponses. We tackle these two problems by taking advantages of a coarse shape prior.

Since our method works in the surface normal domain,
we first convert the coarse depth prior into a surface normal
prior for solving our problem (Step 4). A naı̈ve computationof derivatives over a coarse depth map results in a noisy
normal map; therefore, we use a plane principal component
analysis method introduced in [18] for robustly computing
the surface normal prior. Given the depth map and camera
intrinsics of the reference view, the method first projects thedepth map to 3D points in the world coordinate system. For
each 3D point, the method groups a set of points withina short distance d. For the i-th group that contains qi 3Dpoints, by stacking them in a matrix Q ∈ Rqi ×3 , the surface
normal is computed asñ = argmin k(Q − Q̄)nkF ,

(1)

n
qi ×3

where Q̄ ∈ R
is a matrix containing the centroid of
Q in all the rows. A larger d produces a smoother normal
estimate when the input depth contains more noise. The
calculated normals are then projected back to the image
plane of reference view. This step produces the normal priorÑ ∈ Rp×3 where p is the number of foreground pixels in aregistered image.
B. Uncalibrated photometric stereo
1) Image formation model: We begin Step 5 with aLambertian image formation model under natural lightings.
Given a scene point with Lambertian albedo ρ and surface
normal n = [nx , ny , nz ]> , its radiance r can be written as:Zr=
ρL(ω) max((n> ω), 0)dω,
(2)
Ω
3×1

where ω ∈ R
is a unit vector of spherical directions Ω,and L(ω) is the light intensity from the direction ω. This
integration can be approximated using spherical harmonics
asr = s> l,(3)where s = [s1 , s2 , . . . , sk ]> ∈ Rk×1 are harmonics imagesof surface normal n and albedo ρ, and k is the number of elements determined by the order of spherical harmonics. The
vector l ∈ Rk×1 is the k-dimensional lighting coefficients.
Given p pixels observed under q different illuminations,
we store all these p × q radiance values into a radiancematrix R ∈ Rp×q . By a row-wise stacking of p transposed
harmonics images s> in a shape matrix S ∈ Rp×k and acolumn-wise stacking of q lighting coefficients l in a lighting
matrix L ∈ Rk×q , Eq. (3) can be written in a matrix formas:
R = SL.

(4)

We further include the effect of sensor gains and responses
in the image formation model. Under varying illumination,
the exposure time for each image is likely different for an
uncontrolled sensor. Each exposure time corresponds to ascaling of its lighting coefficient l, which is one column
in the lighting matrix L. For simplicity of notations, we
still use L to represent the scaled lighting coefficient matrix.
In addition, a nonlinear response function transforms theradiance R. Let us denote the camera’s radiometric response
as f . For now, we assume the response function f is thesame for all images. The registered images are vectorized

and stacked together in a column-wise manner to form theobservation matrix I ∈ Rp×q . I can be expressed using the
response function f , which is applied in an element-wisemanner using an operator ◦, as:
I = R ◦ f = (SL) ◦ f.

(5)

Our method approximates the nonlinear response function
f using a high-dimensional linear transformation F ∈ Rq×q
as
I = (SL) ◦ f ≈ SLF.

(6)

The transformation F varies with the response function fand radiance R. We will explain and verify the appropriateness of this approximation in Sec. V. Since our goal isto estimate the shape component S, we rewrite Eq. (6) as
I = SLF by LF = LF so that the illumination component
embeds the transformation caused by response functions.
2) Normal estimation algorithm: Similar to previous approaches [13], [6], we perform the singular value decomposition (SVD) on the observation matrix I to estimate the shapematrix S up to a linear ambiguity B ∈ Rk×k . In other words,
the ambiguous S̃ and L̃F are related to their ground truths
S and LF by S̃B = S and B−1 L̃F = LF , respectively.
As discussed in [6], the surface normal is encoded in thesecond to fourth columns of S̃. Therefore, a k × 3 matrix A
is sufficient for computing normal from S̃ as S̃A.
Given the coarse normal prior Ñ from Step 4, we canestimate A to remove the ambiguity:Â = argmin kS̃A − ÑkF ,(7)

A

By applying Â to the original S̃, we obtain disambiguatednormals N̂ by N̂ = O(S̃Â), where O is a row-wisenormalization operator forcing each row of the matrix tobe a unit vector. In practice, we apply Gaussian smoothing
to both Ñ and the ambiguous shape matrix S̃ before solving
for Â.
Solving Eq. (7), however, can only provide a correct
solution if the object’s albedo is uniform. When a scene
contains variant albedos, the norm of rows of S̃Â varies,
while Ñ only contains unit normal vectors. To explicitly
handle the albedo variations, we further optimize A using
the following objective function:A∗ = argmin kO(S̃A) − ÑkF .

(8)A

The above optimization problem is highly nonlinear, but wecan use the linear solution Â as an initial guess to solve
for A. The optimization is solved using a Matlab buildin function “fminsearch”. While the global optimum is notguaranteed, in our experiments this nonlinear refinement
works well because of the good initialization. The final
surface normal is computed by N∗ = O(S̃A∗ ).
The complete normal estimation method (Step 5) is summarized in Algorithm 1.

Algorithm 1 Normal estimation with shape prior
1: Decompose observation matrix I as I = S̃L̃F ;2: Solve the linear equations for Â using Eq. (7);3: Nonlinear refinement to obtain A∗ using Eq. (8);4: Compute normal by N∗ = O(S̃A∗ ).

C. Depth normal fusion
The shape prior is beneficial not only for surface normal
estimation, but also for surface reconstruction (Step 6) byserving as anchor points for the surface recovery from thenormal map [19], [20], [21]. To estimate the optimal depthZ∗ ∈ Rp×1 by combining the estimated surface normal N∗(Step 5) and a vectorized noisy depth map Z ∈ Rp×1 (Step
1 and 2), we can form a linear system of equations as [19]to reconstruct the surface:




λId
λZ
∗
[Z
]
=
,
(9)
∂N∗
∇2
where ∇2 is a Laplacian operator, Id is an identity matrixand λ is a weighting parameter controlling the contribution
∂ ny
∂ nx
of depth constraint. ∂N∗ is the stacks of − ∂xnz − ∂y nz foreach normal n ∈ N∗ . While it forms a large linear systemof equations, because the left matrix is sparse, it can be
efficiently solved using existing sparse linear solvers (e.g.,
QR decomposition based solvers), or multigrid techniques.
V. L INEAR A PPROXIMATION OF S ENSOR R ESPONSES
A useful byproduct of our pipeline is that by using theshape prior, our method naturally ignores the nonlinear
response functions embedded in Internet images. This is
another challenging issue of using Internet images for photometric stereo, due to that for each input image the nonlinear
response function is unknown and arbitrary. Further, theuncontrolled cameras used for recording Internet images
bring difficulty in performing radiometric calibration using
conventional methods [22]. We address this challenge usinga practical approximation by aligning the shape estimateswith the shape prior and encoding the unknown responses
to a linearly multiplexed lighting component.
A. Intuition
The shape estimation method in Sec. IV relies on a highdimensional linear approximation of nonlinear responses(Eq. (6)), which allows us to separate the effects of unknown
sensor gains and responses from the shape estimation as
I = S(LF). The resulting lighting component LF (= LF)
becomes different from the actual L. However, it is a linear
combination of the original lightings, and it can be viewed as
pseudo multiplexing of natural lightings, which allows us toeffectively account for unknown sensor responses. Intuitively
speaking, for each image, the pseudo multiplexing can be
explained as such a process: A uniform surface is illuminatedby a natural illumination l and captured via a nonlinear

response function f that maps r = Sl to i = r ◦ f ; theobserved image is approximately equal to that of the same
surface illuminated by lF (one column of LF ) and captured
with a linear camera, i.e., i = SlF .
Qualitatively, the linear approximation becomes less accurate when a surface contains more diverse albedos. For
example, consider two surface points that have the same
normal n but different albedos ρ1 and ρ2 (ρ1 6= ρ2 ). Witha little bit abuse of notations, we use s = ρn for simplicity
1. The radiance at these two points are r1 = s>
1l =>
ρ1 n> l and r2 = s>
l
=
ρ
n
l,
respectively.
A
camera22
response function maps these radiance values to f (r1 ) andf (r2 ). Since f is a nonlinear function, generally, the ratiof (ρ1 n> l) : f (ρ2 n> l) becomes different from ρ1 : ρ2 .
However, the linear approximation is limited in representing
this nonlinear effect, and this error becomes more obvious
as the difference between ρ1 and ρ2 becomes larger.
For our method, it is not necessary to estimate themultiplexing matrix F; however, the approximation powerof the linear transformation F is of interest because itis related to the shape estimation accuracy. We therefore
assess the appropriateness of the approximation using thedatabase of measured response functions [15]. Fortunately,
as we will see below, the approximation error is consistently
and sufficiently small for the real world response functions
even for variant albedos, due to the high regularity of real
response functions and good approximation capability ofhigh-dimensional linear multiplexing.
B. Verification
We use synthetic images 2 to assess the approximation ability. We simulate the imaging process where f
is applied to R in two different manners: 1) “Nonlinearfixed”: the same response is applied to all images undervarying illuminations. This case corresponds to a scenario
with an uncontrolled camera. And, 2) “Nonlinear-random”:
each image under one lighting condition is distorted by arandomly selected response function in the database. Thiscase more closely mimics Internet images, where each image
is recorded via a distinct unknown and nonlinear response.
We average the results over all 201 response functions in“Nonlinear-fixed” case, and 201 random trials are performed
and averaged for “Nonlinear-random” case. The test sceneconsists of two spheres with different albedos, whose values
are shown at bottom of Fig. 2.
To assess the approximation ability, we evaluate the reconstruction error of Eq. (6). Given R and I, we solve for F bylinear least squares as F̂ = R+ I, where R+ is the pseudoinverse of R. Then, a reconstruction of RF is computed
1 Strictly speaking, s = ρ, s = ρn , s = ρn , s = ρn , s =
x 3
y 4
z 5
1
2
ρ(3n2z − 1), s6 = ρn2x , s7 = ρnx nz , s8 = ρny nz , s9 = ρ(n2x − n2y ) fora second order spherical harmonics representation.
2 We use a 9D spherical harmonics expansion of both normal and lighting
to create these images.

r

rF

i

1

1

0
Err. = 0.31

Err. = 0.004

r

rF

i

1

Brigtness (I)

0.8

i

f

0.6
0.4
0.2
0

0
Err. = 0.49

00.2

Err. = 0.077

0.4
0.6
Radiance (R)

r

0.8

1

Figure 3. Input radiance r is transformed by a nonlinear response f (shown on rightmost) to i as i = r ◦ f . Top row: two spheres with the same albedo(α = 1); bottom row: two spheres with greatly different albedo (α = 5). Our linearly approximated rF shows very close appearance to i. The errors
below show mean of relative differences of r and rF from i, respectively. Color encoded images are used for better visualization.

trast. The high correlation of different radiometric response
functions makes our method works for both “Nonlinearfixed” and “Nonlinear-random” cases (their reconstruction
errors are always similar), i.e., our approximation is valid
for Internet images.
As an intuitive example, we show linearly approximatedimages of uniform abledo (α = 1, top row) and strong
contrast case (α = 5, bottom row) in Fig. 3. Note that iand rF have very small difference visually, especially for
uniform albedo case, which shows the validity of our linear
approximation.

0.3

 = 1/Nonlinear-fixed = 2/Nonlinear-random = 2/Nonlinear-fixed = 2/Nonlinear-random = 3/Nonlinear-fixed = 3/Nonlinear-random = 4/Nonlinear-fixed = 4/Nonlinear-random = 5/Nonlinear-fixed = 5/Nonlinear-random

Reconstruction error

0.25

0.2

0.15

0.1

0.05

0

VI. E XPERIMENTS
A. Quantitative evaluation
5

10

15

20

25

30

35

40

q

α=1

α=2

α=3

α=4

α=5

Figure 2.
Reconstruction errors of Eq. (6) w.r.t. varying numbers ofimages (q) for scenes containing two spheres with different albedos.
α = {1, 2, 3, 4, 5} indicate that left/right spheres have albedo values of{0.5/0.5, 0.4/0.6, 0.3/0.7, 0.2/0.8, 0.1/0.9}.

as RF = RF̂. The reconstruction error is evaluated as themean of |i − rF |/i, where i is a pixel observation of I andrF is the corresponding element in RF . This is a relative
error (percentage) defined for each observation. We show thereconstruction errors with respect to the varying numbers of
input images q and albedo contrast in Fig. 2. The errors are
pretty low (about 1%) when the number of input images q
becomes q ≥ 9 for the case of uniform albedo (α = 1).
On the other hand, as the albedo contrast becomes greater,
the errors increase accordingly. Except for the extreme case(α = 5) 3 , the reconstruction errors are consistently low(below 5%). Therefore, the approximation generally works
well, except for scenes that exhibit significantly high con3 α = 5 mimics a scene of almost black and white spheres.

We use a synthetic scene, C AESAR, to quantitatively
evaluate our method. The data is synthesized in the same
way as the simulation test did in Sec. V-B. We fix the numberof distinct lightings q = 40 for this test. We evaluate how
the varying albedo contrast and coarseness of shape priorinfluence the normal estimates given nonlinear images, byapplying various real-world response functions [15] to input
images. We found that our method works well for mostalbedo variations except for the highly contrasted one, andeven a severely contaminated shape prior is quite usefulin estimating accurate normal. We provide the complete
analysis in the supplementary material.
6Here, we show
a typical example in Fig. 4. In thisexample, a severely contaminated depth map, quantized to
4 bits with zero-mean Gaussian of standard deviation 0.04
being added, is used as shape prior. The normal estimation
accuracy is evaluated using angular difference (degree) w.r.t.
the true normal. With a noisy normal prior of angularerror 21.3◦ , we obtain a normal map with much smaller
errors, which is about 14◦ smaller than the prior. The
surface reconstruction error is defined as the mean valueof |z0 − z ∗ |/z0 across pixels, where z0 is the true depth
and z ∗ is the depth estimate. The original rough depth(Shape prior) is noisy, but it still provides useful positional

Input images

True surface

True normal

Normal prior Estimated normal

Shape prior

Err. = 21.3o
Normal integ.

Err. = 7.4o
FusionErr. = 0.10

Err. = 0.17

Err. = 0.06

Shape prior

Our result

Shape prior

Our result

Figure 4. Surface normal estimation (top row) and surface reconstruction
result (depth and surface in two bottom rows) using synthetic data. The
numbers on normal maps show the angular errors and the numbers on
surfaces show the reconstruction errors.

information. On the other hand, direct integrating a surface
from the normals results in a distorted reconstruction with alarger bias (Normal integ.). By fusing the normal and depth
information, a more accurate surface can be reconstructed(Fusion), as pointed out by previous work [19].
B. Result using Internet images
In addition to the K AMAKURA B UDDHA data shownin Fig. 1, we show results of three more scenes named
M OTHERLAND C ALLS, S TATUE OF L IBERTY, and M OUNTRUSHMORE in Fig. 5. These four datasets contain 200, 109,
320, and 128 downloaded images respectively, which areused for SfM and MVS. We use 163, 109, 63, and 42 images
in each dataset which roughly have the same viewpoints for
3D warping and normal estimation. We model the natural
illumination using the third order spherical harmonics (k =16) for all the experiments. We compare our method to [11]by using the same input for reference. Our normal estimatesshow more meaningful shape information than the resultsfrom [11], because of the capability of handling naturallightings and variations of camera responses. For example,
in the pedestal of S TATUE OF L IBERTY, our result shows
consistent normals for plane structures and clearer detailsof the bricks. We also show the reconstructed surfaces byfusing our estimated normal and the shape prior from MVS

Figure 6. 3D reconstruction results using Internet images. Close-up views
are indicated by red rectangles.

(as baseline for surface reconstruction comparison) in Fig. 6,
where more details can be observed thanks to the refinedsurface normals by photometric stereo.
When the Internet images have almost the same viewpoints, SfM and MVS will produce degenerated results. But
if a scene contains (partly) regular shapes, we can directly
use this knowledge as the shape prior. We show such an
example of TAJ M AHAL where the shape of the dome has acomprehensive structure. We manually assign a hemispheresurface normal map to the dome part as the shape prior. The
Internet images of TAJ M AHAL are registered to a reference
view via homography using SIFT [23] features in this case,as the 3D shape information is unavailable. The result using
66 images is shown in Fig. 7. The rendered Lambertianshading using the estimated normal under a distant lighting[ √13 , √13 , √13 ]> is shown on the rightmost for verification.
VII. C ONCLUSION
We present a photometric stereo method that works withunorganized Internet images, captured under general unknown illumination, with uncontrolled sensors. We suggest
using shape priors from SfM and MVS to fully remove

Input images

Normal prior

Our results

Results using [11]Figure 5. Surface normal estimation results using Internet images. Four representative images from the input dataset are shown in the leftmost column,
with the image on top left is the reference image to which other images are registered. Other columns show normal prior, normal estimated by our methodand [11], respectively. Close-up views of estimated normal maps are indicated by red rectangles.

the ambiguity in uncalibrated illumination setting, to guidethe normal to surface integration, and to avoid the effectof uncontrolled sensor. The proposed method shows highquality 3D modeling over existing MVS method.
Limitations: In our current solution, cast shadows are
not handled. Due to the shape-light ambiguity, it is difficult
to explicitly calculate the visibility map like [24]. We
have investigated some robust algorithm [25] to handle castshadows as outliers by forcing the input matrix to be rank-k.
However, the result showed almost no improvement in ourcontext, because the ideal rank-k matrix is seldom observedfor Internet images. Properly modeling cast shadows in our
pipeline is left as our future work. With recent progress on
SfM and MVS techniques, it is possible to obtain finer shapepriors. Combining our method with better shape priors using
Internet images, which consider large scale data [9] and view
selection [2], is also an interesting direction.
ACKNOWLEDGEMENTPart of this work was done while Boxin Shi was anintern at Microsoft Research Asia and a Ph.D. candidate
at the University of Tokyo. Boxin Shi is supported by

Registered images

Template normal

Estimated normal

Relighting result

Registered images

Normal prior

Estimated normal

Relighting result

Figure 7.

Surface normal estimation result using Internet images and the known shape as a prior.

SUTD-MIT joint postdoctoral fellowship. Sai-Kit Yeung is
supported by SUTD StartUp Grant ISTD 2011 016, SUTDMIT International Design Center Grant IDG31300106, andSingapore MOE Academic Research Fund MOE2013-T2-1159.
R EFERENCES[1] M. Goesele, N. Snavely, B. Curless, H. Hoppe, and S. Seitz,“Multi-view stereo for community photo collections,” in Proc.
ICCV, 2007.

[12] I. Kemelmacher-Shlizerman and S. Seitz, “Face reconstruction in the wild,” in Proc. ICCV, 2011.
[13] H. Hayakawa, “Photometric stereo under a light source with
arbitrary motion,” Journal of the Optical Society of America,
vol. 11, no. 11, pp. 3079–3089, 1994.
[14] B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan, “Selfcalibrating photometric stereo,” in Proc. CVPR, 2010.
[15] M. Grossberg and S. Nayar, “Modeling the space of camera
response functions,” IEEE Trans. on Pattern Analysis andMachine Intelligence, vol. 26, no. 10, pp. 1272–1282, 2004.

[2] E. Zheng, V. Jojic, E. Dunn, and J.-M. Frahm, “Patchmatch
based joint view selection and depthmap estimation,” in Proc.
CVPR, 2014.

[16] M. Diaz and P. Sturm, “Radiometric calibration using photocollections,” in Proc. ICCP, 2011.

[3] R. Woodham, “Photometric method for determining surface orientation from multiple images,” Optical Engineering,
vol. 19, no. 1, pp. 139–144, 1980.

[17] M. Kazhdan, M. Bolitho, and H. Hoppe, “Poisson surface reconstruction,” in Proc. Eurographics Symposium on Geometry
processing (SGP), 2006, pp. 61–70.

[4] L.-F. Yu, S.-K. Yeung, Y.-W. Tai, and S. Lin, “Shading-based
shape refinement of RGB-D images,” in Proc. CVPR, 2013.

[18] K. Klasing, D. Althoff, D.Wollherr, and M. Buss, “Comparison of surface normal estimation methods for range sensing
applications,” in Proc. ICRA, 2011.

[5] P. Belhumeur, D. Kriegman, and A. Yuille, “The bas-relief
ambiguity,” International Journal of Computer Vision, vol. 35,
no. 1, pp. 33–44, 1999.
[6] R. Basri, D. Jacobs, and I. Kemelmacher, “Photometric stereo
with general unknown lighting,” International Journal ofComputer Vision, vol. 72, no. 3, pp. 239–257, 2007.

[19] D. Nehab, S. Rusinkiewicz, J. Davis, and R. Ramamoorthi,“Efficiently combining positions and normals for precise 3D
geometry,” in Proc. SIGGRAPH (ACM Trans. on Graphics),
2005, pp. 536–543.
[20] N. Joshi and D. Kriegman, “Shape from varying illuminationand viewpoint,” in Proc. ICCV, 2007.

[7] N. Snavely, S. Seitz, and R. Szeliski, “Photo tourism: Exploring photo collections in 3D,” in Proc. SIGGRAPH (ACM
Trans. on Graphics), 2006, pp. 835–846.

[21] T. Higo, Y. Matsushita, N. Joshi, and K. Ikeuchi, “A handheld photometric stereo camera for 3-D modeling,” in Proc.
ICCV, 2009.

[8] Y. Furukawa and J. Ponce, “Accurate, dense, and robust
multiview stereopsis,” IEEE Trans. on Pattern Analysis andMachine Intelligence, vol. 32, no. 8, pp. 1362–1376, 2010.

[22] T. Mitsunaga and S. Nayar, “Radiometric self calibration,” inProc. CVPR, 1999.

[9] Q. Shan, R. Adams, B. Curless, Y. Furukawa, and S. Seitz,“The visual turing test for scene reconstruction,” in Proc.
3DV, 2013.

[23] D. Lowe, “Distinctive image features from scale-invariant
keypoints,” International Journal of Computer Vision, vol. 60,
no. 2, pp. 91–110, 2004.

[10] L. Shen and P. Tan, “Photometric stereo and weather estimation using Internet images,” in Proc. CVPR, 2009.

[24] C. Wu, B. Wilburn, Y. Matsushita, and C. Theobalt, “Highquality shape from multi-view stereo and shading undergeneral illumination,” in Proc. CVPR, 2011.

[11] J. Ackermann, M. Ritz, A. Stork, and M. Goesele, “Removing
the example from example-based photometric stereo,” inProc. ECCV Workshop RMLE, 2010.

[25] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma,“Robust photometric stereo via low-rank matrix completion
and recovery,” in Proc. ACCV, 2010.

Photometric Stereo using Internet Images – Supplementary Material:
Complete Quantitative Evaluation in Sec. 6.1
Boxin Shi1,2 Kenji Inose3

Yasuyuki Matsushita4 Ping Tan5 Sai-Kit Yeung1 Katsushi Ikeuchi3

1 Singapore University of Technology and Design
3 The University of Tokyo

4 Microsoft Research Asia

145
k=9/Nonlinear-fixed
k=9/Nonlinear-random
k=16/Nonlinear-fixed
k=16/Nonlinear-random

4

Angular error (degree)

Angular error (degree)

2 MIT Media Lab
5 Simon Fraser University

32
1

k=9/Nonlinear-fixed
k=9/Nonlinear-random
k=16/Nonlinear-fixed
k=16/Nonlinear-random

12

10

8

6

0

1

α=1

2

α=2

3
Albedo ratio index ( )

α=3

4

α=4

5

α=5

Figure 1.
Normal estimation accuracy (angular error in degrees) w.r.t. different albedo contrasts. α = {1, 2, 3, 4, 5} indicate that left/right half of the object have albedo values of{0.5/0.5, 0.4/0.6, 0.3/0.7, 0.2/0.8, 0.1/0.9}. Two different dimensionsof lighting coefficients k are evaluated.

I. E FFECT OF ALBEDO CONTRAST
We evaluate the effect of albedo contrast to normal
estimation accuracy. To exclude other factors except for
nonlinear sensor responses, in this we use the ground truth
normal as Ñ to remove the ambiguity. Figure 1 shows thenormal estimation accuracy with respect to varying albedo
contrast α for different dimensions of lighting coefficients k.
As we have observed in Sec. 5, the accuracy is affected by
the greater albedo contrast in general, and the errors become
smaller with a larger k. This indicates that the higher-order
lightings make the pseudo multiplexing more effective. The“Nonlinear-fixed” cases show larger errors than “Nonlinearrandom” cases due to large errors caused by some response
functions in unusual shapes that are difficult to approximate.
II. E FFECT OF NOISE IN SHAPE PRIORS
Figure 2 shows the variation of normal estimation errors
with different noise levels in shape priors. The input depth
values are quantized to 3 bits in the worst case, and Gaussian

1

β=1

2

β=2

3
Noise level index ()

β=3

β=4

4

5

β=5

Figure 2. Normal estimation accuracy (angular error in degrees) w.r.t.
varying noise levels in shape priors. β = {1, 2, 3, 4, 5} are labels
to represent corruptions where the clean depth maps are quantized to{8, 6, 5, 4, 3} bits, with zero-mean Gaussian noise of standard deviations{0.02, 0.03, 0.04, 0.04, 0.05} added. Two different dimensions of lighting
coefficients k are evaluated.

7
noise with standard 7deviations
up to 0.05 is further added inorder to simulate the real-world shape priors. The computed
surface normal priors have errors from about 13◦ to 25◦ .
In this test, the albedo is set as uniform to remove theeffect from albedo variations. The errors increase with theroughness of the shape priors. Except for the extreme case(β = 5), the normal estimation accuracy is consistently higheven with nonlinear responses. Under severe noise, a large
k lowers the normal estimation accuracy, because it allows
too much freedom in the ambiguity matrix, which makes thesolution sensitive to noise. In practice, k should be adjusted
according to the tradeoff between nonlinearity of responses(prefers a larger k) and the coarseness of shape prior (prefersa smaller k).