1

Radiometric Calibration by Rank Minimization
Joon-Young Lee, Student Member, IEEE, Yasuyuki Matsushita, Senior Member, IEEE, Boxin Shi,
In So Kweon, Member, IEEE, and Katsushi Ikeuchi, Fellow, IEEE
Abstract—We present a robust radiometric calibration framework that capitalizes on the transform invariant low-rank structure in
the various types of observations, such as sensor irradiances recorded from a static scene with different exposure times, or linear
structure of irradiance color mixtures around edges. We show that various radiometric calibration problems can be treated in a principledframework that uses a rank minimization approach. This framework provides a principled way of solving radiometric calibration problems
in various settings. The proposed approach is evaluated using both simulation and real-world datasets and shows superior performance
to previous approaches.
Index Terms—radiometric calibration, camera response function, rank minimization, low-rank structure

F

1

I NTRODUCTIONTHE radiometric response function relates sensor
irradiance and image brightness values. In mostcameras, there exists a non-linear radiometric response
function for purposes such as compressing the dynamic
range of sensor irradiance or adapting to a display’s
non-linear mappings. The non-linearity is intentionallydesigned by camera manufacturers for enhancing the visual quality of recorded images. However, for computer
vision applications, this non-linearity causes problems,
especially when a linear relationship between the sensor
irradiance and the recorded intensity is explicitly assumed. For example, photometric approaches like shapefrom-shading and photometric stereo methods generally
assume a linear (or affine) mapping from the physicalscene brightness to the image intensity. If the non-linear
mapping is not accounted for in such methods, it leads
to unexpected error in the computed results. Therefore,
radiometric calibration is an important first step forvarious computer vision algorithms that assume a linear
mapping to make them work correctly.
The problem of radiometric calibration is determining
the non-linear camera response function from observations. The observations can be multiple images takenwith different exposure times from a fixed camera location [1], [2], color distributions around image edges [3],or color profiles obtained from varying lightings [4].
In all of these cases, recovering radiometric response
functions without additional knowledge or assumptions
is an under-constrained problem. Since these methods
use different observations, they naturally use different• J.-Y. Lee and I. S. Kweon are with the Department of Electrical Engineering, KAIST, Korea.
E-mail: jylee@rcv.kaist.ac.kr, iskweon@ee.kaist.ac.kr• Y. Matsushita is with Microsoft Research Asia, Beijing, China.
Email: yasumat@microsoft.com• B. Shi and K. Ikeuchi are with Institute of Industrial Science, theUniversity of Tokyo, Japan.
Email: {shi, ki}@cvl.iis.u-tokyo.ac.jp

assumptions and criteria for recovering the response
functions. Our key observation is that a class of radiometric calibration methods rely on a certain linearity
defined in spaces that they use.
Based on this observation, we formulate the radiometric calibration problem in a unified and principledframework. Our method capitalizes on the fundamental
linear dependency of sensor irradiances. Our methodarranges the observations (e.g., multiple images obtainedwith different exposure times, or multiple images withvarying illumination, and etc.) in a matrix form, wherea column vector corresponds to an observation vector.
When the linear relationship is observed in such amatrix, the rank of the matrix ideally becomes one. On
the other hand, when the response function is nonlinear, the matrix generally becomes higher or even fullrank. From this observation, we develop a framework
for radiometric calibration problems using a matrix rank
minimization technique (see Fig. 1).
There are a few benefits from casting the problem inthis manner. First of all, the rank minimization approachcan efficiently avoid the over-fitting problem. While
there exist robust approaches [5] that are not affected by
the over-fitting problem, still the problem is commonly
observed in previous approaches [2] that rely on the`2 -norm minimization. Second, it gives a unified framework for solving the radiometric calibration problem by
explicitly using the linear dependency of the irradiancevectors. With this framework, various radiometric calibration problems can be treated in the same manner, andtherefore gives a better understanding to the problem.
Third, because of the new formulation, we can utilize anadvanced rank minimization technique that can reliably
derive the solution in a principled manner.
In presenting our work, we begin by briefly describing
the radiometric response function and reviewing related
works.

2

0.02

1

0.01

0.8

0.3
0.1

0.6
0
0.4
-0.01

Imaging process

-0.02

-0.1

0.2

image noise0

non-linear response function-0.3

outliers

irradianceoobservation

1

1

1

1

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0

0

0

0

Our method

rank minimization

outlier rejection

Fig. 1. Illustration of our approach. Our method arranges the observations (e.g., multiple images obtained with different
exposure times, or color profiles obtained from image edges, etc.) in a matrix form, where a column vector corresponds
to an observation vector. In the matrix form, the irradiance matrix should be rank-1, while the observations matrix
generally has a higher rank. Our method seeks the inverse response function that transforms the observation matrix
into a low-rank matrix.
1.1

Radiometric response function

The radiometric response function relates the sensor irradiance I and recorded intensity value M (observation)by
M = f (I).

(1)

In the case of color cameras, color channel has an independent mapping. With assumptions of monotonicity
and continuity, the response function f : R → R canbe uniquely inverted to an inverse response functiong = f −1 if we neglect quantization and clipping error.
Using the inverse response function g, the recorded
intensity values can be linearized by g(M ). Since whatis available as input are the observations M , most radiometric calibration methods aim at recovering the inverse
response function g. Once we know the inverse response
function g, we can automatically determine the response
function f with the aforementioned assumptions.
1.2

Prior work

There are various approaches for radiometric calibration.
One widely-used approach takes as input a set of images
recorded with varying exposures from a fixed camera.
The early work by Mann and Picard [6] and Debevec
and Malik [1] show that the response function can berecovered using images acquired by a static camera butwith known exposure ratios. Mann and Picard [6] use agamma correcting function to represent response functions. Debevec and Malik [1] assume a smoothness property of response functions and estimate them in a nonparametric manner. Mitsunaga and Nayar’s method [2]uses a parametric yet flexible representation, i.e., a polynomial representation. Using approximate prior knowledge of the ratios of exposures, their approach iteratively

estimates response functions and accurate ratios of the
exposures. Grossberg and Nayar [7] create a database ofreal-world camera response functions (DoRF) and use it
for better representation of response functions. Instead oftaking a set of images, Nayar and Mitsunaga [8] use an
optical mask to capture an image with spatially varying
pixel exposures.
There are a few prior approaches that allow somecamera movement or scene motion. In Mann and Mann’s
method [9], images are taken by a rotating and zooming
camera and the method simultaneously estimates the
response function with a non-parametric model, as wellas the exposure ratios. Grossberg and Nayar [10] use intensity histograms of two images of different exposuresinstead of finding exact pixel correspondence for determining response functions from image sequences withcamera/scene movement. Litvinov and Schechner [11]propose an approach for simultaneously estimating theradiometric response, spatial non-uniformity, and gain
variations from an image sequence acquired during
camera motion. Kim and Pollefeys’s method [5] allows
for free movement of the camera and some motionin the scene by dynamic programming on the joint
histogram built from correspondences. As pointed out
in [5], [10], [11], without the knowledge of exposure
ratios, the estimate still has an exponential ambiguity.
While not unique, such an estimate is still useful formany applications, such as radiometric alignment, highdynamic range (HDR) image production, and smooth
image stitching.
Instead of using varying exposures, some approaches
use statistical properties embedded in images to achieve
radiometric calibration. Farid [12] presents a blind
gamma estimation by minimizing specific higher-ordercorrelations in the frequency domain. Tsin et al.’s

3method [13] estimates non-parametric response functions using a statistical model of the CCD imaging
process. Pal et al. [14] use probabilistic imaging models
and weak prior models for deriving response functions
to produce high-quality HDR images. Ng et al. [15]present a geometry invariance-based method with ageneralized gamma curve representation for estimating
response functions. Matsushita and Lin [16] propose to
use the symmetric property of image noise by observing
noise distributions contained in images. Takamatsu etal. [17] improve the noise-based method with a probabilistic intensity similarity measure, which requires afewer number of images. In [18], Takamatsu et al. derive
the relationship between the response function and noise
variances, and estimate response functions using theimage noise variance. In a more general setting, Lin etal. [3] and Lin and Zhang [19] propose methods that takeonly a single image as input. Their methods use edges
for obtaining color or gray-scale histogram distributions,and the optimal inverse response function is determinedby transforming the observed non-linear distributions
into linear distributions. A similar idea is introducedby Wilburn et al. [20]. Wilburn et al. [20] use temporal
irradiance mixtures caused by motion blur, instead of
spatial mixtures at edges.
There are some other methods which use varying illumination conditions [4], [21], [22], [23]. Manders et al. [21]use the superposition property of light, and Shafiqueand Shah [22] use illumination and geometric invariants.
Kim et al. [23] present a method to use natural lighting
conditions to obtain exposure changes. Shi et al. [4]perform calibration by linearizing a set of color values
at the same coordinate from images taken under varying
lighting conditions in a photometric stereo setting.
The idea of utilizing low-rank structures in real-world
scenes has been used in various other vision problems.
For example, the low-rank representation is shown to
be useful for face recognition [24], structure from motion [25], and photometric stereo [26]. In more recent
works, low-rank structure is widely explored in various novel applications with an emphasis on advanced
optimization techniques. Peng et al. [27] align multiple
images of objects of interest to a fixed canonical template
to make the measurement on image similarity moreaccurate for recognition or classification purposes. Thesealigned images naturally bear the low-rank property.
Paladini et al. [28] propose a variant shape model called a
3D-implicit low-rank shape model to represent the shapeand apply it to non-rigid structure from motion. Zhang etal. [29] propose a transform invariant low-rank texture
to deal with images containing camera projection deformations such as affine or projective transformations.
Wu et al. [30] robustly solve photometric stereo in thepresence of large and sparse errors such as specularity
by extracting low-rank structure from the input matrix.
Zhang et al. [31] present a method for estimating acamera’s intrinsic and extrinsic parameters and lensdistortion from low-rank textures that appear in images.

All these current research works show the variety andvalidity of low-rank structures that exist in the realworld.
A preliminary version of this paper that focuses onradiometric calibration from multiple-exposure images
appeared in [32]. We extend the work by defining aunified framework that treats a class of radiometric calibration problems in the same manner. With the proposed
framework, we show that outliers can be efficiently
handled with a simple optimization algorithm. Also, we
perform extensive experiments with carefully redesigned
simulations and in-depth analysis for analyzing the behavior of the proposed method.

2 L OW- RANK S TRUCTURE
O BSERVATIONSINC ALIBRATEDOur method casts the radiometric calibration problem
as a low-rank structure recovery problem. As describedin the previous section, there are various approaches forradiometric calibration depending on the type of observations, e.g., a set of images with varying exposures, a setof images with varying illuminations, or a single image
that contains irradiance mixtures around edges. We firstshow that these radiometric calibration problems that
use the linearity in sensor irradiance can be formulated
as a rank minimization problem.
Suppose we are given a set of observation vectors{vec(M1 ), . . . , vec(Mn )} where each observation vector
is defined as a vector of measured intensities vec(Mi ) ∈Rm×1 . The numbers m and n represent the number of
elements in each observation vector vec(Mi ) and thenumber of observation vectors, respectively. We create
an observation matrix D ∈ Rm×n by stacking all theobservation vectors vec(Mi ) as
D = [vec(M1 ) | · · · | vec(Mn )] .

(2)

If we know the inverse response function g(= f −1 ),the observation matrix D can be transformed into an
irradiance matrix A with a scaling ambiguity byg ◦ D = A = [vec(I1 ) | · · · | vec(In )] ,(3)

where ◦ is an operator that describes the element-wise
mapping, and vec(Ii ) ∈ Rm×1 represents an irradiancevector. For the problems that we discuss, because theirradiance vectors vec(Ii ) should be linearly dependent,
the rank of the irradiance matrix A becomes one. Therefore, the problem of radiometric calibration can be formulated as determining the inverse response function
g that minimizes the rank of the transformed matrixg ◦ D(= A):ĝ = argmin rank (A)s.t.

A = g ◦ D.

(4)g

This formulation is fairly general and capable of treating
various types of observation vectors in a unified framework. The only requirement is to arrange the observation
vectors in a way that they exhibit linear-dependency

4after transformation by the correct inverse response
function g. In the following, we describe the usage of thisframework with different types of observation vectors.
2.1 Multiple-exposure case
One of the most widely used radiometric calibration
methods uses multiple images of a static scene takenfrom a fixed camera with various exposures, such as Mitsunaga and Nayar’s method [2]. In this setting, because
the scene radiance L is constant, the sensor irradiance I
becomes proportional to the exposure time e as
(5)

I = kLe ∝ e,with a constant scaling k. For each exposure, an observation vector vec(Mi ) is created by vectorizing the image ofm pixels, e.g., by arranging recorded intensities in a scanline order. Given n different exposures, the observation
matrix D can be created by stacking n observation
vectors as Eq. (2). While we cannot directly observethe (vectorized) sensor irradiance vec(Ii ), we know thatthe irradiance is proportional to the exposure time from
Eq. (5). In other words, we know the irradiance matrix
A = [vec(I1 ) | · · · | vec(In )] becomes rank-1. Therefore,
we can directly use the rank minimization formulationof Eq. (4) for the radiometric calibration problem inthis setting, i.e., finding the inverse response function gthat best minimizes the rank of transformed observation
matrix g ◦ D.
2.2 Varying illumination case
Recently, Shi et al. [4] showed that radiometric calibrationfrom photometric stereo images can be performed bylinearizing the color profiles. The color profile is defined
as an ordered set of color values in the RGB color
space. Their method capitalizes on the property thatcolor profiles form straight lines in the RGB color space if
the response function is linear, while non-linear response
functions f bend color profiles to be non-linear. Theirmethod therefore seeks the inverse response function gthat linearizes all the color profiles that are observed.
It is straightforward to cast their problem into our lowrank computation scheme. By treating each color profile
as an observation matrix, we can cast the linearization
problem as a matrix rank minimization problem. The
color profile matrix should become rank-1 when thecorrect inverse response function is applied.
Given images captured under n different lighting conditions, we can form a color profile matrix D ∈ R3×nfrom each pixel location. The number of matrix rows
corresponds to the number of color channels (RGB). For
each of m pixels we have such a matrix Di , i = 1, . . . , m.
If a correct inverse response function g is found, ittransforms the observation matrix Di to an irradiancematrix Ai , where the rank of Ai is one. Using all the
pixels, we have an objective function written as:ĝ = argmin
g

m
X
i=1

rank (Ai )s.t.

Ai = g ◦ Di .

(6)

Notice that it computes a single global inverse response
function g by minimizing the rank of all the color profile
matrices. Thus, the radiometric calibration problem inthis setting can be solved in the same computationmanner as discussed in the previous section.
2.3

Color mixture case

There are radiometric calibration approaches that use
color/irradiance mixtures in an observed image. Lin etal. use spatial color mixtures around edges [3], andWilburn et al. use temporal irradiance mixtures of motion
blur [20]. The color/irradiance mixtures are observed
around image edges where two different radiances are
mixed together due to the limited spatial resolution of
the image array or temporal motion blur. In edge regions, RGB colors are linearly mixed together, however,observed intensities have non-linearity due to non-linear
response functions. These methods use this property for
estimating the inverse response function by transforming
non-linear edge color distributions to linear ones.
Our method can also be applied to this setting in asimilar manner with the previous varying-illumination
setting. The difference is that a color mixture forms aline passing through two different color values in the
RGB space, while a color profile always passes throughthe origin. To apply our rank-minimization approach, wetranslate the color mixture profile so as to pass throughthe origin by simply subtracting the RGB value of either
color from the color mixture profile. This operation is
to avoid offsets in the data; in this manner, we convert
the affine problem into a linear problem. After the coordinate transform of the observation matrix, the exactly
same formulation as Eq. (6) can be applied in this setting.

3

C ALIBRATION A LGORITHM

To efficiently compute the rank minimization problem, previous approaches use an approximate solution
method that minimizes a nuclear norm [33], which is.
definedas the sum of the singular values ||A||∗ =Pn
i=1 σi (A). The nuclear norm minimization turned out
to be an effective computation tool that can be applied
to problems when the transformation does not alter themagnitude of matrix elements.
In our case, however, nuclear norm minimization cannot be directly employed because the function g significantly alters the absolute values of matrix elements,
and it therefore results in variations of the singularvalue magnitudes. When nuclear norm minimization is
directly applied in this setting, the optimization prefers adegenerate solution where the inverse response function
g tries to turn matrix elements into zero so that the
singular values become smaller. To resolve this issue, we
instead use condition numbers κi , i.e., a ratio of singularvalues, defined as.
κi (A) = σi (A)/σ1 (A), i = 2, 3, . . . , n.

(7)

g

1

1

0.8

0.8

0.6

0.4

0.2

0

0.6

0.4

0.2

0
0

0.2
0.4
0.6
0.8
Normalized measurement1

Normalized irradiance10.8

Normalized irradiance10.8

Normalized irradianceWith the condition numbers κi , we can work with therelative magnitudes of the singular values that are unaffected by their absolute magnitudes. Now we approximate the energy functional (4) as a minimization of thesummation of condition numbers as:nXĝ = argmin
κi (A) s.t. A = g ◦ D.
(8)

Normalized
irradiance
Normalized irradiance

5

0.6

0.4

0.2

0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
Normalized measurement
Normalized measurement
Normalized observation

0

(a)

(b)

0.6

0.4

0.2

0
1

0

0.2
0.4
0.6
0.8
Normalized measurement

(c)

1

(d)

i=2

There are two main factors causing rank variations:non-linearity of response function and image noise. To
see these effects, let us take synthetic response functions(RF 1 − 4) shown in Fig. 2 as an example. The top halfof Table 1 shows condition numbers κ2 , . . . , κ5 obtained
from synthetic data that have five different exposure
times for each response function. κall indicate the sumof condition numbers. Depending on the shape of the
response function, condition numbers vary in the table. With the ideal linear response function RF 1, allcondition numbers become zero, which indicates rank-1
structure. With other more general non-linear response
functions, however, the condition number κ2 becomessignificantly larger than other condition numbers. Thisis commonly observed because of the smooth (or lowfrequency) characteristics of response functions.
In addition, the effect of noise also changes the condition numbers. To investigate of the effect of noise, we
add image noise to the irradiance using a Poisson noisemodel [34]. In [34], the Poisson noise ηe is modeled as.
ηe = Ie − µ s.t.

µ(Ie ) e−µ,
Ie ∼ fP (Ie |µ) =(Ie )!

(9)

where Ie is a noisy observation that follows a Poisson
distribution with a variance µ. The Poisson distribution
is shifted by µ, therefore the mean of image noise ηe iszero. The variance µ is determined for each irradiance I
asµ(I) = RI I + DI ,(10)
where RI is a scalar parameter for signal-dependent shotnoise, and DI refers to dark current noise. We represent
both parameters RI and DI as functions of the camera
gain Cg and empirically determine them as
RI = 0.035Cg ,DI = 0.1Cg2 ,(11)using images captured from a Point Grey Scorpion camera for this test. The response functions are applied after
adding the noise with the camera gain Cg = 3 to theirradiances. These signals are quantized at 256 levels,and then the condition numbers are re-computed. The
average condition numbers of 100 trials are shown inthe bottom half of Table 1. The image noise makes thecondition numbers similar across diverse response functions. It becomes more obvious in the higher condition
numbers, e.g., κ3 ...κ5 , because the noise almost evenly
affects all the condition numbers while the signals in
the higher condition numbers have smaller values. In
the absence of noise, we can achieve the same result
from minimizing κall and minimizing κ2 . However, inFig. 2. Reference response functions used for Table 1.
(a) Ideal linear response function (RF 1). (b-d) Non-linear
response functions (RF 2 − 4).
condition number
RF 1
without
RF 2
noise
RF 3
RF 4
RF 1
with
RF 2
noise
RF 3
RF 4

κ2
0.00000.02840.11710.00280.01780.03160.11810.0199κ30.00000.00350.01860.00000.01270.01570.02580.0107κ40.00000.00090.00250.00000.00970.01400.01500.0065κ50.00000.00040.00000.00000.00780.01230.00580.0044

κall
0.00000.03320.13830.00280.04800.07370.1646
0.0416TABLE 1
Variations of condition numbers w.r.t. response functions
illustrated in Fig. 2. Top half shows the condition
numbers of the response functions RF 1 − 4. Bottom halfshows the case when noise is added.

the presence of noise, signal-to-noise ratios of the higher
condition numbers rapidly degrade, therefore minimizing κall as Eq. (8) may cause an over-fitting result. For
example, it shows that κall of RF 4 is smaller than κall of
RF 1 in Table 1 with noise. From these observations, we
use only the condition number κ2 instead of using all
the condition numbers. Therefore, the objective function
is simplified from Eq. (8) toĝ = argmin κ2 (A)s.t.

A = g ◦ D.

(12)

g

This objective function is different from recent sparseand low-rank decompositions [27], [29], [30], which are
based on the robust principal component analysis (PCA)
algorithm [33]. The goal of sparse and low-rank decomposition is to find the minimum rank automaticallyin the presence of sparse outliers, while our objectivefunction finds the best rank-1 matrix in the presenceof noise that is not necessarily sparse. In other words,
we explicitly use the prior knowledge about the rank1 structure for efficiently solving the problem in thissetting.
3.1

Monotonicity constraint

In general, response functions as well as inverse response functions are monotonic functions. To avoid trivial solutions, we enforce this property during the computation by putting a monotonicity constraint, which isrepresented as ∂g/∂D > 0. The constraint can be implemented as either hard or soft constraints. We empirically6employ a soft constraint because it can be well solved bya simple optimization method. With this constraint, thefinal form of the objective function becomesĝ = argmin κ2 (A) + λ

X

g

t



∂g(t)
H −
∂D
s.t.

A = g ◦ D, (13)where H(·) is the Heaviside step function (H(x) = 1
when x ≥ 0, and H(x) = 0, otherwise). The derivatives∂g(t)/∂D are assessed by sampling g(t) at various t in
the range of [0, 1].
3.2

Handling outliers

While we see that our method is robust against highfrequency noise in the previous section, systematic outliers can still cause drifting during the optimization. In
radiometric calibration, there exist various outliers due
to misalignment of images, small motions in a scene,
shadows, and so on. We develop a method based on
robust PCA [33] for reliably rejecting these outliers.
Robust PCA decomposes the input data matrix X as
X = Z + S,

(14)

where Z is a low-rank matrix, and S is a sparse error
matrix. The solution is obtained by solving
minimize||Z||∗ + λs ||S||1

subject to

Z + S = X,

(15)

where λs is a regularization weight.
However, straightforward application of the robust
PCA algorithm is unsuitable for our data because ourdata have different characteristics in the matrix thanwhat robust PCA algorithms address. Because of therank-1 structure of irradiance matrices and the lowfrequency nature of response functions, our data has one
dominant singular value. This introduces confusion torobust PCA in determining the rank of our data, therefore it is difficult to set a proper λs for balancing thesetwo costs. A larger λs favors more outliers remaining
in the low-rank matrix Z, while a smaller λs prohibits
rank variations, and even non-linear response functions
can be treated as an outlier.
For these reasons, we make a simple modification
to the original robust PCA for removing outliers. We
iteratively decompose the matrix X into rank-1 andhigh-rank matrices, and apply thresholding to the highrank matrix to find outliers according to the statisticsof the high-rank matrix. The statistics are calculatedwith an assumption of the Gaussian distribution of rank
deviations from the rank-1 approximation.
We show the outlier rejection process in Algorithm 1.
In the algorithm table, R is a rank-1 matrix, N is ahigh-rank matrix, σi is the ith singular value, ui and vi
are ith columns of U and V , respectively, and Ψ is a

Algorithm 1 Outlier Rejection Process
1: procedure O UTLIER R EJECTION (X, ρ)
2:
Z0 = X
3:
while not converged do
4:(U, Σ, V ) = svd(Zk )5:
R = U Σ1 V T = σ1 u1 v1T
6:
N = Zk − R
7:dE = Ψρ [N ]8:
Zk+1 = Zk − dE
9:end while
10:
return Z ∗11: end procedurethresholding operator, defined as(ni if |ni − µN | >  + ρsNΨρ [ni ] =0 otherwise,(16)

where  is small epsilon value, and µN and sN are
the mean and standard deviation of N , respectively.
The operator is extended to matrices by applying it toeach element. Therefore, only outlier components out ofρ standard deviations remain in N .
3.3

Ambiguities

The estimated inverse response function ĝ may sufferfrom ambiguities that are specific to the calibration setting. For example, a certain class of response functions
f where f (x)/f (ax) = const. holds cannot be handled in
the varying illumination setting [4]. In the color mixture
setting, we cannot determine response functions whenboth measured color distribution lies along the R=G=B
line and RGB channel responses are identical [3]. In
practice, these cases are rare while theoretically they
exist.
One of the major ambiguities is exponential ambiguityin the multiple-exposure setting discussed by Grossberg
and Nayar [10], i.e., if g is a solution for I = g(M ),then g γ can also become a solution as I γ = g γ (M )for arbitrary γ. Recovery of the exposure ratios and the
response function is impossible without either making
assumptions on the response function or having roughestimates on the exposure ratios.
Without making any assumptions about exposure
times, our method can robustly recover the response
function up to the exponential ambiguity from at least
two images. However, our method also suffers fromthis exponential ambiguity like other methods when the
exposure times are completely unknown, because therank-1 structure is retained even after any exponential function is applied. In short, if rank(I) = 1, thenrank(I γ ) = 1 for any γ.
To resolve the ambiguity, we use only one known
exposure ratio. Using the known exposure ratio, we canestimate γ by solving the least-square problem described7asγ̂ = argmin
γ

X

γ

γ

2

[ĝ (Mi ) − ri,j ĝ (Mj )] ,

(17)

i,j

where ri,j is the exposure ratio ei /ej of measurement
pairs Mi and Mj .
3.4

Representation of response functions

The empirical model of response function (EMoR) [7]and polynomial representation [2] are the most widelyused representations of response functions. As mentioned in various prior works [3], [5], [16], [17], [18], [19],[23], EMoR can represent various response functions
with few coefficients. However, it is difficult to computederivatives for EMoR either analytically or numericallybecause EMoR is a highly non-smooth function. Therefore, EMoR is not suitable for gradient-based convex
optimization. On the other hand, the polynomial representation is more appropriate for gradient-based convex
optimization because of its smoothness. Also, recent
research shows that most typical response functions canbe represented well with only a 5th order polynomial
function [35]. For these reasons, we choose to use thepolynomial representation as done in [2].
We normalize the irradiance and observation domains
in the range of [0, 1] by enforcing g(0) = 0 and g(1) = 1.
We explicitly embed this boundary condition in theparametric representation of g, and the polynomial representation becomesg(M ) = M + M (M − 1)n−1
X

ci M n−i−1 ,(18)where n is the order of a polynomial function, and ci
are the coefficients to estimate. With this, the degree of
freedom is reduced to only 4 for the 5-order polynomial
representation.

p(g) =KX

αi N (g; µi , Σi ).

(19)

i=1

The likelihood of the rank constraint is modeled asp(D|g) =1
exp (−λE (g; D)) ,Z(20)

where Z is a normalization constant, λ is a regularization
parameter, and E is a rank error function. The optimal
response function ĝ is obtained by solving the following
MAP problem:ĝ = argmax p (g; D) = argmax p (D|g) p (g) .

(21)

g

Substituting Eq. (19) and Eq. (20) into the logarithmic
form of Eq. (21), we obtainĝ = argmin λE (g; D) − log p(g).

(22)

g

The use of the prior information is demonstrated in Sec.
4.3.

Iterative optimization

We use the rank minimization and outlier rejection algorithms in an iterative manner to solve the problem. The
overall process is summarized in Algorithm 2. Once theoptimal solution of inverse response function g is obtained, if exponential ambiguity exists in the solution in
the multiple-exposure setting, we resolve the exponential
ambiguity γ after convergence using Eq. (17).
3.6

Using an Expectation Maximization (EM) algorithm, theprior model is constructed with K kernels as

g

i=1

3.5

Algorithm 2 Radiometric Calibration Process
1: procedure R ANK C ALIB (D)2:initialize g3:while not converged do
4:
A=g◦D
5:
Z = OutlierRejection(A, ρ)
6:
Calculate error using Eq. (13)7:
Update g
8:
end while
9:
return ĝ
10: end procedure

Prior model

In most cases, we can estimate a right solution from
Eq. (13). However certain methods like color mixturecase [3] have a difficulty to estimate a unique solution
due to insufficient information. For that case, we canuse prior information to regularize the solution. The
prior model is formed by fitting a multivariate Gaussian
mixture model to the DoRF database [7] as done in [3].

4E XPERIMENTS

We evaluate the proposed method in three different
settings: multiple-exposure, varying illumination, andcolor mixture settings. In particular, we make a thoroughanalysis of the proposed method using the multipleexposure scenario in Sec. 4.1. Sec. 4.2 and Sec. 4.3
describe the result of varying illumination and color
mixture cases, respectively.
The optimization of Eq. (13) can be performed using
various minimizers. Though Eq. (13) is not convex, it
converges to the correct solution with a broad range of
input data as we will see below. In our experiments, we
use the Levenberg-Marquardt (LM) method. We always
set the weighting factor λ = 1e10 in Eq. (13) and ρ = 3
in Algorithm 1 throughout the experiments.

8

All datasets

D3

D4

Likelihood

All datasets

3000

3000

2500

2500

2000

2000

1500

1000

MN − 3rd# of scenes < RMSE

D2

# of scenes < RMSE

D1

1500

1000

Ours − 3rd

MN − 4th

Ours − 4th

MN − 5th
th

500

0

0.2

0.4

0.6

0.8

MN − 6Ours − 6th

MN − 7th

1

Normalized radiance

Ours − 5th500

Ours − 7th

MN − 8th0

Ours − 8th
0

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Fig. 3. Scene radiance distributions of four synthetic
datasets (D1-D4).

Fig. 4. Effect of numbers of coefficients of the polynomial
representation.

4.1

simply starting from a linear initial guess, and the otheruses five different gamma curves as initial guesses of
the inverse response function g and picks the optimal
solution. In the latter setting, the range of γ is set to0.8−1.6 with a step of 0.2. In the experiment, we observethat our method as well as the MN method are notsusceptible to the initial guess. The solution using ourrank minimization approach may result in various local
minima due to the exponential ambiguity. However, forradiometric calibration purposes, it is sufficient to findsuch a local minimum for our method because it is oneof the true solutions with the exponential ambiguity
in most cases. We use the linear response function as
an initial guess for the rest of the evaluations for its
simplicity.
Number of coefficients. One of the most importantparameters for both methods is the number of coefficients in representing inverse response functions. In thepolynomial representation, it is the degree of polynomial
functions. To find the best number of coefficients, we
perform performance evaluation with different numbersof coefficients, varying from 3rd to 8th -order polynomials.
Fig. 4 shows performance variations w.r.t. the polynomial orders for the MN method and our method. The
MN method always produces a smaller cost value whenhigher-order polynomials are used. However, it does
not necessarily mean a better estimation. Instead, it isdue to an over-fitting problem, and it leads to worse
estimates. We observed that for the MN method a 5thorder polynomial works well with various noise levels.
This result agrees with the observation in [35], where
it is said that a 5th order polynomial representation
has a good balance between accuracy and flexibility inrepresenting response functions.
On the other hand, normally our method shows better
results as the polynomial order increases. In our method,
the smaller cost values coincide well with the precisionof estimates. We decide to use a 6th order polynomial
for our method because we observe that there is no
performance gain from a polynomial order higher than
the 6th order.
Condition numbers and image noise. In Sec. 3, we

Multiple-exposure caseTo evaluate the proposed method, we perform experiments using both simulation and real-world data. We
compare our results with results of Mitsunaga andNayar’s method [2] (MN method, henceforth). In our
implementation, we use a polynomial function to represent response functions as described in Sec. 3.4 as donein [2]. Once the response function is estimated up to
the exponential ambiguity, we resolve the ambiguity todetermine the final response function using Eq. (17). This
optimization is also performed using the LM method.
4.1.1 Simulation
In this section, we use a synthetic dataset for quantitativeevaluation. The synthetic dataset is generated using
the DoRF database [36], which contains 201 measured
response functions. We synthetically generate scene radiances in the range of [0, 1] and create 1000 observations(which forms a 1000-pixel image) using each response
function with five exposure times. The step of exposure
ratios is set to 0.5 because this is a commonly available
setting in commercial cameras. Fig. 3 shows four scene
radiance distributions (D1-D4) of the synthetic datasetthat we have used. We further add five different magnitudes of noise for each radiance distribution. Therefore,
the comparison with the MN method is performed using
4020(= 201 × 4 × 5) synthetic data in total. Poissonnoise distributions in Eq. (9) with camera gains Cg (={0, 1, 3, 6, 9}) are used for diversifying the datasets.
In presenting the simulation result, we plot cumulative
histograms that show the number of results below acertain root mean squared error (RMSE) for each configuration. The greater area under the curve shows thesuperior performance. On the top of each figure, wedenote the radiance distribution used for the experiment.
“All dataset” indicates that all the radiance distributions
are used.
Initial guess. We first assess the effect of initial guesses
to the solution. As mentioned, our method does nothave a guarantee to converge to the globally optimal
solution, like any other non-linear optimization problems. For this experiment, we use two settings; one is9D1

All datasets

All datasets

3000
MN − 0%
MN − 1%
MN − 3%RANK − 0%
RANK − 1%
RANK − 3%
Ours − 0%
Ours − 1%
Ours − 3%

2500

2500

2000

2000

80

κall, Cg=160

κall, Cg=6

κall, Cg=3
κall, Cg=9

# of scenes < RMSE

100

# of scenes < RMSE

# of scenes < RMSE

120

1500

1000

1500

1000

2000

1500

1000

κ2, Cg=140

κ2, Cg=3

MN
κall

500

κ2, Cg=620

κ2, Cg=902500# of scenes < RMSE

140

All datasets

3000

3000
160

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

500

Ours − 0% outliers
Ours − 1% outliers
Ours − 3% outliers
Ours − 5% outliers
Ours − 7% outliers

500

κ2
0

0

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Fig. 5. (Left) Evaluation of cost functions of Eq. (8)and Eq. (12) using synthetic scenes D1. “κall ” indicatesminimization of all the condition numbers using Eq. (8),while “κ2 ” corresponds to minimization of 2nd condition
number using Eq. (12). (Right) Performance comparisonof the MN method, Eq. (8) and Eq. (12)observed the relation between condition numbers andimage noise. Here, we perform experiments for further
validation. For this experiment, we compare the resultsof two cost functions, Eq. (8) and Eq. (12). Eq. (8)minimizes the sum of all condition numbers, while
Eq. (12) minimizes only the 2nd condition number. The
experimental results are shown in Fig. 5. In the left plot,
it shows that minimizing the 2nd condition number ismuch more robust than minimizing the sum of all condition numbers, especially when there is greater noise.
The right figure shows the cumulative histogram of the
MN method, Eq. (8) and Eq. (12) for all of datasets. In
this figure, the result of Eq. (8) shows worse result than
the result of the MN method, while Eq. (12) shows thebest result among the others.
Outliers. In practice, recorded data have outliers due
to misalignment of images or small motion in a scene.
To validate our outlier rejection procedure, we perform
experiments by adding outliers in the dataset. The outliers are added by randomly changing selected pixel
values in the range of [0, 1] of the original dataset. Fig. 6
shows results. While our method with outlier rejection
procedure is not affected by outliers, performances ofother methods, i.e., straightforward rank minimizationand the MN method, are highly degraded even for only
one percent outliers.
Results with synthetic dataset. Now we assess theperformance of our method in comparison with the
MN method. Fig. 7 shows results with the synthetic
dataset. With very small noise, the MN method works
better than ours. However, the performance of the MN
method rapidly degrades as the noise becomes larger
while our method is still robust against the greater noiselevel. As shown in the figure, the estimation accuracy
depends on not only the noise level but also the radiancedistribution. With uniformly distributed scene radiance,
the estimation becomes quite stable for both methods.
However, when the radiance distribution becomes bi-

0
0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Fig. 6. (Left) Effect of outliers on different methods. MN
indicates the MN method, RANK is a straightforward rank
minimization without outlier rejection, and Ours is ourcomplete method. (Right) Performance comparison of our
method with different amounts of outliers.
ased, which is very common, it introduces performancedegradation particularly for the MN method as shownin the D4 result. The biased radiance distribution introduces biased observations, thus it yields an over-fitting
problem in the MN method. In contrast, our method
suffers less from such an effect.
We summarize the quantitative results in Table 2. The
DoRF database contains some response functions which
are difficult to represent with 5th order polynomials.
These functions give trivial results for the MN method;
therefore, we only use the best 150 fitting results for eachof the algorithms to compute the mean and standard deviation of RMSE and disparity (the maximum deviationfrom the ground truth). The table shows that our method
performs well for most datasets. It also shows that our
method is robust against increasing noise levels, while
the performance of the MN method rapidly degrades.
4.1.2

Real-world experiment

We also perform experiments using real-world cameras
and scenes. In our experiments, we use three different
cameras: Canon 550D, Nikon D1X, and Sony DSC-828.
Each dataset is collected by capturing a static scene withdifferent exposure times. Fig. 8 shows estimation resultswith real-world scenes. For presenting the results, weconvert images to irradiances as done in [37] and compare the results both with and without calibration. We
estimate response functions from a set of input images
and linearize the input by applying estimated inverse
response functions. For the estimation of response functions, we first reject observations in the under-exposurelevel (e.g., 5 in 0-255 range) and saturation level (e.g., 250
in 0-255 range) to avoid estimation bias. There are whitebalancing and color space transformation to convert
camera color space into the linear sRGB. The transformation that includes both the white balancing and the color
space transformation can be represented by a 3×3 matrix
and can be obtained in the least-square manner. The
error maps are calculated by comparing the estimated

10

D1

D2

D3

D4

160

160

160

160

140

140

140

140

120

120

120

120

MN, Cg=1
MN, Cg=3
MN, Cg=6

MN, Cg=180

MN, Cg=3
MN, Cg=660

100
MN, Cg=180

MN, Cg=3
MN, Cg=660

MN, Cg=9

Ours, Cg=6

0

0.005

0.01

0.015
RMSE

0.02

0.025

MN, Cg=3
MN, Cg=660

0.03

Ours, Cg=1

40

Ours, Cg=3

Ours, Cg=620

Ours, Cg=9
0

0

0.005

0.01

0.015
RMSE

0.02

0.025

Ours, Cg=3
Ours, Cg=6100

Ours, Cg=9
80
60
40

Ours, Cg=3

Ours, Cg=620

Ours, Cg=90

MN, Cg=1

80

Ours, Cg=1

MN, Cg=9

Ours, Cg=1

40

Ours, Cg=320

100

MN, Cg=9

Ours, Cg=1

40

# of scenes < RMSE

100

# of scenes < RMSE

# of scenes < RMSE

# of scenes < RMSE

MN, Cg=920

Ours, Cg=9
0

0.03

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

0

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Fig. 7. Cumulative histogram of the number of successful cases w.r.t. RMSE.
Noise level (Cg )
MN
RMSE
Ours
D1
MN
Disparity
Ours
MN
RMSE
Ours
D2
MN
Disparity
Ours
MN
RMSE
Ours
D3
MN
Disparity
Ours
MN
RMSE
Ours
D4
MN
Disparity
Ours

00.00580.00710.01420.01570.00790.00950.01800.01960.00760.00800.01700.01730.00860.01440.01950.028510.00540.00700.01390.01640.00770.01000.01840.02140.00740.00860.01760.01840.01180.01310.02620.0275Mean
30.00620.00750.01520.01740.01340.01050.02990.02280.01110.01050.02500.02180.02330.01380.05070.029360.01090.00840.02380.01940.02080.0111
0.04680.02400.01740.01160.03980.02450.03930.01540.08220.032790.01600.01020.03540.02250.02870.0132
0.06340.02850.02340.01280.05210.0272
0.05650.01890.11550.039300.00540.00540.01340.01240.00760.00820.01720.01630.00720.00660.01640.01420.00790.01380.01840.0270Standard Deviation
1
3
60.00570.00560.00630.00600.00670.00670.01400.01390.01470.01430.01530.01640.00810.01070.01330.00860.00830.00940.01900.02390.03040.01830.01860.02070.00800.00920.01170.00690.00680.00860.01800.02070.02820.01390.01370.01820.00970.01520.02520.01340.01240.01180.02210.0331
0.04860.02760.02730.027590.00750.00700.01940.01790.01740.00980.03700.02270.01410.00850.03090.01930.03440.0126
0.06350.0308TABLE 2
Quantitative results using the synthetic dataset in comparison with Mitsunaga and Nayar’s method (MN) in terms of
the RMSE and disparity over various datasets.

RAW images with the ground truth RAW images. The
error in the estimated RAW images with calibration isdue to a gamut mapping [37] at saturated pixels anddemosaicing/compression around edges. For comparison, we also estimate the response function using the
MN method from the same input data. With reasonable
input, both methods present similar estimation results.
To evaluate the performance using diverse response
functions, we additionally perform an experiment witha Point Grey Flea3 camera. With the camera, we canchoose any gamma function as a response function.
Three different gamma functions (γ = {0.6, 1.0, 2.2}) areused, and the result is shown in Fig. 9. For various typesof gamma functions, our method performs well while
the MN method fails to estimate the right solution for
γ = 0.6.
Radiometric alignment. To perform radiometric alignment without any exposure information, all we need isa simple constraint to anchor the solution. In our case,
we derive a response function, which satisfies that theobservation 0.5 is mapped to an irradiance 0.5. Note thatthe constraint is arbitrary for the radiometric alignment
application as long as a unique radiometric response

function, which may differ from the true response function by the exponential ambiguity, is determined. Basedon this, we estimate the (pseudo) exposure ratios of input
images and perform the radiometric alignment.
Here we show a radiometric alignment result using
the estimated response function in Fig. 10. Without
calibration, the aligned image with proper intensity normalization using the ratio of (R+G+B) intensity gives an
RMSE of 0.0657. After calibration with our method, theerror is reduced to 0.0225.
4.2

Varying illumination caseWe apply our method to the varying illumination caseand compare with the method proposed by Shi et al. [4].
We use two real-world scenes S HEEP and F ROG recorded
by a Canon 20D and Nikon D70, respectively. The
ground truth is obtained by the MN method using a
Macbeth Color Chart beforehand.
For each scene, we randomly select samples to produce color profiles and compute the average of theresults of ten trials. In practice, the accuracy of Shi etal.’s method has dependency on the selected samples,
therefore we take more trials and remove obvious failure

11
0.2

Inverse response function - Red
1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2

0.05

0.1

(I) sRGB (Input)(III) RAW (ground truth)

0(V) Error map (without calibration)
8001000

0

0.1

0.2

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

0.9

1

Inverse response function - Green1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2

0.05

0.1

(II) sRGB (Calibrated)(IV) RAW (estimated)0(VI) Error map (with calibration)(a) Canon 550D

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

8001000
0.2

Inverse response function - Red
1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2

0.05

0.1

(I) sRGB (Input)(III) RAW (ground truth)

0(V) Error map (without calibration)
8001000

0

0.1

0.2

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

Inverse response function - Green1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2

0.05

0.1

(II) sRGB (Calibrated)(IV) RAW (estimated)0(VI) Error map (with calibration)(b) Nikon D1X

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

8001000
0.2

Inverse response function - Red
1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2

0.05

0.1

(I) sRGB (Input)(III) RAW (ground truth)

0(V) Error map (without calibration)
8001000

0

0.1

0.2

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

Inverse response function - Green1
Mitsunaga & Nayar
Ours

0.9
0.8

0.1

0.05

Normalized irradiance

0.15

0.7
0.6
0.5
0.4
0.3
0.2
0.1

(II) sRGB (Calibrated)(IV) RAW (estimated)0(VI) Error map (with calibration)(c) Sony DSC-F828

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.88001000

Fig. 8. Results of our calibration method. From top to bottom, the results of different cameras are shown: (a) Canon
500D, (b) Nikon D1X, and (c) Sony DSC-F828. For each camera, several images are shown: (I) one of the input
images, (II) the calibrated image by applying estimated inverse response functions, (III) the ground truth RAW image,(IV) the estimated RAW image from the calibrated sRGB image, (V) error map between the ground truth RAW image
and the estimated RAW image without calibration, (VI) error map between the ground truth RAW image and theestimated RAW image with calibration. The white pixels on the error maps indicate saturated pixels having a value of
255 in any of the channels. The right column shows estimated inverse response functions.

12

Inverse response function - Green

Inverse response function - Green

1

0.9
0.8

0.7

Normalized irradiance

Normalized irradiance

0.8

0.6
0.5
0.4
0.3

0.5
0.4
0.3

0.1

0.1
0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

0.8

0.6

0.2

0

1

Groundtruth (g=2.2)
Mitsunaga & Nayar
Ours

0.9

0.7

0.2

0

1
Groundtruth (g=1.0)
Mitsunaga & Nayar
Ours

Normalized irradiance0.9

0Inverse response function - Green

1
Groundtruth (g=0.6)
Mitsunaga & Nayar
Ours

0.7
0.6
0.5
0.4
0.3
0.2
0.1

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

0

1

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

1

Fig. 9. Results of our calibration method with a dataset captured from a Point Grey Flea3 camera.

Inverse response function

Inverse response function1

1
Ground truth
Shi et al.
Ours

0.9

0.6
0.5
0.4
0.3

modes. In this experiment, for the purpose of directlycomparing our rank constraint with Shi et al.’s nonlinearity constraint, we do not apply the outlier rejection
process.
Fig. 11 shows the results. In the figure, the median
response curve of ten trials is used for plotting. We
also summarize the quantitative results in Table 3. For
both datasets, our method shows more accurate andstable results. With a small number of samples, theresult of the original method fluctuates because of theheuristic non-linearity constraint, but our method robustly converges well because of the more principled
computation scheme. With more samples, both methods
give more stable and robust estimation. As shown in
the table, our method is about 1000 times faster thanthe original method that uses Nelder-Mead simplex
optimization [38]. Since our method has a computational
complexity of O(n), we can use more samples to achieve
far more robust estimation.
4.3

Color mixture caseWe also perform radiometric calibration from a single
image, which utilizes the linearity of color mixtures
around texture edges [3]. This method relies on priorinformation of response functions because color mixtures
usually only span a limited range of observations, which

0.7
0.6
0.5
0.4
0.3

0.2

0.2

0.1

0.1

0

Fig. 10. Radiometric alignment of real-world images. Top:
input images, Bottom: radiometrically aligned images.

0.8

0.7

Normalized irradiance

Normalized irradiance

0.8

0

0.1

0.2

Ground truth
Shi et al.
Ours

0.9

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

1

0

0

S HEEP (12)

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

1

F ROG (20)Fig. 11. (Top) example images of S HEEP and F ROG
dataset. (Bottom) result of radiometric calibration in thevarying illumination case. S HEEP dataset is recorded by a
Canon 20D (12 images) and F ROG data is recorded by a
Nikon D70 (20 images).
is not sufficient to uniquely determine the response
function. Therefore, Lin et al. [3] infer the solution using a
Bayesian approach with prior information of real-world
camera response functions [7]. In [3], the linearity ofcolor mixtures is measured using the distance from each
mixed color to the line defined by the two colors in the
RGB color space.
We also use prior information to regularize the solution as explained in Sec. 3.6. There is a trade-off inusing the prior information between estimation bias androbustness. Therefore, a strong prior might introduce
bias in the estimation result. In our implementation, we
use five kernels (K = 5) to construct the prior modeland λ is set empirically to 103 . The rank error function
E is defined asE(g; D) =m
X

wi κ2 (Ai )s.t.

Ai = g ◦ Di ,

(23)

i=1

where m is the number of color mixtures, and wi is a

13

DatasetsamplesmethodShi et al. [4]50Ours
Shi et al. [4]
100
Ours

S HEEP (Canon 20D)
RMSE (std.)
Disparity (std.)
0.0488 (0.0227)
0.0846 (0.0429)
0.0226 (0.0006)
0.0447 (0.0061)
0.0418 (0.0172)
0.0729 (0.0345)
0.0218 (0.0003)
0.0430 (0.0021)

Time(s)433.660.84
914.64
1.11

F ROG (Nikon D70)
RMSE (std.)
Disparity (std.)0.0279 (0.0050)
0.0456 (0.0064)
0.0191 (0.0010)
0.0419 (0.0051)
0.0293 (0.0038)
0.0441 (0.0056)
0.0176 (0.0005)
0.0396 (0.0032)

Time(s)697.050.80
1316.88
0.86

TABLE 3
Radiometric calibration results in the varying illumination case. S HEEP/F ROG dataset captured by Canon 20D and
Nikon D70 are used. The average RMSE of ten trials are shown. The value inside the brackets is the standard
deviation of the ten trials.

Inverse response function - Green

the proposed approach is verified using both simulation
and real-world experiments.

1
Mitsunaga & Nayar
Lin's method
Ours

0.9

Normalized irradiance

0.8
0.7
0.6

5.1

0.5
0.4
0.3
0.2
0.1
0

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized observation

0.8

0.9

1

Fig. 12. Radiometric calibration from a single image.
Lin et al.’s result is obtained from their original paper. Data
courtesy of Stephen Lin.

Limitations

In our current implementation, we use a simple optimizer for the rank minimization, and it is sufficient
to find the right solution as shown in experiments formost of the cases. However, sometimes the optimization
fails to converge to the global minimum, or converges
to a trivial solution due to insufficient input data or
great noise. This could be improved by using a better
optimization method for the rank minimization.

ACKNOWLEDGMENTS
weighting parameter. The parameter wi is defined as the
distance between the two irradiances of a color mixture
in the RGB space because a color mixture that has
longer distance is affected more by response functions.
The coefficients of g is initially set to the mean inverse
response function of the DoRF database [7].
Fig. 12 shows the result for the input image used in [3].
The estimation result is compared with the ground truththat is obtained by MN method beforehand. While thiscase requires prior information, the result shows that
low-rank constraint can be applied well in color mixture
based methods.

The first author is partially supported by National Strategic R&D Program for Industrial Technology, Korea, andby the National Research Foundation of Korea(NRF)
grant funded by the Korea government(MEST) (No.20110018250)

R EFERENCES
[1]
[2]
[3]

5

C ONCLUSIONSIn this paper, we have introduced a robust radiometriccalibration algorithm that uses the low-rank structure ofirradiance vectors. The problem is formulated as oneof rank minimization and solved by minimization of
the condition number of the input matrix. We further
developed an algorithm for outlier rejection. These methods are combined in an iterative manner for achievingrobust estimation. The advantage of our method is thatit capitalizes on the fundamental linear dependency of
sensor irradiances, therefore it can handle the radiometric calibration problem in a unified and principledframework. We show that our framework can be easily
applied to other radiometric calibration problems thatuse linearity of irradiance by taking Shi et al.’s [4] andLin et al.’s [3] methods as examples. The effectiveness of

[4]

[5]
[6]

[7]

[8]

P. Debevec and J. Malik, “Recovering high dynamic range radiance maps from photographs,” in Proceedings of ACM SIGGRAPH,
1997, pp. 369–378. 1, 2
T. Mitsunaga and S. K. Nayar, “Radiometric self calibration,”in Proceedings of IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 1999, pp. 374–380. 1, 2, 4, 7, 8
S. Lin, J. Gu, S. Yamazaki, and H.-Y. Shum, “Radiometric calibration from a single image,” in Proceedings of IEEE Conference on
Computer Vision and Pattern Recognition (CVPR), 2004, pp. 938–945.
1, 3, 4, 6, 7, 12, 13
B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan, “Self-calibrating
photometric stereo,” in Proceedings of IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2010, pp. 1118–1125. 1, 3,
4, 6, 10, 13
S. J. Kim and M. Pollefeys, “Robust radiometric calibration andvignetting correction,” IEEE Transactions on Pattern Analysis andMachine Intelligence, vol. 30, no. 4, pp. 562–576, 2008. 1, 2, 7
S. Mann and R. Picard, “On being ‘undigital’ with digital cameras: extending dynamic range by combining differently exposed
pictures,” in Proceedings of IS & T, 48th annual conference, 1995, pp.
422–428. 2
M. Grossberg and S. K. Nayar, “Modeling the space of camera
response functions,” IEEE Transactions on Pattern Analysis andMachine Intelligence, vol. 26, no. 10, pp. 1272–1282, 2004. 2, 7,
12, 13S. K. Nayar and T. Mitsunaga, “High dynamic range imaging:spatially varying pixel exposures,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2000, pp.
472–479. 2

14

[9]

S. Mann and R. Mann, “Quantigraphic imaging: Estimating thecamera response and exposures from differently exposed images,”in Proceedings of IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2001, pp. 842–849. 2[10] M. Grossberg and S. K. Nayar, “Determining the camera response
from images: what is knowable?” IEEE Transactions on PatternAnalysis and Machine Intelligence, vol. 25, no. 11, pp. 1455–1467,
2003. 2, 6
[11] A. Litvinov and Y. Y. Schechner, “Addressing radiometric nonidealities: a unified framework,” in Proceedings of IEEE Conference on
Computer Vision and Pattern Recognition (CVPR), 2005, pp. 52–59.
2[12] H. Farid, “Blind inverse gamma correction,” IEEE Transactions on
Image Processing, vol. 10, no. 10, pp. 1428–1433, 2001. 2[13] Y. Tsin, V. Ramesh, and T. Kanade, “Statistical calibration of
CCD imaging process,” in Proceedings of International Conference
on Computer Vision (ICCV), 2001, pp. 480–487. 3[14] C. Pal, R. Szeliski, M. Uyttendale, and N. Jojic, “Probability
models for high dynamic range imaging,” in Proceedings of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2004,
pp. 173–180. 3[15] T.-T. Ng, S.-F. Chang, and M.-P. Tsui, “Using geometry invariants
for camera response function estimation,” in Proceedings of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2007,
pp. 1–8. 3[16] Y. Matsushita and S. Lin, “Radiometric calibration from noise
distributions,” in Proceedings of IEEE Conference on Computer Vision
and Pattern Recognition (CVPR), 2007, pp. 1–8. 3, 7
[17] J. Takamatsu, Y. Matsushita, and K. Ikeuchi, “Estimating camera
response functions using probabilistic intensity similarity,” inProceedings of IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2008, pp. 1–8. 3, 7
[18] ——, “Estimating radiometric response functions from image
noise variance,” in Proceedings of European Conference on Computer
Vision (ECCV), 2008, pp. 623–637. 3, 7
[19] S. Lin and L. Zhang, “Determining the radiometric response
function from a single grayscale image,” in Proceedings of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2005,
pp. 66–73. 3, 7
[20] B. Wilburn, H. Xu, and Y. Matsushita, “Radiometric calibration
using temporal irradiance mixtures,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2008,
pp. 1–7. 3, 4[21] C. Manders, C. Aimone, and S. Mann, “Camera response function
recovery from different illuminations of identical subject matter,”in Proceedings of International Conference on Image Processing (ICIP),
2004, pp. 2965–2968. 3[22] K. Shafique and M. Shah, “Estimation of the radiometric response
functions of a color camera from differently illuminated images,”in Proceedings of International Conference on Image Processing (ICIP),
2004, pp. 2339–2342. 3[23] S. J. Kim, J.-M. Frahm, and M. Pollefeys, “Radiometric calibration
with illumination change for outdoor scene analysis,” in Proceedings of IEEE Conference on Computer Vision and Pattern Recognition(CVPR), 2008, pp. 1–8. 3, 7
[24] M. Turk and A. Pentland, “Eigenfaces for recognition,” Journal of
Cognitive Neuroscience, vol. 3, no. 1, pp. 71–86, 1991. 3[25] C. Tomasi and T. Kanade, “Shape and motion from image streams
under orthography: a factorization method,” International Journal
on Computer Vision, vol. 9, no. 2, pp. 137–154, 1992. 3[26] 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. 3[27] Y. Peng, A. Ganesh, J. Wright, W. Xu, and Y. Ma, “RASL: Robustalignment by sparse and low-rank decomposition for linearlycorrelated images,” in Proceedings of IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2010, pp. 763–770. 3, 5[28] M. Paladini, A. Bartoli, and L. Agapito, “Sequential nonrigid structure-from-motion with the 3D-implicit low-rank shapemodel,” in Proceedings of European Conference on Computer Vision(ECCV), 2010, pp. 15–28. 3[29] Z. Zhang, A. Ganesh, X. Liang, and Y. Ma, “Tilt: Transforminvariant low-rank textures,” in Proceedings of Asian Conference on
Computer Vision (ACCV), 2010, pp. 314–328. 3, 5[30] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma,“Robust photometric stereo via low-rank matrix completion andrecovery,” in Proceedings of Asian Conference on Computer Vision(ACCV), 2010, pp. 703–717. 3, 5[31] Z. Zhang, Y. Matsushita, and Y. Ma, “Camera calibration withlens distortion from low-rank textures,” in Proceedings of IEEE
Conference on Computer Vision and Pattern Recognition (CVPR), 2011,
pp. 2321–2328. 3[32] J.-Y. Lee, B. Shi, Y. Matsushita, I.-S. Kweon, and K. Ikeuchi, “Radiometric calibration by transform invariant low-rank structure,”in Proceedings of IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2011, pp. 2337–2344. 3[33] J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma, “Robust
principal component analysis: Exact recovery of corrupted lowrank matrices by convex optimization,” in Proceedings of Neural
Information Processing Systems, 2009. 4, 5, 6[34] 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, 2011. 5[35] A. Chakrabarti, D. Scharstein, and T. Zickler, “An empiricalcamera model for internet color vision,” in Proceedings of British
Machine Vision Conference (BMVC), 2009. 7, 8[36] M. Grossberg and S. K. Nayar, “What is the space of camera
response functions?” in Proceedings of IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2003, pp. 602–609. 8[37] H. T. Lin, S. J. Kim, S. Susstrunk, and M. S. Brown, “Revisitingradiometric calibration for color computer vision,” Proceedings ofInternational Conference on Computer Vision (ICCV), 2011. 9, 10[38] J. Nelder and R. Mead, “A simplex method for function minimization,” Computer Journal, vol. 7, pp. 308–313, 1965. 12

Joon-Young Lee received the B.S degree inElectronic Engineering from Yonsei University
in 2008, and the M.S degree in Electrical Engineering from KAIST in 2009. He is currently
working toward the Ph.D. degree in ElectricalEngineering at KAIST. His research interests
include photometric methods in computer visionand 3D modeling. He is a student member of theIEEE.

Yasuyuki Matsushita received his B.S., M.S.
and Ph.D. degrees in EECS from the Universityof Tokyo in 1998, 2000, and 2003, respectively.
He joined Microsoft Research Asia in April 2003.
He is a Lead Researcher in Visual Computing
Group of MSRA. His major areas of research arecomputer vision (photometric techniques, suchas radiometric calibration, photometric stereo,
shape-from-shading), computer graphics (image
relighting, video analysis and synthesis). Dr.
Matsushita is on the editorial board memberof IEEE Transactions on Pattern Analysis and Machine Intelligence(TPAMI), International Journal of Computer Vision (IJCV), IPSJ Journalof Computer Vision and Applications (CVA), The Visual Computer
Journal, and Encyclopedia of Computer Vision. He served/is serving asa Program Co-Chair of PSIVT 2010, 3DIMPVT 2011 and ACCV 2012.
He is appointed as a Guest Associate Professor at Osaka University(April 2010-) and National Institute of Informatics, Japan (April 2011-).
He is a senior member of IEEE.

15

Boxin Shi received the B. Eng. degree from Beijing University of Posts and Telecommunicationsin 2003 and M. Eng. degree from Peking University in 2007. He joined the graduate schoolof information science and technology at the
University of Tokyo in Oct. 2010, where currently
he is a Ph.D. candidate under the supervisionof Prof. Katsushi Ikeuchi. His research interests
include photometric methods in computer vision,reflectance and illumination modeling.

In So Kweon received the BS and MS degrees
in mechanical design and production engineering from Seoul National University, Seoul, Korea, in 1981 and 1983, respectively, and the
PhD degree in robotics from the Robotics Institute, Carnegie Mellon University, Pittsburgh,
PA, 1990. He worked for Toshiba R&D Center,
Japan, and joined the Department of Automation
and Design Engineering, KAIST, Seoul, Korea,
in 1992, where he is now a Professor with the
Department of Electrical Engineering. He is arecipient of the best student paper runner-up award at the IEEE Conference on Computer Vision and Pattern Recognition (CVPR 09). His
research interests are in camera and 3-D sensor fusion, color modelingand analysis, visual tracking, and visual SLAM. He was the program
co-chair for the Asian Conference on Computer Vision (ACCV 07) andis the general chair for the ACCV 12. He is also on the editorial boardof the International Journal of Computer Vision. He is a member of theIEEE and the KROS.

Katsushi Ikeuchi received the BE degree from
Kyoto University in 1973 and the PhD degreefrom the University of Tokyo in 1978. After working at the Massachusetts Institute of Technology(MIT) Artificial Intelligence (AI) Laboratory forthree years, ETL for five years, and Carnegie
Mellon University (CMU) Robotics Institute for 10
years, he joined the University of Tokyo in 1996
and is currently a full professor. His researchinterest spans computer vision, robotics, andcomputer graphics. In these research fields, he
has received several awards, including the David Marr Prize in computational vision for the paper “Shape from Interreflection,” and IEEE
Robotics and Automation Society K. S. Fu memorial best transaction
paper award for the paper “Toward Automatic Robot Instruction from
PerceptionMapping Human Grasps to Manipulator Grasps.” In addition,in 1992, his paper “Numerical Shape from Shading and Occluding
Boundaries” was selected as one of the most influential papers to haveappeared in the Artificial Intelligence Journal within the past 10 years.
His IEEE activities include general chair of IROS’95, ITSC’00, IV’01;
program chair of the CVPR’96 and ICCV’03; associate editor of IEEE
Transactions on Robotics and Automation, IEEE Transactions on PatternAnalysis and Machine Intelligence; and a distinguished lecturer of the
Signal Processing Society in 2000-2002 and Robotics and Automation
Society in 2004-2006. He was elected an IEEE fellow in 1998. He is
the editor-in-chief of the International Journal of Computer Vision. He
received the Computer Vision Significant Researcher Award in ICCV2011.