Radiometric Calibration for Internet Photo Collections
Zhipeng Mo1 Boxin Shi2∗ Sai-Kit Yeung1 Yasuyuki Matsushita2,3
1
Singapore University of Technology and Design2
Artiﬁcial Intelligence Research Center, National Institute of AIST
3
Osaka University
Abstract

computer vision problems, when they are to be applied toInternet photos, e.g., intrinsic image decomposition [18]and photometric stereo [29], require input images to be radiometrically linearized. However, a commercial camera
usually maps the scene radiance to its pixel values in anonlinear manner for compressing the dynamic range andaesthetic purpose. Such a mapping is unknown in mostcases and treated as business secrets by camera manufactures. The goal of radiometric calibration is to estimate thecamera’s (inverse) radiometric response function so that theobserved pixel values linearly relate the scene radiance.

Radiometrically calibrating the images from Internet
photo collections brings photometric analysis from lab datato big image data in the wild, but conventional calibration methods cannot be directly applied to such image data.
This paper presents a method to jointly perform radiometriccalibration for a set of images in an Internet photo collection. By incorporating the consistency of scene reﬂectancefor corresponding pixels in multiple images, the proposed
method estimates radiometric response functions of all theimages using a rank minimization framework. Our calibration aligns all response functions in an image set up to thesame exponential ambiguity in a robust manner. Quantitative results using both synthetic and real data show theeffectiveness of the proposed method.

Classic radiometric calibration approaches capture astatic scene under various exposure times [25]. For outdoor scene with illumination changes, given an image sequence from a ﬁxed viewpoint, the problem could also be
solved by physically modeling the image formation and using the consistent scene albedo across images as the key
constraint [15]. Radiometrically calibrating Internet photos
shares a similar spirit to [15], but it is much more challenging due to that 1) the scenes are captured from multiple diverse viewpoints, 2) each image is captured by an unknown
camera with unknown settings (e.g., white balance and exposure time), and 3) a batch of different response functions
need to be estimated simultaneously.

1. Introduction
For a popular landmark, millions of pictures are captured
and shared through the Internet. Such Internet images andcommunity photo collections provide comprehensive image
resources to computer vision research, because they contain images captured from different viewpoints, at different
time, under different illumination, and using different typesof cameras and settings. The cameras capturing these images naturally organize a “network”. By exploring such a
network, geometric analyses such as geometric camera calibration and 3D reconstruction, which are generally infeasible using a single image, become tractable by establishing correspondences for recovering 3D points using multiview observations. Recent progress on structure from motion (SfM) [32] and multi-view stereo (MVS) [9, 10] show
successful applications using Internet photos.
Photometric analysis is an important complement problem for analyzing images from Internet photo collectionsorganized by such a camera network, and radiometric calibration is a key prerequisite for photometric analysis. ManyIn this paper, we propose a new method to jointly perform radiometric calibration to all cameras using a collection of Internet photos for the same scene. We utilize thegeometric information calculated from SfM and MVS, andassume the scene reﬂectance (albedo) is the same for corresponding pixels in all images. Our key observation is thatthe ratio of albedo values for pixels in the same image withthe same surface normal should be equal across the image
set, only when each image is applied with a correctly calibrated inverse radiometric response function. Such a ratio operation cancels the inﬂuence of different white balance settings, exposure times, and environment illumination
conditions which are usually arbitrary for each image (andcamera). The albedo ratios are stacked as a matrix whose
rank should be 1. We minimize the rank of such a matrix
to achieve robust estimation that aligns all response func-∗ Corresponding author: boxin.shi@aist.go.jp.

1tions up to a global exponential ambiguity. Figure 1 illustrates key ideas and the pipeline of our method. The performance of our approach is evaluated using various synthetic
and real-world datasets in a quantitative manner.

2. Related Work
Our work is related to conventional radiometric calibration methods and computer vision applications to Internet
photo collections.
Radiometric calibration. A Gretag Macbeth twenty-four
patch color checker, whose reﬂectance value for each patch
is known, is a commonly used tool for radiometric calibration [2]. Single image calibration without calibration
chart is possible by assuming linear radiance distribution
at edges [21, 22]. Using multiple images with different exposure times is a popular and practical approach for radiometric calibration. Classic methods include Debevec andMalik’s approach ﬁtting a non-parametric, smooth response
function [6], and Mitsunaga and Nayar’s method adopting
a polynomial model [25]. The problem in [25] can be formulated as a rank minimization one to achieve superior robustness [19]. The response functions can also be represented using a more realistic model by using the databaseof measured response functions (DoRF) [11]. Such a representation can also be used in log-space to deal with moving cameras [16], varying illumination conditions [15], anddynamic scenes in the wild [1]. In addition to multiple
exposure constraint, other constraints such as the statistical model of CCD imaging [33], geometry-invariant property [26], symmetric distribution of noise [24], temporal
mixture of motion blur [34], and multiple directional lighting [30], are proposed for different imaging setups and applications. Note in all these approaches only one camera
is calibrated using one image or multiple images, and thecamera is controlled to adjust its settings (e.g., manual modewith ﬁxed white balance, ISO, but varying exposure times)
for calibration purpose.
There are existing works that perform radiometric calibration for Internet photos. Kuthirummal et al. [17] explore
priors on large image collections. By assuming the same
camera model has a consistent response function and someradiometrically calibrated images of that camera model areavailable, a camera-speciﬁc response function could be estimated for the image collection according to the deviationfrom statistical priors. Due to the improvement of 3D reconstruction techniques on Internet-scale image sets, radiometric calibration becomes feasible by using the scene geometry estimated from SfM [32] and MVS [9, 10]. Diazand Sturm [7, 8] jointly solve for the albedos, response
functions, and illuminations by using nonlinear optimization and priors from DoRF. A more recent work by Li andPeers [20] assume a local smoothness of image patch appearances and a 1D linear relationship over corresponding
image patches, so that the radiometric calibration for multiview images under varying illumination could be recast as
the classic multiple exposure one [16], which could potentially be applied to Internet photos.
Computer vision meets Internet photos. Various computer vision problems could be extended to deal with
Internet photos. Successful applications include scene
completion [12], virtual tourism [31], weather estimation [28], composing pictures from sketches [3], image
restoration [5], generating face animations [14], image colorization [4], intrinsic image decomposition [18], color consistency [27], photometric stereo [29], synthesizing timelapse video [23], and so on and so forth. In many of theseapplications, the arbitrary nonlinear response functions for
all cameras are simply approximated as a global Gamma
correction, or even completely ignored as a linear one. So
we believe a radiometric calibration solution for Internet
photos is a very important technique that could potentially
beneﬁt miscellaneous application scenarios relying on photometric analysis.

3. Image Formation Model
We assume the scene reﬂectance follows the Lambertianmodel, and we know the scene geometry (surface normal)
from SfM and MVS. The correspondence between 3D scene
points and 2D pixels in all images are also obtained from 3Dreconstruction. We take the j-th image in the image collection as an example, in which the scene is illuminated under
the j-th natural lighting Lj (ω). Then the scene radiance of
the i-th 3D point is determined by the interaction of lightingwith its surface normal ni ∈ R3×1 scaled by Lambertianalbedo ρi asRij =υij (ω)ρi Lj (ω) max((n(1)i ω), 0)dω,Ωwhere ω ∈ R3×1 is a unit vector of spherical directions Ω,and Lj (ω) is the environment map for the j-th image which
encodes the light intensity from the direction ω. υij (ω) isthe visibility function which is set to 0 if the illuminationfrom direction ω is not visible for the i-th 3D point projected to the j-th image or 1 otherwise. For any ni with visibility function being equal to 1, it receives the light from its
visible hemisphere Ωi , and the integration over the visible
hemisphere is simpliﬁed as
Rij = ρi (n
i l̄j ),

(2)


where l̄j = Ωi Lj (ω)dω.
When a scene is captured by the j-th camera, the image
irradiance for the i-th pixel in the k-th color channel (e.g.,Pixel ratios

…

…

g2

…

gQ

P
212

. .. Q…

…

…
P

Pixel pairs with the same normal but
different albedos

2 12. .. Q

g2

Observation

Rank-1 matrix

Image Q

Image 2

Irradiance

g1

…

…
…
Image 1

g1
Observation…
Image 2

Irradiance

Pixel ratios and rank minimization

Irradiance

3D reconstruction from a photo collection

gQ
Observation

Inverse radiometricresponse functions

Figure 1. Pipeline of our method. We estimate the (inverse) radiometric response functions for each image {g1 , g2 , · · · , gQ } by rank
minimization over the stacks of pixel pairs. The “ ◦ ” operator applies each inverse response function g to both the numerator anddenominator of ratio terms in the same row. The correct g transforms each row of the matrix to the same vector (up to a scale) to make thematrix rank-1.

RGB) can be represented ask= ckj tj ρki (nIij
i l̄j ),

(3)

where ckj is a white balance scale and tj is an exposure time
for the j-th camera. Due to a nonlinear mapping of the camera radiometric response function f (·), the observations are
distorted by f askk
= fjk (Iij
) = fjk (ckj tj ρki (nBij
i l̄j )).

(4)

The response function is a monotonic function, so there exists a unique inverse function g = f −1 to map the observation values to image irradiance values. By applying theinverse response function g to both sides of Eq. (4), we obtaink) = ckj tj ρki (ngjk (Biji l̄j ).

(5)

Radiometric calibration could be performed for three different color channels independently, so we drop the k-related
terms thereafter. Denote qj = cj tj as an image-dependentscaling factor, then Eq. (5) is simpliﬁed asgj (Bij ) = qj ρi (n
i l̄j ).

(6)

In the context of radiometric calibration for Internet photos,
each image in the photo collection has its own g. Our goal
is to simultaneously estimate g for all the images in a photocollection.

4. Radiometric Calibration Method
4.1. Formulation
With the scene 3D information available, it is possible
to ﬁnd points with the same surface normal, receiving thesame amount of light, but with different albedos in each image. We assume such pixels are identiﬁed for now, and the
pixel selection method will be introduced in Sec. 4.4. Leta pair of such 3D points have normal n and lighting l̄, andtheir albedo values be ρm and ρn (ρm = ρn ) respectively.
Substituting these two points into Eq. (6) and taking the ratio between them, we obtaingj (Bmj )qj ρm (n l̄)ρm=.
=gj (Bnj )ρn
qj ρn (n l̄)(7)

The albedo ratio consistency above is the key constraint
we employ for radiometric calibration. Given a sufﬁcientnumber of observation values B that cover a broad intensity
range, we can build a system of equations solved by nonlinear optimization [25]. Recent progress in radiometric calibration shows that the rank minimization could solve sucha problem in a more robust manner and effectively avoid
overﬁtting [19]. Therefore, we formulate our problem in amatrix form, whose minimum rank corresponds to the correct estimates of inverse response functions.
Denote P as the total number of pixel pairs with thesame normal (and lighting) but different albedos, and Q as
the total number of images, or equivalently pairs of lighting
conditions and cameras. We could arrange these pixels according to the ratio format of Eq. (7) and stack them as thefollowing matrix:⎞⎛ g (B )g1 (B21 )1
11
1 (BP 1 )· · · g1g(B
g1 (B01 )g1 (B11 )(P −1)1 )⎟⎜ g2 (B12 )g2 (B22 )
2 (BP 2 )⎟⎜· · · g2g(B
g2 (B12 )⎜ g2 (B02 )(P −1)2 ) ⎟AQ×P = ⎜⎟ . (8)
..
..
..
..
⎟
⎜
.
.
.
.
⎠
⎝
gQ (B1Q )gQ (BP Q )gQ (B2Q )· · · gQ (B(P −1)Q )gQ (B0Q )gQ (B1Q )The optimal inverse response function gj for each row
transforms each pixel ratio in the matrix to its corresponding albedo ratio, so that each row of A becomesρ1 ρ2
ρP
1
ρ0 , ρ1 , · · · , ρ(P −1) , which obviously makes A a rank-1

matrix. Thus, our radiometric calibration problem becomesthe following rank minimization one:∗{g1∗ , g2∗ , · · · , gQ}=

argminrank(A).

(9){g1 ,g2 ,··· ,gQ }The above operation and process is illustrated in Fig. 1.

4.2. Optimization
We solve the above rank minimization using a similar
approach as in [19], which is represented by the condition
number as∗}={g1∗ , g2∗ , · · · , gQ

σ2 (A),{g1 ,g2 ,··· ,gQ } σ1 (A)argmin(10)

where σi (A) is the i-th singular value of A.
We choose to use the polynomial representation for g assuggested by [19]. The main consideration is that polynomial representation is more appropriate for gradient-based
convex optimization because of its smoothness. Both irradiance and observation values are normalized in the rangeof 0 to 1. Then the polynomial representation of g becomes
S−1

pi B S−i−1 ,g(B) = B + B(B − 1)

(11)

i=1

where {p1 , p2 , · · · , pS−1 } are the polynomial coefﬁcientsto be estimated. Such an expression uses only S − 1 unknowns to represent an S-order polynomial. The end point
constraints for inverse response functions are explicitly enforced, since Eq. (11) satisﬁes g(0) = 0 and g(1) = 1.
Note that we only borrow the optimization strategy
from [19] to solve Eq. (10). In fact, our problem is much
more challenging than [19] due to the joint estimation ofmany different response functions, and the structure of thematrix whose rank needs to be minimized is completely different due to pixel ratios. We ﬁnd such a problem cannot
be directly solved like [19] as directly solving for all g simultaneously is quite unstable, because each gj transforms
one row of A independently and this signiﬁcantly increases
the search space for minimum rank. So we solve this issue
by using a pairwise optimization followed by a global
ref inement.
The pairwise optimization means we select two rows as“base” image pair and align all the other rows to the basein an incremental manner. The base image pair is selected
as the two rows of A with the minimum difference after
applying the estimated inverse response functions, through1 For easy representation, we only show one option for arranging the
pixel pairs (there could be ratio terms like ρρ2 , ρρP , etc.). Given p + 1
0
12
different ρ values, there are CP+1 possible combinations of taking theratio.

Algorithm 1 Radiometric calibration algorithm
INPUT: Input images, with pixels selected and stacked
as the matrix of Eq. (8).
// Pairwise optimization:
for all pairwise combinations using two rows of A do
Solve two g using Eq. (10);Apply these two g to the corresponding rows in A;
end for0
Select gm
and gn0 that make corresponding rows of A
have the minimum difference;for k = {1, 2, · · · , Q} ∧ k = {m, n} do
Build a matrix with the {m, n, k}-th rows of A and0
and gn0 ﬁxed;
solve gk0 using Eq. (10), with gm
end for// Global reﬁnement:
Solve Eq. (10) for all g simultaneously using0} as initial values;{g10 , g20 , · · · , gQ∗}.
OUTPUT: Inverse response functions {g1∗ , g2∗ , · · · , gQ

2solving Eq. (10) for all CQ
submatrices composed by tworows of A. Then we add one row at a time to solve forthe remaining Q − 2 rows for submatrices with three rows.
The estimated inverse response functions here are denoted
as g 0 . The global reﬁnement takes g 0 as initial values tosolve for all g simultaneously using Eq. (10). Our complete
algorithm is summarized as Algorithm 1.

4.3. Exponential ambiguity
Similar to existing radiometric calibration methods [1,11, 15, 16, 19, 25], our method also suffers from the exponential ambiguity. This is because for a set of optimized solution g, g γ for any unknown γ also keeps the ratio consistent in Eq. (7) and makes A rank-1. Note there only exists
one γ for all images. The uncertainty of γ brings challenges
to joint estimation of all g and this is part of the reason that
we have to start from a pairwise optimization approach in
Algorithm 1. Although inﬁnite many solutions exist for one
set of images without additional constraint, our method naturally uniﬁes all response functions to the same ambiguity
and such results could be useful for tasks such as high dynamic range imaging and panorama image stitching.
We could remove such a one-parameter ambiguity by inserting only one radiometrically calibrated image into thephoto collection. Without losing generality, we put pixel
pairs from this image in the ﬁrst row and ﬁx g1∗ as a linear
function during the optimization. Then the optimized re∗} automatically align all the othersults for {g2∗ , g3∗ , · · · , gQ
rows of A with the linearized values to form a set of solution free of exponential ambiguity.

4.4. Implementation details
3D reconstruction. To build the matrix in Eq. (8), we
need to extract corresponding pixels in all images that havethe same surface normal, under the same lighting condition,
but with different albedos. We ﬁrst perform 3D reconstruction (SfM [32] and MVS [10]) using the input photo collection. The 3D points with the same surface normal are
selected and projected onto 2D images. We then calculate
pairwise ratio for these pixels as initial guess of albedo ratios. The selected pixels in each pair should receive thesame amount of environment illumination, if their visibility
function υij deﬁned in Eq. (1) were the same. However,the visibility information cannot be accurately estimatedthrough sparse 3D reconstruction and unknown environment illumination, which results in noise in real data suchas cast shadow and local illumination variations. Therefore,
we propose a simple outlier rejection approach to deal withthis issue. We ﬁnd that the majority of such initially selected pixel pairs show similar ratio values, and any noisy
pixels appearing in either numerator or denominator cause
the ratio signiﬁcantly different from others. Such outliers
could be easily identiﬁed and discarded by a line ﬁtting using RANSAC. Finally, remaining pixel pairs observed in all
images are stacked as the matrix in Eq. (8) for optimization.
Details of optimization. In practice, we only require
dozens of images as input for the 3D reconstruction by SfM
and MVS. According to Algorithm 1, it is not necessary to
optimize Q response functions simultaneously, since we use
an incremental approach to estimate all response functions
except for the two that are selected as base. Given a complete set of images, we divide it into several subgroups (e.g.,
10 images in each subgroup), and solve for each subgroupusing Algorithm 1. We empirically ﬁnd such divide-andconquer strategy gives a more stable solution, and this property allows the parallel processing of large amount data.
We use the Matlab build-in function “lsqnonlin” tosolve our nonlinear optimization. The initial guess for inverse response functions are chosen as a linear function for
all images in the pairwise optimization step. A monotonicity constraint is added to penalize non-monotonic estimates
similarly as adopted in [19]. We further add a second-orderderivative constraint by assuming most response functions
have either concave or convex shapes. There are response
functions with more irregular shapes according to [11], butthey are rarely observed in common digital cameras.
Degenerate case. One obvious degenerate case for our
problem is a scene with uniform albedo. Because the uniform albedo causes ρm = ρn in Eq. (7), and every element
in A becomes one and thus A becomes rank-1. Therefore,
we need at least two different albedo values in the scene.

However, if pixels with two different albedos are all on thesame plane (assuming there is no shadow in the scene, sothat the same surface normal receives the same amount ofenvironment lighting), it falls into a similar degenerate case.
Therefore, the minimum requirement to avoid the degenerate case is a pair of different albedos on two different planes.
This is because even if the albedo ratio is the same after optimization, if two pixel pairs are scaled by different shading terms n l and then nonlinearly mapped by the same
response function, their ratios become different before optimization. Fortunately, a wild scene may contain much morevariations in either albedo or normal than our minimum requirement. Since the problem formulation is highly nonlinear, it is non-trivial to provide analytical proof for thenumber of different albedo or surface normal required, but
we will experimentally analyze such an issue in the next
section.

5. Experimental Results
We perform quantitative evaluation of our method using both synthetic and real data. The error metric used forevaluation is the rooted mean square errors (RMSE) of theestimated inverse response function w.r.t. the ground truth
and the disparity, i.e., the maximum absolute difference between the estimated and the ground truth curves.

5.1. Experiment with synthetic data
Number of pixel pairs vs. order of polynomial. Our
method is expected to be more stable and accurate given
more diverse values of pixel pairs (albedo and normal variations) and ﬁtted with higher order polynomials. We use
synthetic data to verify the accuracy under different number of albedo values and polynomial orders by testing 6
types of pixel pairs and 6 different polynomial orders. The 6
groups of input data are generated by changing the numberof different albedo values multiplied by different normals as{1×2, 2×3, 3×4, 4×6, 5×10, 6×15}. Here, 1×2 means
one pair of different albedo values on two different planes.
We then apply 10 different lighting conditions and 10 different response functions from the DoRF database [11] togenerate our observation values. The inverse response functions are estimated using our methods, and the RMSE anddisparity are summarized in Fig. 2.
As expected, the average errors show a row-wise decreasing tendency due to more diverse input data variations.
It is interesting to note that only 24 pairs of points from each
image (four pairs of albedo values on six different planes)
produce reasonably small error (RMSE around 0.01) for a
joint estimation of 10 different response functions. From
Fig. 2, we can also see our method is not sensitive to the
choice of polynomial order and in general a polynomial order larger than 5 works well. We ﬁx the polynomial orderto 7 in all experiments as a tradeoff between accuracy and30.0761
0.06620.14580.07420.06280.0665

40.04330.02240.01670.01440.01590.0194# pixel pairs

60.05660.01710.01940.01300.00660.007970.05540.01490.01700.00940.00950.006780.05220.01680.01490.00750.00750.0073Polynomial order

Disparity
1ൈ2
2ൈ3
3ൈ4
4ൈ6
5ൈ10
6ൈ15

5
0.05390.01850.01620.01440.00760.00943
0.13340.11620.26810.13340.11170.117440.07130.04120.03200.02880.03250.039850.09400.03510.03120.02880.01750.0206

6
0.10140.03290.03770.02600.01520.0175

7
0.09900.02870.03390.01950.02060.0166

8
0.09250.03190.02880.01690.01690.0180RMSE/Disparity values# pixel pairs

1ൈ2
2ൈ3
3ൈ4
4ൈ6
5ൈ10
6ൈ15

Noise level (Cg):

Polynomial order

RMSE

Figure 3. Evaluation under various noise levels and comparison
between our method and Diaz13 [8]. The box-and-whisker plotshows the mean (indicated as “×”), median, the ﬁrst and third
quartile, and the minimum and maximum values for RMSE anddisparity for 200 (20 × 10) estimated response functions.
5HIOHFWDQFH

Figure 2. The average RMSE/disparity w.r.t. number of pixel pairs(row-wise) and order of polynomials (column-wise). “Red” means
larger and “blue” means smaller errors.

6KDGLQJcomplexity.

Results on single-view images. Single-view synthetic
test is performed to provide an intuitive example with quan-

1RUPDOL]HGLUUDGLDQFH












 



(VWLPDWHG 506('LVSDULW\
*URXQG7UXWK




 




:LWK5)

:LWK,5)

2ULJLQDO















          
1RUPDOL]HGREVHUYDWLRQ

1RUPDOL]HGLUUDGLDQFH

Performance with various noise. We ﬁrst add quantization to mimic the 8-bit image formation, and then we add
Poisson noise as suggested in [13, 19], which describes theimaging noise in a more realistic manner by considering
signal-dependent shot and dark current noise. The noiselevel is controlled by the camera gain parameter Cg , andlarger Cg means more severe noise. Please refer to Equations (9)-(11) in [19] for the noise model representation. We
perform 20 trials and each test contains 10 randomly selected response functions from DoRF with 4 × 6 pixel pairs.
We evaluate another photometric image formation based
method [8] (denoted as “Diaz13”) implemented by ourselves using the same data. We ﬁnd a joint estimation toall variables (albedos, lighting and response function coefﬁcients) produces unreliable results, due to the nonlinear optimization over too many variables. Therefore, we
provide the ground truth lighting coefﬁcients in our implementation of Diaz13 and use this as the stable performance of Diaz13 for comparison. The results under various Cg = {0, 0.1, 0.2, 0.5} (where Cg = 0 means only
quantization noise) for both methods are plotted in Fig. 3.
Our method outperforms Diaz13 for most noise levels, but
Diaz13 shows more stable but less accurate performance
under these noise levels. When the noise is large, ourmethod shows degraded performance partially due to thatratio operation magniﬁes the noise. Note that in real case,Diaz13 requires the dense reconstruction for lighting estimation, while we can only work on a few selected pixel
pairs.








(VWLPDWHG 506('LVSDULW\
*URXQG7UXWK








          
1RUPDOL]HGREVHUYDWLRQ:LWK5)

:LWK,5)

2ULJLQDO









Figure 4. Radiometric calibration results using a synthetic dataset.
The upper row shows the ground truth reﬂectance, shading images,
and the ten selected pixel pairs (yellow numbers) with the same
normal but different albedo values. Two example results of the estimated inverse response functions and the ground truth curves areplotted, with the RMSE and disparity values shown in the legend.
The nonlinear observations (“With RF”), linearized images (“WithIRF”) and their absolute difference maps w.r.t. the “Original” images are shown next to the inverse response curve plots.

titative analysis, which is free of errors from 3D reconstruction. We use the data from [18], which is generated using a

0.5
0.4
0.3
0.2
0.1

2

3

Normalized irradiance

1

0.7
0.6
0.5
0.4

0.5
0.4
0.3
0.2

0.3
0.2

0.8
0.7
0.6
0.5
0.4
0.3
0.2

0.8
0.7
0.6
0.5
0.4

Ours (RMSE/Disparity)
0.0558/0.1034
Diaz13 (RMSE/Disparity)
0.0942/0.1659
Lin04 (RMSE/Disparity)
0.1049/0.1412
Calibrated
3

0.3
0.2
0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

1
0.9

Ours (RMSE/Disparity)
0.0179/0.0315
Diaz13 (RMSE/Disparity)
0.0855/0.1645
Lin04 (RMSE/Disparity)
0.1210/0.2218
Calibrated
2

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation
1
Ours (RMSE/Disparity)
0.9 0.0374/0.0668
0.8
0.7
0.6
0.5
0.4

Diaz13 (RMSE/Disparity)
0.0746/0.1210
Lin04 (RMSE/Disparity)
0.0887/0.1577
Calibrated
3

0.3
0.2
0.1

0.1

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation
1
Ours (RMSE/Disparity)
0.9 0.0155/0.0402

00
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation
Normalized observation
1
1
Ours (RMSE/Disparity)
Ours (RMSE/Disparity)
0.9 0.0137/0.0320
0.9 0.0168/0.0303
0.8 Diaz13 (RMSE/Disparity)
0.8 Diaz13 (RMSE/Disparity)
0.7 0.0458/0.0771
0.7 0.0714/0.1348
Lin04 (RMSE/Disparity)
Lin04 (RMSE/Disparity)
0.6
0.6
0.0991/0.1483
0.1086/0.2143
0.5
0.5 Calibrated
3
Calibrated
2
0.4
0.4

0.8
0.7
0.6
0.5
0.4

Diaz13 (RMSE/Disparity)
0.0727/0.1540
Lin04 (RMSE/Disparity)
0.1185/0.2186
Calibrated
1

0.3
0.2

Normalized irradiance3

Normalized irradiance2

0.8

Ours (RMSE/Disparity)
0.0240/0.0473
Diaz13 (RMSE/Disparity)
0.0738/0.1585
Lin04 (RMSE/Disparity)
0.1133/0.2057
Calibrated
1

Normalized irradiance

1
0.9

0.6

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation1

0.7

Normalized irradiance

0.6

0.8

1
0.9

Normalized irradiance

0.7

Ours (RMSE/Disparity)
0.0075/0.0158
Diaz13 (RMSE/Disparity)
0.0835/0.1366
Lin04 (RMSE/Disparity)
0.0956/0.1685
2
Calibrated

0.3
0.2

Normalized irradiance3

Normalized irradiance

2

0.8

1
0.9

Ours (RMSE/Disparity)
0.0117/0.0361
Diaz13 (RMSE/Disparity)
0.1173/0.2338
Lin04 (RMSE/Disparity)
0.1396/0.2354
1
Calibrated

Normalized irradiance

1
0.9

1

0.1

0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

00 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

Figure 5. Estimated inverse radiometric response functions using an image collection mixed with Internet photos and captured images (inred box). The results compared with “Diaz13” [8] and “Lin04” [21] on RGB channels of three images and calibrated ground truth are
shown. The RMSE and Disparity are in the legend of each plot.

physics-based renderer. Given the ground truth reﬂectanceand shading images, as show in the upper part of Fig. 4, we
manually select 10 pairs of pixels with the same surface normal but different albedo values (labeled using yellow numbers on the reﬂectance and shading images). We randomly
apply 10 different response functions to produce 10 observation images, and then use the selected pixels as input toperform radiometric calibration. Two typical results of theestimated inverse radiometric response functions w.r.t. theground truth curves are shown in Fig. 4. We further apply
the inverse response functions to the observation images,
and the close appearances between the linearized images
and the original images show the correctness of our radiometric calibration method.

5.2. Real data experiment
To perform quantitative evaluation using real data, wecreate datasets containing mixture of Internet photos andimages captured using controlled cameras for three different
scenes. The total numbers of images are 55 for the dataset
used in Fig. 5, and 44, 31 respectively for the other two

datasets used in Fig. 6. We use three controlled cameras(x Sony Alpha7, y Nikon D800, and z Canon EOS M2)
for all datasets, and calibrate their response functions using multiple exposure approach [19]. We perform 3D reconstruction and radiometric calibration for all images ineach dataset. Two (out of three) captured images are usedas the base image pair for pairwise optimization, and thecalibrated response functions are used to remove the exponential ambiguity for all estimates from our method. We
also run our own implementation of “Diaz13” [8]2 and [21](denoted as “Lin04”)3 on the same dataset for comparison.
The quantitative results (estimated inverse radiometricresponse functions and their RMSE/Disparity w.r.t. theground truth) on RGB channels of three images (capturedby three different camera models) from the ﬁrst scene areplotted in Fig. 5, and we show similar plots only using the
R channel of the other two scenes in Fig. 6. Except for the
2 This time the lighting coefﬁcients are also part of the optimization
which are initialized randomly.
3 Implemented by Jean-Franois Lalonde and available from: https://github.com/jflalonde/radiometricCalibration

0.6
0.5
0.4
0.3
0.2
0.1

3

0.7
0.6
0.5
0.4
0.3
0.2

1
0.9

0.1

0.8
0.7
0.6
0.5
0.4

Ours (RMSE/Disparity)
0.0196/0.0697
Diaz13 (RMSE/Disparity)
0.1442/0.2634
Lin04 (RMSE/Disparity)
0.1617/0.2854
Calibrated
3

0.3
0.2
0.1

0
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation
Normalized observation
1
1
Ours (RMSE/Disparity)
Ours (RMSE/Disparity)
0.9 0.0214/0.0366
0.9 0.0451/0.0951
0.8 Diaz13 (RMSE/Disparity)
0.8 Diaz13 (RMSE/Disparity)
0.7 0.1199/0.2444
0.7 0.0609/0.0996
Lin04 (RMSE/Disparity)
Lin04 (RMSE/Disparity)
0.6
0.6
0.1279/0.1993
0.0291/0.0918
0.5 Calibrated
0.5 Calibrated
1
2
0.4
0.4
0.3
0.3

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation
1
Ours (RMSE/Disparity)
0.9 0.0132/0.0387

0.2

0.2

Normalized irradiance2

Normalized irradiance1

0.8

Ours (RMSE/Disparity)
0.0205/0.0656
Diaz13 (RMSE/Disparity)
0.1617/0.3227
Lin04 (RMSE/Disparity)
0.1831/0.3435
Calibrated
2

Normalized irradiance

0.7

1
0.9

0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

0.2

Normalized irradiance3

Normalized irradiance2

0.8

Ours (RMSE/Disparity)
0.0239/0.0504
Diaz13 (RMSE/Disparity)
0.0664/0.1293
Lin04 (RMSE/Disparity)
0.0779/0.1315
Calibrated
1

Normalized irradiance

1
0.9

1

0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

0.8
0.7
0.6
0.5
0.4

Diaz13 (RMSE/Disparity)
0.1541/0.3026
Lin04 (RMSE/Disparity)
0.1605/0.2885
Calibrated
3

0.3
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Normalized observation

Figure 6. Estimated inverse radiometric response functions using another two image collections mixed with Internet photos and captured
images (in red box). The results compared with “Diaz13” [8] and “Lin04” [21] on R channel of three images and calibrated ground truth
are shown. The RMSE and Disparity are in the legend of each plot.

second example in the bottom row of Fig. 6 where “Lin04”[21] shows better result than ours, all our estimates show
the closest shape (with minimum RMSE/Disparity) to thecalibrated ground truth. The high accuracy of our method ispartly due to that we use calibrated ground truth to resolve
the exponential ambiguity, but the experiment here provesthat our method is able to jointly estimate the inverse response functions for a set of images up to the same exponential ambiguity reliably.

time settings of individual images, we can only recover thescene radiance up to a linearly scaled value (for each color
channel) through the radiometric calibration. Explicitly inferring the white balance and exposure-time settings of Internet photos is an interesting future topic. Current method
cannot handle images with severe noise, and an even morerobust solution that compensates the noise magniﬁcation issue caused by ratio operation is desired.

6. Discussions

Acknowledgement

We present a method to perform radiometric calibration
for images from Internet photo collections. Compared tothe conventional radiometric calibration problem, the Internet photos have a wide range of unknown camera settings
and response functions, which are neither accessible nor adjustable. We solve this challenging problem by using thescene albedo ratio, which is assumed to be consistent in all
images. We develop an optimization method based on rank
minimization for jointly estimating multiple response functions. The effectiveness of the proposed method is demonstrated using both synthetic and real-world data.
Currently, we need to assume the 3D reconstruction
is sufﬁciently reliable for extracting points that share thesame surface normal. With radiometric calibration problem
solved, it will be interesting to combine with photometric
3D reconstruction using Internet photos [29] to further improve the quality of 3D reconstruction. Although our solution is naturally invariant to white balance and exposure-We thank Joon-Young Lee for providing the code of [19].
Boxin Shi and Yasuyuki Matsushita are supported by aproject commissioned by the New Energy and Industrial
Technology Development Organization (NEDO). Boxin
Shi is also supported by JSPS KAKENHI Grant Number JP17K12722. Sai-Kit Yeung is supported by Singapore MOE Academic Research Fund MOE2016-T2-2154; Heritage Research Grant of the National Heritage
Board, Singapore; SUTD Digital Manufacturing and Design (DManD) Centre which is supported by the National
Research Foundation (NRF) of Singapore; the NRF, Prime
Ministers Ofﬁce, Singapore under its IDM Futures Funding Initiative and its Environmental & Water Technologies Strategic Research Programme and administered by the
PUB, Singapores National Water Agency; this material isbased on research/work supported in part by the National
Research Foundation under Virtual Singapore Award No.
NRF2015VSG-AA3DCM001-014.

References[1] A. Badki, N. K. Kalantari, and P. Sen. Robust radiometriccalibration for dynamic scenes in the wild. In Proc. of Internatoinal Conference on Computational Photography, 2015.
2, 4[2] Y.-C. Chang and J. F. Reid. RGB calibration for color image analysis in machine vision. IEEE Transactions on Image
Processing, 5(10):1414–1422, 1996. 2[3] T. Chen, M.-M. Cheng, P. Tan, A. Shamir, and S.-M. Hu.
Sketch2photo: Internet image montage. ACM Transactions
on Graphics (Proc. of ACM SIGGRAPH Asia), 28(5):124:1–
124:10, 2009. 2[4] A. Y.-S. Chia, S. Zhuo, R. K. Gupta, Y.-W. Tai, S.-Y. Cho,
P. Tan, and S. Lin. Semantic colorization with internet images. ACM Transactions on Graphics (Proc. of ACM SIGGRAPH Asia), 30(6):156:1–156:8, 2011. 2[5] K. Dale, M. K. Johnson, K. Sunkavalli, W. Matusik, andH. Pﬁster. Image restoration using online photo collections.
In Proc. of Internatoinal Conference on Computer Vision,
2009. 2
[6] P. E. Debevec and J. Malik. Recovering high dynamic range
radiance maps from photographs. In Proc. of ACM SIGGRAPH, 1997. 2[7] M. Diaz and P. Sturm. Radiometric calibration using photocollections. In Proc. of Internatoinal Conference on Computational Photography, 2011. 2[8] M. Diaz and P. Sturm. Estimating photometric properties
from image collections. Journal of Mathematical Imaging
and Vision, 43:93–107, 2013. 2, 6, 7, 8
[9] Y. Furukawa, B. Curless, S. M. Seitz, and R. Szeliski. Towards internet-scale multi-view stereo. In Proc. of Computer
Vision and Pattern Recognition, 2010. 1, 2
[10] Y. Furukawa and J. Ponce. Accurate, dense, and robust multiview stereopsis. IEEE Transactions on Pattern Analysis andMachine Intelligence, 32(8):1362–1376, 2010. 1, 2, 5
[11] M. D. Grossberg and S. K. Nayar. Modeling the spaceof camera response functions. IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(10):1272–1282,
2004. 2, 4, 5[12] J. Hays and A. A. Efros. Scene completion using millions ofphotographs. ACM Transactions on Graphics (Proc. of ACM
SIGGRAPH), 26(3), 2007. 2[13] Y. Hwang, J.-S. Kim, and I. S. Kweon. Difference-based
image noise modeling using Skellam distribution. IEEE
Transactions on Pattern Analysis and Machine Intelligence,34(7):1329–1341, 2012. 6[14] I. Kemelmacher-Shlizerman, E. Shechtman, R. Garg, andS. M. Seitz. Exploring photobios. ACM Transactions on
Graphics (Proc. of ACM SIGGRAPH), 30(4):61:1–61:10,
2011. 2[15] S. J. Kim, J. Frahm, and M. Pollefeys. Radiometric calibration with illumination change for outdoor scene analysis. In
Proc. of Computer Vision and Pattern Recognition, 2008. 1,2, 4[16] S. J. Kim and M. Pollefeys. Radiometric self-alignment ofimage sequences. In Proc. of Computer Vision and Pattern
Recognition, 2004. 2, 4
[17] S. Kuthirummal, A. Agarwala, D. B. Goldman, and S. K.
Nayar. Priors for large photo collections and what they re-veal about cameras. In Proc. of European Conference on
Computer Vision, 2008. 2[18] P.-Y. Laffont, A. Bousseau, S. Paris, F. Durand, and G. Drettakis. Coherent intrinsic images from photo collections.
ACM Transactions on Graphics (Proc. of ACM SIGGRAPH),
31(6):202:1–202:11, 2012. 1, 2, 6
[19] J.-Y. Lee, Y. Matsushita, B. Shi, I. S. Kweon, and K. Ikeuchi.
Radiometric calibration by rank minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence,35(1):144–156, 2013. 2, 3, 4, 5, 6, 7, 8
[20] H. Li and P. Peers. Radiometric transfer: Example-based
radiometric linearization of photographs. Computer Graphics Forum (Proc. of Eurographics Symposium on Rendering),
34(4):109–118, 2015. 2[21] S. Lin, J. Gu, S. Yamazaki, and H.-Y. Shum. Radiometriccalibration from a single image. In Proc. of Computer Vision
and Pattern Recognition, 2004. 2, 7, 8
[22] S. Lin and L. Zhang. Determining the radiometric response
function from a single grayscale image. In Proc. of Computer
Vision and Pattern Recognition, 2005. 2[23] R. Martin-Brualla, D. Gallup, and S. M. Seitz. Time-lapse
mining from internet photos. ACM Transactions on Graphics(Proc. of ACM SIGGRAPH), 34(4):62:1–62:8, 2015. 2[24] Y. Matsushita and S. Lin. Radiometric calibration fromnoise distributions. In Proc. of Computer Vision and Pattern
Recognition, 2007. 2[25] T. Mitsunaga and S. K. Nayar. Radiometric self-calibration.
In Proc. of Computer Vision and Pattern Recognition, 1999.
1, 2, 3, 4[26] T.-T. Ng, S.-F. Chang, and M.-P. Tsui. Using geometry invariants for camera response function estimation. In Proc. ofComputer Vision and Pattern Recognition, 2007. 2[27] J. Park, Y.-W. Tai, S. N. Sinha, and I. S. Kweon. Efﬁcient androbust color consistency for community photo collections. In
Proc. of Computer Vision and Pattern Recognition, 2016. 2[28] L. Shen and P. Tan. Photometric stereo and weather estimation using internet images. In Proc. of Computer Vision and
Pattern Recognition, 2009. 2[29] B. Shi, K. Inose, Y. Matsushita, P. Tan, S.-K. Yeung, andK. Ikeuchi. Photometric stereo using internet images. In
Proc. of International Conference on 3D Vision, 2014. 1, 2,
8[30] B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan. Selfcalibrating photometric stereo. In Proc. of Computer Vision
and Pattern Recognition, 2010. 2[31] N. Snavely, R. Garg, S. M. Seitz, and R. Szeliski. Finding paths through the world’s photos. ACM Transactions
on Graphics (Proc. of ACM SIGGRAPH), 27(3):15:1–15:11,
2008. 2[32] N. Snavely, S. M. Seitz, and R. Szeliski. Photo tourism:
Exploring photo collections in 3D. ACM Transactions on
Graphics (Proc. of ACM SIGGRAPH), 25(3):835–846, 2006.
1, 2, 5[33] Y. Tsin, V. Ramesh, and T. Kanade. Statistical calibration of
CCD imaging process. In Proc. of Internatoinal Conference
on Computer Vision, 2001. 2[34] B. Wilburn, H. Xu, and Y. Matsushita. Radiometric calibration using temporal irradiance mixtures. In Proc. of Computer Vision and Pattern Recognition, 2008. 2