Radiometric Calibration by Transform Invariant Low-rank Structure
Joon-Young Lee1∗ Boxin Shi2 Yasuyuki Matsushita3 In-So Kweon1 Katsushi Ikeuchi2
1
Robotics and Computer Vision Lab, KAIST, Korea
2
The University of Tokyo, Japan
3
Microsoft Research Asia, China
jylee@rcv.kaist.ac.kr, shi@cvl.iis.u-tokyo.ac.jp, yasumat@microsoft.com,
iskweon@ee.kaist.ac.kr, ki@cvl.iis.u-tokyo.ac.jp

Abstract

1. Introduction
In most cameras, there exists a radiometric response
function that relates sensor irradiance and output brightness
values. The function is typically designed to be non-linear
for purposes such as compressing the dynamic range of sensor irradiance or for adapting to the display’s non-linear
mapping. While many computer vision algorithms assume
linear (or affine) relationship between the sensor irradianceand the recorded intensity, the radiometric response functions are generally non-linear and unknown. Moreover, theshape of the functions may even vary with camera parameter settings. Therefore, radiometric calibration is an important first step for various computer vision algorithms thatassume the linear relationship between the irradiance andobservation to make them work correctly.
When multiple images are recorded from a static scene
with different exposure times, the sensor irradiances at corresponding pixel locations are linearly related by the ratioof exposure times. If the response function is linear, the∗ This work was done while Joon-Young Lee was an intern at Microsoft
Research Asia.

Linear-dependent

Observation

We present a robust radiometric calibration method that
capitalizes on the transform invariant low-rank structure of
sensor irradiances recorded from a static scene with different exposure times. We formulate the radiometric calibration problem as a rank minimization problem. Unlike previous approaches, our method naturally avoids over-fitting
problem; therefore, it is robust against biased distribution of
the input data, which is common in practice. When the exposure times are completely unknown, the proposed method
can robustly estimate the response function up to an exponential ambiguity. The method is evaluated using both simulation and real-world datasets and shows a superior performance than previous approaches.

Our method
Response function

Linear-independent
Irradiance

Irradiance

Imaging process

Observation

Figure 1. Illustration of our approach. Our method takes a set ofimages taken from a fixed view point with varying exposure times.
In the matrix form, by putting images as column vectors, the irradiance matrix should be rank 1, while the observations matrix has
the higher rank. Our method seeks the inverse response functionthat transforms the observation matrix into a low-rank matrix.

recorded pixel values (observations) show the linear relationship across the images. However, when the response
function is non-linear, the linear relationship breaks down.
In this paper, we present a new radiometric calibration
method that capitalizes on the fundamental linear dependency of sensor irradiances. Our method arranges the observations obtained from multiple exposure times in a matrix
form, where a column corresponds to an image vector, anda row corresponds to a pixel location. When the linear relationship is observed in such a matrix, the rank of the matrix
becomes one. On the other hand, when the response function is non-linear, the matrix generally becomes full rank.
From this observation, we formulate the estimation of radiometric response functions as a rank minimization problem.
Contributions We show that the radiometric calibration
problem can be robustly and efficiently solved by a rank
minimization approach. Unlike previous approaches, we
explicitly use the linear dependency of the irradiance vectors. Our rank minimization approach naturally avoids theover-fitting problem that is common in previous approaches
that rely on the `2 -norm minimization. In addition, the pro-posed method can recover the response function up to an
exponential ambiguity even when no information is available about the exposure times. Such recovered functions
are useful for applications such as high-dynamic range image production and image stitching. In addition, we showthat the same formulation can be applied to previous radiometric calibration approaches that use the linearity of theirradiance. As an example, we show the radiometric calibration from photometric stereo input images where thelighting directions vary across the images [15].
In presenting our work, we begin by briefly describing the radiometric response function and reviewing related
works.

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)

With assumptions of monotonicity and continuity, the response function f can be uniquely inverted to an inverse
response function g = f −1 . With the inverse response
function g, the recorded intensity values can be linearizedby g(M ). Since what is available as input are the observations M , the radiometric calibration aims at recovering theinverse response function g.

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 of Mann and Picard [10] used a gamma
correcting function to represent response functions. With
known exposure ratios, their method can successfully recover the inverse response function in the parametric form.
Mitsunaga and Nayar used a different parametric representation, i.e., a polynomial representation. Using an approximate knowledge of relative exposure times, their approach
estimates response functions. Grossberg and Nayar [4] created a database of response functions (DoRF) and used it for
better representation of response functions. Debevec andMalik [1] assumed a smoothness property of the response
functions and estimated them in a non-parametric manner.
There are a few prior approaches that allow some camera
movement or scene motion. In Mann’s [9] method, images
are taken by a rotating and zooming camera. Kim and Pollefeys’s method [5] allows for free movement of the cameraand some motion in the scene.
As pointed out in [3, 5, 8], without the knowledge of
exposure ratios, the estimate still has an exponential ambiguity. While not unique, such an estimate is still usefulfor many applications, such as radiometric alignment, HDR
image production, and image stitching.

Instead of using varying exposure times, some approaches use statistical properties embedded in images toachieve radiometric calibration. Tsin et al.’s method [17]estimates non-parametric response functions using a statistical model of the CCD imaging process. Pal et al. [14]used probabilistic imaging models and weak prior models
for deriving response functions to produce high-quality highdynamic range images. Matsushita and Lin [11] proposed
to use the symmetric property of image noise by observing noise distributions contained in images. Takamatsu etal. [16] improved the noise-based method with a probabilistic intensity similarity measure, which requires fewer
number of images. Lin et al. [6] and Lin and Zhang [7]proposed methods that take only a single image as input.
Their methods use edges for obtaining color or gray-scale
histogram distributions, and the optimal inverse response
function is determined by transforming obtained non-linear
distributions into linear distributions. A similar idea is recently introduced by Shi et al. [15] that achieves calibrationfrom images taken under varying lighting conditions.

2. Low-rank Structure of Calibrated Observations
Our method casts the radiometric calibration problem asa low-rank structure recovery problem. Suppose we are
given multiple images taken from a static scene with various
exposures. The sensor irradiance I and observed intensity
value M are related by a radiometric response function f by
Eq. (1). Because the images are taken from a static scene,
the scene radiance L is constant, and the sensor irradiancebecomes proportional to the exposure time e with a constantscaling k:
I = kLe ∝ e.
(2)
Consider an observation matrix D ∈ Rm×n , where mand n are the numbers of pixels and images, respectively,
created by stacking all the vectorized images vec(Mi ) as
D = [vec(M1 ) | · · · | vec(Mn )] ,(3)

T

where vec(Mi ) = [Mi (1), · · · , Mi (m)] for i =(1, . . . , n). If we know an inverse radiometric response
function g = f −1 , the observation matrix D can be transformed into an irradiance matrix A with a scaling ambiguitybyg ◦ D = A = [vec(I1 ) | · · · | vec(In )] ,(4)
where ◦ is an operator that describes the element-wise mapping. Since the irradiance vectors vec(Ii ) should be lineardependent, the rank of the irradiance matrix A becomes one.
Therefore, our method seeks out the inverse response function g that minimizes the rank of the transformed matrixg ◦ D(= A) asĝ = argmin rank (A)gs.t.

A = g ◦ D.

(5)

1

ĝ = argmingnXκi (A)s.t.

A = g ◦ D.

0.6
0.6
0.4

0.4

0.2

0.8
0.6
0.4

0

0

0.2
0.4
0.6
0.8
Normalized measurement

0

1

0.6
0.4
0.2

0.2

0.2
0

0.8

0

0

1

1

0.8

0.8

Normalized irradiance1

Normalized irradianceTo efficiently compute the rank minimization problem ofEq. (5), previous approaches use an approximate solution
method using
a nuclear norm (sum of the singular values;. Pn||A||∗ = i=1 σi (A)) minimization [18].
In our case, however, the nuclear norm minimization
cannot be directly employed because the function g alters
the absolute values of matrix elements, and it therefore results in the variations of the singular value magnitudes. To
resolve this issue, we instead use condition numbers, i.e., aratio of singular values. We define the condition numbers.
as κi (A) = σi (A)/σ1 (A), i = 2, 3, . . . , n. With the condition numbers, we can work with the relative magnitudesof the singular values that are unaffected by their absolutemagnitudes. Now we approximate the energy functional (5)as a minimization of the summation of condition numbersas:Normalized
irradiance
Normalized irradiance

1

0.8
0.8

Normalized irradiance

3. Calibration Algorithm

Normalized irradiance

1

0.6
0.4
0.2

0
0
0.4
0.8
0.2
0.40 0.6 0.20.8
1
0 0.60.2
0.4
0.61 0.8
Normalized measurement
Normalized measurement
Normalized measurement

0.2
0.4
0.6
0.8
Normalized measurement

1

0.6
0.4
0.2

0

0

0.2
0.4
0.6
0.8
Normalized measurement

1

1

(a)
(b)
(c)
(d)
Figure 2. Reference response functions used for Table. 1. (a) Ideal
linear response function (RF 1). (b,c,d) General non-linear response functions (RF 2-4).

κ20.00000.02890.11630.00290.00790.02930.11200.0094condition number
RF 1
without
RF 2
noise
RF 3
RF 4
RF 1
with noiseRF 2(σ = 0.005) RF 3
RF 4

(6)

κ3
0.00000.00340.01770.00000.00770.00590.01750.0087κ40.00000.00100.00250.00000.00750.00500.00540.0085κ50.00000.00050.00010.00000.00720.00480.00460.0081

Table 1. Condition number variations

i=2

There are two main factors causing rank variations: nonlinearity of response function and image noise. The top halfof Table. 1 shows condition numbers from synthetic datathat have five different exposure times. We assume the zeromean Gaussian characteristics for image noise. Response
functions in Fig. 2 are used as reference and we averaged
the results of 100 trials for noisy cases. Fig. 2 (a) is an ideal
linear response function. With this linear response function,
all condition numbers are almost zero which means rank1 structure. With other more general non-linear response
functions, however, the condition number κ2 becomes comparably larger than other condition numbers. This is commonly observed because of the monotonic and smooth characteristics of response functions. On the other hand, theeffect of noise also changes the condition numbers. As
shown in the bottom half of Table. 1, image noise makes
all the condition numbers pretty even because of its highfrequency nature. From the above observations, our method
only uses the condition number κ2 instead of using all thecondition numbers. Now the objective function is changed
from Eq. (6) toĝ = argmin κ2 (A)s.t.

A = g ◦ D.

(7)g

To achieve a better convergence, we add more constraints tothe objective functions as described in the next subsection.

3.1. Additional constraints
In general, response functions as well as inverse response
functions are monotonic and smooth. To enforce the monotonicity of inverse response functions g, we put in a monotonicity constraint, which is represented as ∂g/∂M > 0.

In addition, we use the smoothness constraint for the inverse response function g as done in [1] to avoid the effect of noise. This smoothness constraint is represented
as min ∂g 2 /∂ 2 M . The smoothness constraint, however,
could cause under-fitting results in the case of complex response functions. Since the summation of the condition
numbers is highly correlated to noise level as discussed
above, we use the condition numbers for balancing theweight of the smoothness constraint. With these constraints,
the final form of the objective function to estimate response
function ĝ is defined asX  ∂g(t) ĝ = argmin κ2 (A) + λ1
H −
∂D
g
t
+ λ2

n
X
i=2

κi (A)

X ∂ 2 g(t)
t

∂D2

s.t.

A = g ◦ D, (8)where H(·) is the Heaviside step function (H(x) = 1
when x ≥ 0, and H(x) = 0, otherwise). The derivatives∂g(t)/∂D and ∂ 2 g(t)/∂D2 are assessed by sampling g(t)at various t in the range of [0, 1]. The boundary conditions
defined as g(0) = 0 and g(1) = 1 are embedded in parametric representation of g. More specifically, we represent the
response function by a polynomial representation as donein [12] with the explicit boundary conditions:g(t) = t + t(t − 1)n−1
X

ci tn−i−1 ,

(9)

i=1

where n is the order of a polynomial function, and ci are thecoefficients to estimate. To obtain a greater representation
power, our method can alternatively use other representations, such as EMoR [4].

3.2. Exponential ambiguity

γi,jwhere ri,j is the exposure ratio ei /ej of measurement pairs
Mi and Mj .

4. Experiments
To evaluate the proposed method, we perform experiments using both simulation and real-world data. We
compare our results with results of Mitsunaga and Nayar’s
method [12] (MN method, here-forth) as both methods takethe same input. In our implementation, we used a 6th order polynomial function to represent response functions as
done in [12].
The optimization of Eq. (8) is performed using the
Nelder-Mead simplex method [13]. We choose this method
due to its simplicity, which can be implemented easily usinga Matlab built-in function “fminsearch”. In all the experiments, we use the linear function as the initial guess, and
it is confirmed that the proposed method has a good convergence property even with such simple initialization. We
always set λ1 = 1 and λ2 = 1 in Eq. (8) during the experiment. Once the response function is estimated up to the
exponential ambiguity, we estimate the final response function using Eq. (10). This optimization is also performed
using the Nelder-Mead simplex method [13].

4.1. Simulation
In this section, we use a synthetic dataset for quantitative evaluation. The synthetic dataset is generated using the
DoRF database [4], which contains 201 measured response

D2

D3

D4

Likelihood

As discussed by Grossberg and Nayar in [2], recoveryof the exposure ratios and the response function is impossible without making assumptions on the response functionor having rough estimates on the exposure ratios. When
exposure times are completely unknown, there remains an
exponential ambiguity in the estimate. This exponential ambiguity means that if g is a solution of I = g(M ), then g γ
can also become a solution as I γ = g γ (M ) for any constant value γ. Our method also suffers from this exponential ambiguity like other methods when the exposure times
are completely unknown, because the rank-1 structure is retained even after any exponential function is applied. In
short, if rank(I) = 1, then rank(I γ ) = 1 for any γ.
However, 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. To resolve the ambiguity, our method only requires
at least one exposure ratio. Using the known exposure ratio, we can estimate γ by solving the optimization problem
described asX2γ̂ = argmin[ĝ γ (Mi ) − ri,j ĝ γ (Mj )] ,
(10)

D1

0

0.2

0.4

0.6

0.8

1

Normalized radiance

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

functions. We synthetically generate scene radiances in therange of [0, 1] and create 1000 observations (which forms
a 1000-pixel image) using each response function with fiveexposure times. The step of exposure ratios is set to 0.5
because this is a commonly available setting in commercial
cameras. Fig. 3 shows scene radiance distribution of thesynthetic dataset. To perform an evaluation, we add five different magnitudes of noise for each radiance distributions;
therefore, a comparison with the MN method is performed
using 1005 synthetic scenes in total. Gaussian noise withthe standard deviation σ(= 0, 0.0025, 0.0050, 0.0075, 0.01)
is used for diversifying the datasets.
Results with synthetic dataset Fig. 4 shows results withthe synthetic dataset. In this figure, we plot cumulative histograms that show the number of successful cases, i.e., successful estimation of response functions under a certain root
mean squared error (RMSE) with each configuration. As
shown in the figure, the estimation performance is dependent not only on the noise level but also on the radiance distribution. With uniformly distributed scene radiance, the estimation becomes quite stable. However, when the radiancedistribution is biased, which is very common, it introduces
performance degradation as shown in the D4 result.
We summarize quantitative results in Table. 2. The
DoRF database contains some response functions which are
difficult to represent with 6th order polynomials. These
functions give trivial results for both our method and the
MN method; therefore, we only use the best 150 fitting results for each of the algorithms to compute the mean andstandard deviation of RMSE and disparity (the maximum
deviation from the ground truth). The tables shows our
method performs well for biased datasets (D2-D4), while
the MN method works better with the uniform data (D1).
Our method is robust against the increasing noise level,
while the performance of the MN method rapidly degrades.

4.2. Real-world experiment
We also perform experiments using real-world cameras
and scenes. In our experiments, we used four different cameras: Canon 20D, Nikon D70, Sony A200 and Olympus E1. Each dataset is collected by capturing a static scene withdifferent exposure times. Fig. 7 shows estimation results

160

140

140

120

120

120

120

100
80

MN,  =0
MN,  =0.0025

60

MN,  =0.005
MN,  =0.0075
Ours,  =0

40

100
80

MN,  =0
MN,  =0.0025

60

MN,  =0.005
MN,  =0.0075
Ours,  =040

Ours,  =0.0025
20

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

80

MN,  =0
MN,  =0.0025

60

MN,  =0.005
MN,  =0.0075
Ours,  =0

40

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

80

MN,  =0
MN,  =0.0025

60

MN,  =0.005
MN,  =0.0075
Ours,  =040

Ours,  =0.0025
20

Ours,  =0.005

Ours,  =0.0075
0100

Ours,  =0.0025
20

Ours,  =0.005

Ours,  =0.0075
0100

Ours,  =0.0025
20

Ours,  =0.005# of scenes < RMSE

160

140

# of scenes < RMSE

160

140

# of scenes < RMSE

# of scenes < RMSE

160

Ours,  =0.005

Ours,  =0.0075
0

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Ours,  =0.0075
0

0

0.005

0.01

0.015
RMSE

0.02

0.025

0.03

Figure 4. Cumulative histogram of the number of successful cases w.r.t. RMSE. From left to right, D1-D4 datasets are shown.

Sdv. of noise (σ)
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

0.00000.00420.00310.01130.00820.00690.00600.01630.01380.00650.00570.01540.01290.00830.00760.01970.01650.00250.00400.00560.01110.01440.00700.00790.01730.01780.00620.00620.01520.01410.01150.01070.02630.0219Mean
0.00500.00440.00810.01180.01960.01180.00880.02660.01980.00850.00720.01980.01640.02200.01350.04960.02840.00750.00560.00930.01390.02200.01770.01050.04020.02400.01340.00840.03030.01910.03500.01780.07670.03740.010
0.00830.01190.01880.02700.02440.01270.05540.02910.01860.00940.04230.02130.05280.02110.11270.04390.00000.00420.00230.01130.00650.00670.00480.01610.01110.00620.00430.01500.01000.00880.0071
0.02140.0150Standard Deviation
0.0025 0.0050 0.0075
0.0041 0.0042 0.0039
0.0041 0.0063 0.0070
0.0108 0.0105 0.0110
0.0103 0.0140 0.0161
0.0076 0.0104 0.0139
0.0060 0.0064 0.0082
0.0182 0.0230 0.0317
0.0138 0.0145 0.0188
0.0065 0.0076 0.0114
0.0047 0.0050 0.0055
0.0154 0.0168 0.0249
0.0110 0.0122 0.0137
0.0107 0.0166 0.0249
0.0082 0.0096 0.0116
0.0237 0.0362 0.0518
0.0164 0.0203 0.0247

0.01000.00450.00830.01220.01860.01770.00980.03960.02290.01320.00610.03020.01470.03580.01260.07090.0284

Table 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.

# of samples
MN
RMSE
Ours
MN
Disparity
Ours

100
0.02610.0141
0.05480.02951000
0.00650.00420.01390.008910000
0.00220.00130.00470.0029

Table 3. Result of stability test using the blue channel of the scene
Fig. 7 (a).

with real-world scenes. In the figure, the left column shows
one of the input images, and the middle and right columns
show the estimated response functions. For comparison, we
also estimate the response function using MN method fromthe same data input. RMSE and disparity on the figures arecalculated by comparing our result and the MN method’s
result by taking one of them as the ground truth. With reasonable data input, both methods present similar estimation
results. For the evaluation of the stability with respect tothe varying input data, we estimate the response function ofFig. 7 (a) using randomly sampled pixels (100, 1000, and10000 pixels). For each sampled pixel number, we test 100
trials and calculate the RMSE and disparity of both methods. For this test, we use the median of the estimation results with 10000 samples as the ground truth because theexperiment intends to test stability of the algorithms withdifferent input data. The result using the blue channel issummarized in Table. 3. As shown in the table, our method
produces stable estimates even when the number of sampled
pixels is small.
Radiometric alignment One of the advantages of ourcalibration algorithm is that our method does not require
any knowledge of exposure times when recovering the response function up to exponential ambiguity. Such estimates are useful for applications like radiometric alignment
and image stitching. Here we show a radiometric alignment result using the estimated response function in Fig. 5.
Without calibration, the aligned image with proper intensity
normalization using the ratio of (R+G+B) intensity gives an
RMSE of 0.2019. After calibration with our method, theerror is reduced to 0.0385.

Inverse response function

Inverse response function1

1

Ground truth
Shi et al.
Our Method

0.9

5. Application to other radiometric calibration
methods
The low-rank structure of the image irradiance matrix
is a powerful cue for recovering response functions. Our
method can be easily extended to deal with other radiometric calibration settings. In this section, we show an exampleof the extension by reformulating Shi et al. [15] for robustly
deriving the solution using the rank minimization scheme.
Shi et al. [15] showed that the radiometric calibration forphotometric stereo images can be performed by linearizing
the color profiles. The color profile is defined as an ordered
set of RGB color values in the RGB color space. Theirmethod capitalizes on the property that color profiles form
straight lines in the RGB color space if the response function is linear, while non-linear response functions f bendcolor profiles to non-linear shapes. Their method therefore
seeks the inverse response function g that linearizes all thecolor profiles.
It is straightforward to cast their problem into our lowrank computation scheme. If the color profiles are stacked
as row vectors to form a matrix, the linearization of color
profiles becomes equivalent to the matrix rank minimization
problem. More specifically, the color profile matrix should
become rank-1 when the correct inverse response function
is applied.
Following the notations of previous sections, we have an
image stack matrix D ∈ Rm×n , and we still use m andn as the numbers of pixels and images. One difference isthat the images are captured under n different lighting conditions. For each pixel, we have an n-dimensional row vector vec(Di ), i = 1, . . . , m, and this time we extend eachvec(Di ) into a 3×n vector to account for RGB color values.
0
We denote this observation matrix as Di . If a correct inverse0
response function g is given, it transforms the matrix Di toan irradiance matrix Ai , where rank(Ai ) = 1. Involving all
the pixels, we have the objective function written as:ĝ = argmin
g

m
X
i=1

rank (Ai )

s.t.

0

Ai = g ◦ Di .

(11)

0.8

0.7
0.6
0.5
0.4
0.3

0.7
0.6
0.5
0.4
0.3

0.2

0.2

0.1

0.1

0

Figure 5. Radiometric alignment of real-world images. Top: input images, Bottom: radiometrically aligned images. For display,
gamma correction(γ = 0.5) is applied.

0.9

Normalized irradiance

Normalized irradiance

0.8

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

Ground truth
Shi et al.
Our Method

0.8

0.9

1

0

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

1

S HEEP (12)F ROG (20)Figure 6. Images of the real-world scene with the number of images used and radiometric calibration results by linearizing thecolor profiles.

Thus, the problem can be solved in exactly the same way as
discussed in previous sections.
To verify the robustness and effectiveness of our method,
we performed experiments on 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 and the experiment
was performed on an Intel Core2 Duo E6750 (2.66GHz)CPU.
For each configuration, we randomly select samples to
produce color profiles and compute the average of the results of ten trials. In practice, the accuracy of the Shi etal.’s method has dependency on the selected samples, therefore we take more trials and remove obvious failure modes.
Fig. 6 shows the results. The median response curve often trials is used for plotting in the figure. We summarize the quantitative results in Table. 4. For both datasets,
our method shows more accurate and stable results. With asmall number of samples, the result of the original method
fluctuates because of the heuristic nonlinearity constraint,
but our method robustly converges well because of the 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 than the original
method. Since our method has a computational complexityof O(n), we can use more samples to achieve robust estimation.

6. Conclusions
In this paper, we have introduced a robust radiometriccalibration algorithm that uses the low-rank structure of irradiance vectors. The problem is formulated as a rank minimization and solved by minimization of the condition number of the input matrix. Our method can robustly estimate
response functions; when the exposure times are completely
unknown, it can estimate the response functions up to an exponential ambiguity. The advantage of our method is that it
can avoid over-fitting because our method does not rely on
the traditional least-square fitting. Additionally, we show

Dataset
Methodsamples50
Shi et al.
100
50
Our method
100

S HEEP (Canon 20D)
RMSE (std.)
Disparity (std.)
0.0488 (0.0227) 0.0846 (0.0429)
0.0418 (0.0172) 0.0729 (0.0345)
0.0219 (0.0002) 0.0429 (0.0027)
0.0218 (0.0002) 0.0423 (0.0022)

Time(s)284.01595.420.27
0.40

F ROG (Nikon D70)
RMSE (std.)
Disparity (std.)0.0279 (0.0050) 0.0456 (0.0064)
0.0293 (0.0038) 0.0441 (0.0056)
0.0183 (0.0007) 0.0393 (0.0027)
0.0176 (0.0004) 0.0390 (0.0023)

Time(s)478.67891.970.36
0.60

Table 4. Radiometric calibration results by linearizing the color profiles. S HEEP/F ROG dataset captured by Canon 20D and Nikon D70 are
used. Ten trial results are averaged. The value inside the brackets is the standard deviation of ten trials.

that our framework can be applied to other radiometric calibration problems that use linearization of observations bytaking Shi et al.’s method as an example. The effectivenessof the proposed approach is verified using both simulation
and real-world experiments.
Limitations Since our method is based on low-rank structure of irradiances, the result of our method always converge
to the direction of minimal rank. As we mentioned before,rank variations can be observed due to both nonlinear response function and image noise. Therefore, the rank could
be more affected by noise when the non-linearity of the response function is negligible compared with the noise level,
and in such a case, our method could produce inaccurate
estimation.
In current implementation, we use a simple optimizer forthe rank minimization, and it is usually sufficient to find theright solution as shown in experiments. However, sometimes the optimization fails to converge to the global solution, 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.

Acknowledgements
The first author is partially supported by MKE(The Ministry of Knowledge Economy), Korea, under the Human Resources Development Program for Convergence Robot Specialists support program supervised by the NIPA(National
IT Industry Promotion Agency) (NIPA-2010-C7000-10010007). The authors would like to thank Zhengdong Zhang
for helpful discussions.

References[1] P. Debevec and J. Malik. Recovering high dynamic range
radiance maps from photographs. In Proc. of ACM SIGGRAPH, pages 369–378, 1997. 2, 3
[2] M. Grossberg and S. Nayar. What can be known about theradiometric response from images? In Proc. of European
Conf. on Comp. Vis. (ECCV), pages 189–205, 2002. 4[3] M. Grossberg and S. Nayar. Determining the camera response from images: What is knowable?
IEEE Trans.
on Pattern Analysis and Machine Intelligence, 25(11):1455–
1467, 2003. 2[4] M. Grossberg and S. Nayar. What is the space of camera
response functions? In Proc. of Comp. Vis. and Pattern Rec.
(CVPR), pages 602–609, 2003. 2, 3, 4
[5] S. Kim and M. Pollefeys. Robust radiometric calibration andvignetting correction. IEEE Trans. on Pattern Analysis andMachine Intelligence, 30(4):562–576, 2008. 2[6] S. Lin, J. Gu, S. Yamazaki, and H. Shum. Radiometric calibration from a single image. In Proc. of Comp. Vis. andPattern Rec. (CVPR), pages 938–945, 2004. 2[7] S. Lin and L. Zhang. Determining the radiometric response
function from a single grayscale image. In Proc. of Comp.
Vis. and Pattern Rec. (CVPR), pages 66–73, 2005. 2[8] A. Litvinov and Y. Schechner. Addressing radiometric nonidealities: a unified framework. In Proc. of Comp. Vis. andPattern Rec. (CVPR), pages 52–59, 2005. 2[9] S. Mann. Comparametric imaging: Estimating both the unknown response and the unknown set of exposures in a plurality of differently exposed images. In Proc. of Comp. Vis.
and Pattern Rec. (CVPR), pages 842–849, 2001. 2[10] S. Mann and R. Picard. On being ‘undigital’ with digital
cameras: extending dynamic range by combining differently
exposed pictures. In Proc. of IS & T, 48th annual conference,
pages 422–428, 1995. 2[11] Y. Matsushita and S. Lin. Radiometric calibration fromnoise distributions. In Proc. of Comp. Vis. and Pattern Rec.
(CVPR), pages 1–8, 2007. 2[12] T. Mitsunaga and S. Nayar. Radiometric self-calibration. In
Proc. of Comp. Vis. and Pattern Rec. (CVPR), pages 374–
380, 1999. 3, 4[13] J. Nelder and R. Mead. A simplex method for function minimization. Computer Journal, 7:308–313, 1965. 4[14] C. Pal, R. Szeliski, M. Uyttendale, and N. Jojic. Probability
models for high dynamic range imaging. In Proc. of Comp.
Vis. and Pattern Rec. (CVPR), pages 173–180, 2004. 2[15] B. Shi, Y. Matsushita, Y. Wei, C. Xu, and P. Tan. Selfcalibrating photometric stereo. In Proc. of Comp. Vis. andPattern Rec. (CVPR), pages 1118–1125, 2010. 2, 6
[16] J. Takamatsu, Y. Matsushita, and K. Ikeuchi. Estimatingcamera response functions using probabilistic intensity similarity. In Proc. of Comp. Vis. and Pattern Rec. (CVPR), pages
1–8, 2008. 2[17] Y. Tsin, V. Ramesh, and T. Kanade. Statistical calibration of
ccd imaging process. In Proc. of Int’l Conf. on Comp. Vis.
(ICCV), pages 480–487, 2001. 2[18] J. Wright, A. Ganesh, S. Rao, Y. Peng, and Y. Ma. Robust
principal component analysis: Exact recovery of corrupted
low-rank matrices by convex optimization. In Proc. of Neural Information Processing Systems, 2009. 3

Inverse response function - Red

Inverse response function - Green1

1
Mitsunaga & Nayar
Our Method

0.9
0.8

0.8
Normalized irradiance

Normalized irradiance

0.9

RMSE = 0.0104
Disparity = 0.0185

0.7
0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3
0.2

0.1

0.1
0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

0

1

RMSE = 0.0076
Disparity = 0.0213

0.7

0.2

0

Mitsunaga & Nayar
Our Method

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

0.8

0.9

1

(a) Canon 20D
Inverse response function - Red

Inverse response function - Green1

1
Mitsunaga & Nayar
Our Method

0.9
0.8

0.8
Normalized irradiance

Normalized irradiance

0.9

RMSE = 0.0091
Disparity = 0.0243

0.7
0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3
0.2

0.1

0.1
0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

0

1

RMSE = 0.0220
Disparity = 0.0518

0.7

0.2

0

Mitsunaga & Nayar
Our Method

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement(b) Nikon D70
Inverse response function - Red

Inverse response function - Green1

1
Mitsunaga & Nayar
Our Method

0.9

0.9
0.8

RMSE = 0.0015
Disparity = 0.0032

0.7

Normalized irradiance

Normalized irradiance

0.8

0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3
0.2

0.1

0.1
0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

0

1

RMSE = 0.0091
Disparity = 0.0220

0.7

0.2

0

Mitsunaga & Nayar
Our Method

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement(c) Sony A200
Inverse response function - Green

Inverse response function - Blue1

1
Mitsunaga & Nayar
Our Method

0.9

0.9
0.8

RMSE = 0.0032
Disparity = 0.0069

0.7

Normalized irradiance

Normalized irradiance

0.8

0.6
0.5
0.4
0.3

0.6
0.5
0.4
0.3
0.2

0.1

0.1
0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement

0.8

0.9

1

RMSE = 0.0066
Disparity = 0.0178

0.7

0.2

0

Mitsunaga & Nayar
Our Method

0

0

0.1

0.2

0.3 0.4 0.5 0.6 0.7
Normalized measurement(d) Olympus E-1
Figure 7. Results of our calibration method. Left column: One of the input images. Middle/right column: Estimated inverse response
functions. The RMSE and disparity are calculated by comparing our estimation and the MN method’s estimation. From top to bottom, theresults of different cameras are shown: (a) Canon 20D, (b) Nikon D70, (c) Sony A200, and (d) Olympus E-1.