Estimating radiometric response functions
from image noise variance
Jun Takamatsu1? , Yasuyuki Matsushita2 , and Katsushi Ikeuchi1??
11

Microsoft Institute for Japanese Academic Research Collaboration (MS-IJARC)
j-taka@is.naist.jp, ki@cvl.iis.u-tokyo.ac.jp
2
Microsoft Research Asia
yasumat@microsoft.com
Abstract. We propose a method for estimating radiometric response
functions from observation of image noise variance, not profile of its
distribution. The relationship between radiance intensity and noise variance is affine, but due to the non-linearity of response functions, thisaffinity is not maintained in the observation domain. We use the nonaffinity relationship between the observed intensity and noise variance toestimate radiometric response functions. In addition, we theoretically derive how the response function alters the intensity-variance relationship.
Since our method uses noise variance as input, it is fundamentally robust
against noise. Unlike prior approaches, our method does not require images taken with different and known exposures. Real-world experiments
demonstrate the effectiveness of our method.

Introduction

Many computer vision algorithms rely on the assumption that image intensity islinearly related to scene radiance recorded at the camera sensor. However, thisassumption does not hold with most cameras; the linearity is not maintained in
the actual observation due to non-linear camera response functions. Linearizationof observed image intensity is important for many vision algorithms to work,
therefore the estimation of the response functions is needed.
Scene radiance intensity (input) I and observed image intensity (output) O
are related by the response function f as O = f (I). Assuming it is continuousand monotonic, the response function can be inverted to obtain the inverse
response function g (= f −1 ), and measured image intensities can be linearizedusing I = g(O). Since only observed output intensities O are usually available,most estimation methods attempt to estimate the inverse response functions.
1.1

Prior work

Radiometric calibration methods assume known characteristics of radiance ata camera to estimate unknown response functions. One class of these methods
uses information about the ratio of input radiance intensities. The Macbeth
color checker is used for estimating the response functions by several researchgroups such as [1]. Nayar and Mitsunaga [2] use an optical filter with spatially?

Nara Institute of Science and Technology; part of this work was done during visiting
Microsoft Research Asia.
??
Institute of Industrial Science, University of Tokyo2

varying transmittance; the variation corresponds to the ratio. Instead of using
such special equipment, some methods use a set of input images taken withdifferent exposures from a fixed viewpoint so that the ratio becomes available.
These methods are divided into parametric and non-parametric approaches.
In the parametric framework, Mann and Picard [3] propose a method thatassumes the response functions can be approximated by gamma correction functions. Mitsunaga and Nayar use a polynomial function for the representation [4].
Grossberg and Nayar apply principal component analysis (PCA) to a database ofreal-world response functions (DoRF) and show that the space of response functions can be represented by a small number of basis functions [5]. Mann [6] uses
the comparametric equation that defines the relationship between two images
taken with different exposure times.
Debevec and Malik estimate response functions with a non-parametric representation using a smoothness constraint [7]. Tsin et al . estimate non-parametricforms of response functions using the statistical analysis [8]. Pal et al . propose
to use a Bayesian network consisting of probabilistic imaging models and priormodels of response functions [9]. While non-parametric approaches have greater
descriptive power, the large number of unknown variables often lead to computationally expensive or unstable algorithms.
A few authors have proposed more general estimation methods that allow
camera motion or scene motion between input images. Grossberg and Nayar [10]estimate the function from images with scene motion using brightness histograms,
which avoid difficult pixelwise image registration. Kim and Pollefeys [11] propose
a method that allows free camera motion by finding the correspondence among
images. These methods can handle more general cases, but still require images
taken with multiple exposures.
Another class of the estimation methods is based on the physics of the imaging process. Lin et al . estimate the function from the color blending aroundedges in a single image [12]. For grayscale images, a 1D analogue of the 3D color
method is presented in [13]. Wilburn et al . propose to use the temporal blending of irradiance [14]. Matsushita and Lin propose an estimation method fromprofiles of image noise distributions [15]; their method relies on the assumptionthat profiles of noise distributions are symmetric. More recently, Takamatsu etal . [16] propose to use a probabilistic intensity similarity measure [17] for theestimation. Unlike these approaches, our method only uses noise variance, whichis considered a lower level of information.
1.2

Proposed method

In this paper, we propose a method for estimating camera response functions
from observation of image noise variance. Figure 1 illustrates the overview ofthe proposed method. As shown in the figure, noise variance in the input domain has an affine relationship with input intensity. This relationship, however,is not maintained in the output domain due to the non-linearity of camera response functions. By estimating a function that projects the non-affine intensityvariance relationship to an affine one, our method derives the response functions.

3

Observation

Imaging process
t
u
p
t
u
O

e
c
n
a
i
r
a
v
e
s
i

Response fRadiometric Calibration

t
u
p
n
I

e
c
n
a
i
r
a
v
e
s
i

Inverse response g

o
N

o
N

Input

Input

e
c
n
a
i
r
a
v
e
s
i

o
N

OutputOutputInput

Response f

Input: Noise variance in input

OutputResponse f

OutputOutputFig. 1. The noise variance in the input domain has an affine relationship with input
intensity level. Due to the non-linearity of the response function, the affine relationshipis lost in the output domain. The proposed method estimates the response functions
by recovering the affinity of the measured noise variances.
Response f

Input
Input
: Noise variance in output

Fig. 2. The relationship between response function and noise variances in input andoutput domains. The magnitude of output noise variance (height of the filled region)
varies with the slope of the response function with a fixed input noise variance.

Figure 2 depicts the relationship between the response function and noise variances in input and output domains. The ratio of the input/output noise varianceis closely related to the slope of the response function.
Unlike Matsushita and Lin’s work [15], our method does not rely on profilesof noise distributions which are more difficult to obtain in practice. The proposed
method has the following two advantages over their method. First, it uses themore reasonable assumption; the affinity of intensity-variance relationship in the
input domain. We derive this relationship from the nature of photon (electron)
noise [18, 19]. Second, the proposed method uses only noise variances rather thannoise profiles, therefore it only requires less amount of information about noise.
Practically, it makes the procedure of noise measurement easier.
This paper has two major contributions. It first theoretically derives how theradiometric response function alters the intensity-variance relationship, and viceversa. While it has been pointed out that these two quantities are related [8, 20],
the exact alternation has not been explicitly described. Second, it introduces
a method that has a wide range of applicability even with noisy observations.
While many existing algorithms break down when the noise level is high, weshow that the proposed method is not sensitive to the noise level because it uses
the noise as information.

2

Noise Variance and Response Functions

This section provides a theoretical relationship between the response function
and noise variances in the input and output domains. We show that the levelsof noise variance in the input and output domains are related by the square of
the first derivative of the camera response function. This result is later used todevelop the estimation algorithm in Section 3.

4

2.1

Derivation of the relationship between noise variance andresponse functions

We treat the relationship between input or output intensity and noise variance2
as the function similar to Liu et al . [20]. The noise variance function σO(Õ) inthe output domain, where the noise-free output intensity is Õ, can be described
as
Z
Z
2
˜
σO
(Õ) = (O − µO )2 p(O|Õ)dO = (f (I) − µO )2 p(I|I)dI.
The conditional density function (cdf) p(O|Õ) represents the noise distribution
in the output domain, i.e., the probability that the output intensity becomes O
˜ represents the noise
when the noise-free intensity is Õ. Likewise, the cdf p(I|I)˜ The
distribution in the input domain when the noise-free input intensity is I.
˜function f is the response function, and Õ and I˜ are related by f as Õ = f (I).
µO (= µO (Õ)) is the expectation of the output intensity with the cdf p(O|Õ).
Using the Taylor series of f (I) around I˜ and assuming that the second- andhigher-degree terms of the series are negligible (discussed in Section 2.3), we
obtain2
˜ I2 (I),˜σO(Õ) ' f 02 (I)σ(1)
2 ˜
where σI (I) is the function of noise variance in the input domain when the noise˜ This equation shows the relationship between the noise
free input intensity is I.
variance and the response functions. This relationship has been pointed out in [8,20], however, the correctness of the equation was not thoroughly discussed.
To simplify the notation in derivation of Eq. (1), we define:Z˜˜˜µI = µI (I) = Ip(I|I)dI,
Id = µI − I,Z˜ = (I − µI )n p(I|I)dI
˜
Mn = Mn (I)(n ≥ 0).
(2)
˜ Id denotesµI is the expectation, or mean, of input intensity I with the cdf p(I|I).
˜ Mn isthe difference between the expectation µI and the noise-free intensity I.
the n-th moment about the mean. Note that we do not put any assumptions
about the profile and model of the noise distributions.
The expectation of the output intensity µO where the noise-free output intensity is Õ can be written asZ∞X˜f (j) (I)˜,
(3)µO = µO (Õ) = f (I)p(I|I)dI= Õ +
Njj!
j=1
where Nj is the j-th moment about I˜ defined as follows.
ZNj =˜ j p(I|I)dI˜(I − I)=Z˜((I − µI ) + Id ) p(I|I)dI
=
j

j  
X
j j−k
I
Mk . (4)k dk=0Note that the response function f is represented using its Taylor series.

5
2
From the definition, the noise variance function σO(Õ) where the noise-free
output intensity is Õ can be derived as
!2
∞ X
∞
∞
X
X
˜
˜
˜ f (k) (I)f (j) (I)f (j) (I)2
Lj,j + 2
Lj,k ,(5)σO (Õ) =
j!
j!
k!
j=1
j=1
k>j

where Lj,k is defined asZ˜
Lj,k = ((I − µI )j − Nj )((I − µI )k − Nk )p(I|I)dI= Nj+k − Nj Nk .

(6)

As can be seen in the above equation, Lj,k is commutative (Lj,k = Lk,j ).
˜ L1,2 = 2Id σ 2 (I)˜ + M3 , · · · .
A detailed calculation gives us L1,1 = M2 = σI2 (I),I
Substituting these terms into Eq. (5), we obtain:2
˜ 1,1 + f 0 (I)f
˜ 00 (I)L˜ 1,2 + · · ·σO(Õ) = f 02 (I)L2 ˜
˜ 2 (I)˜ + f 0 (I)f
˜ 00 (I)(2I˜= f 02 (I)σd σ (I) + M3 ) + · · · .
I

I

(7)

Eq. (7) is the exact form of the relationship between the response function andnoise variances in the input and output domains. By discarding the second- andhigher-degree terms of Eq. (7), Eq. (1) is obtained. We discuss the validity ofthis approximation in Section 2.3.
2.2 Noise variance function
Eq. (1) shows the relationship between the response function and noise variancefunctions in the input and output domains. The input to our method is themeasured noise variance in the output domain. This section models the noise
˜ in the input domain so that the estimation algorithm
variance function σI2 (I)for the inverse response function g can be developed.
Input intensity I with camera noise can be written as
I = aP + NDC + NS + NR ,(8)

where a is a factor of photon-to-electron conversion efficiency with amplification,
and P is the number of photons. NDC , NS , and NR indicate dark current noise,
shot noise, and readout noise, respectively [21]3 . The noise-free input intensity I˜equals to aP .
Now we consider the noise variance function in the input domain. We assume
the different noise sources are independent. From Eq. (8), the noise variance˜ can be written asfunction in the input domain, σI2 (I),˜ = Iσ˜ 2 + σ2 + σ2 ,
σI2 (I)
S
DC
R

(9)

where σ∗2 denotes the variances of the different noise sources [8]. Eq. (9) can be
written in a simplified form as
˜ = AI˜ + B,
σI2 (I)(10)22
where A = σS2 and B = σDC+ σR. This equation clearly shows the affine˜relationship between the noise-free input intensity I˜ and the noise variance σI2 (I).
3

The effect of the fixed-pattern noise is included in the term P .

6

2.3

Validity of the approximation

It is important to consider the validity of the approximation in Eq. (7). In
this section, we show it in the following steps. First, we show that Li,j becomes
exponentially smaller as i+j increases. Second, the second largest term in Eq. (7),L1,2 , is small enough to be negligible compared with L1,1 through a detailed
calculation. Hereafter, we normalize the input and output intensity ranges from 0
to 1.
Relationship between Li,j and i + j By assuming independence of different noise
sources, the second- and higher-order moments can be computed by summingup the moments of noises from different sources. Three types of noise sourcesmust be considered: dark current noise, shot noise, and readout noise [21]. We
do not consider low-light conditions [22], so the effect of the dark current noise
becomes small.
The probability density of the readout noise can be considered to have a˜
normal distribution with a mean value equal to the noise-free input intensity I.
The moment can be written as0(i is odd)
NRi = Qi/2
(11)2
j=1 (2j − 1)σR (i is even).
2
 1, so the i-th moment of the readout noise NRi about
In Eq. (11), (2j − 1) σR
the noise-free input intensity becomes exponentially smaller as i increases.
Shot noise is modeled as Poisson distribution [21]. From the theory of generalized Poisson distribution [23], the moment MSi about the mean of the distribution is defined as
MSi ' ai−2 σS2 Ĩ + O(σS4 Ĩ ) (i ≥ 2),(12)

since the minimum unit of the distribution equals to a (See Eq. (8)). σS2 is theĨ
˜ By substituting
variance of shot noise where the noise-free input intensity is I.
this into Eq. (4) yields
NSi ' Idi +i  Xi i−j j−2 2
I a σSĨ .
j d
j=2

(13)


Even in the worst case where ji is overestimated as 2i , NSi becomes exponentially smaller, since Eq. (13) is rewritten as
NSi ≤ Idi +i
X(2Id )i−j (2a)j−2 (2σSĨ )2 ,

(14)

j=2

and we know that 2Id  1, 2a  1, and (2σSĨ )2  1. This equation shows that
NSi exponentially decreases as i increases.

7

The term Li,j is defined as Li,j = Ni+j − Ni Nj in Eq. (6). Because the i-th
moment of image noise Ni can be computed as the sum of the readout and shotnoise as Ni = NRi + NSi , it also becomes exponentially smaller as i increases.
From these results, we see that the term Li,j becomes exponentially smaller as
i + j increases.
Ratio of L1,1 to L1,2 Now we show that L1,2 is small enough to be negligible
compared with L1,1 . A detailed calculation gives us L1,2 = 2Id M2 + M3 . The
third moment of shot noise MS3 can be computed from Eq. (12). Also, the third
moment of readout noise can be obtained using Eq. (4) as
2
MR3 = NR3 − 3Id NR2 − Id3 = −3Id σR− Id3 .

(15)

From these results, the following equation is obtained:2
L1,2 = 2Id M2 − 3Id σR− Id3 + aσS2 Ĩ .

(16)2
+ σS2 , a  Id , if M2 6 Id2 , the order of L1,2 is roughly the
Since M2 ' σRĨsame as the order of Id M2 . Since Id is the difference between the noise-free
input intensity and the mean which can be naturally considered very small, it is
implausible to have cases where M2  Id2 .
From these results, the order of L1,2 is roughly equivalent to the order of
Id L1,1 , and Id is small because it is computed in the normalized input domain,e.g., in the order of 10−2 (' 1/28 ) in 8-bit image case. Therefore, L1,2 is about
10−2 times smaller than L1,1 .
To summarize, L1,2 is sufficiently small compared with L1,1 , and Li,j decreases exponentially as i + j increases. Also, because response functions are
˜ 6 f 00 (I).
˜ Therefore, Eq. (7) can be well approximated by Eq. (1).
smooth, f 0 (I)3Estimation Algorithm

This section designs an evaluation function for estimating inverse response functions g, using the result of the previous section.
3.1 Evaluation function
2
From Eqs. (1) and (10), the noise variance σO(Õ) of the output intensity O is222
˜ I2 (I)˜ + σQ˜ I˜ + B) + σQσO(Õ) ' f 02 (I)σ= f 02 (I)(A
.

(17)

2
σQ
is the variance of the quantization noise, which affects after applying the response function. Using the inverse response function g, Eq. (17) can be rewrittenas122(Ag(Õ) + B) + σQ
.
(18)
σO
(Õ) =0g (Õ)2
2
2
The variance of the quantization noise σQ
becomes σQ= l2 /12, where l is thequantization interval. Since its distribution is uniform, the following equation
holds:
Z 2l
Z l
1 2 2
l2
2
2
σQ =
x p(x)dx =x dx =.
(19)l − 2l
12
− 2l

8
2
In the following, we use σO(Õ) to represent the measured noise variance tom2
discriminate from the analytic form of the noise variance σO(Õ). Using Eq. (18)2
and the measured noise variances σO
(
Õ),
our
methodestimates
the inversem
response function g that minimizes the following evaluation function:
X
2
2
2
2
E1 (g; σO(Õ)) = minσO(Õ) − σO(Õ) .
(20)
m
m
A,B

Õ

Eq. (20) involves the estimation of A and B, which can be simply solved bylinear least square fitting, given g.
To make the algorithm robust against the measuring errors, namely the erroneous component in the measured noise, we use weighting factors. Eq. (20) ischanged to
2E2 (g; σO(Õ)) = min P
m
A,B

X
2
1
2
2
w(Õ) σO(Õ) − σO(Õ) ,
m
w(Õ)

(21)

Õ

where the weight function w(Õ) controls the reliability on the measured noise
2variance σO(Õ) at the intensity level Õ. We use a Cauchy distribution (Lorentzianmfunction) for computing the weight function w(Õ):w(Õ) =

1,e2 + ρ

(22)22(Õ). A damping factor ρ controls the(Õ) − σO
where e is defined as e = σOmrelationship between the difference e and weight w(Õ). As ρ becomes smaller,
the weight w(Õ) decreases more rapidly as the difference e increases.
We also add a smoothness constraint to the evaluation function, and theevaluation function becomes1
1 X 002E3 (g; σO(Õ)) = P 2
g (Õ)2 ,(23)
E2 + λs
m
n
σ
(
Õ)
Õ
Õ O
m

Õ

where nÕ is the number of possible noise-free output intensity levels, e.g., 256 in
8-bit case. λs isfactor that controls the effect of the smoothness
Pa regularization
2
constraint. 1/ Õ σO
(
Õ)
is
a
normalization
factor that makes E2 independent
mof the degree of noise level.
Our method estimates the inverse response function ĝ by minimizing Eq. (23)2given the measured noise variance σO
(Õ):
m

2
ĝ = argmin E3 g; σO(Õ) .
m

(24)

g

3.2

Representation of inverse response functions

To reduce the computational cost, we represent the inverse response functions
using a parametric model proposed by Grossberg and Nayar [5]. In their method,

9

principal component analysis (PCA) is performed on the database of real-world
response functions (DoRF) to obtain a small number of eigenvectors that canrepresent the space of the response functions. As done by Lin et al . [12, 13],
we compute the principal components of the inverse response functions using
the DoRF. Using the principal components, we represent the inverse response
function g as g = g0 +Hc, where g0 is the mean vector of all the inverse response
functions, H is a matrix in which a column vector represents an eigenvector, andc is a vector of PCA coefficients. Following Lin et al . [12, 13], we use the first
five eigenvectors. Using this representation, the number of unknown variables issignificantly decreased, e.g., from 256 to 5 in the case of 8-bit images.
3.3

Implementation

In our implementation, we set the damping factor ρ to the variance of the difference e in Eq. (22). The regularization factor λs is set to 5 × 10−7 from ourempirical observation. Minimization is performed in an alternating manner. We
perform the following steps until convergence:
1. minimize the evaluation function in Eq. (23) with fixing the weight function w(Õ)
2. recompute the values of the weight function w(Õ) using the current estimation result
We use the Nelder-Mead Simplex method [24] as the minimization algorithm
implemented in Matlab as a function fminsearch. The values of the weightfunction w(Õ) are set to one for every Õ at the beginning. During the experiments, we used five initial guesses for the inverse response function g as the
input to the algorithm. The converged result that minimizes the energy score is
finally taken as the global solution.

4

Experiments

We used two different setups to evaluate the performance of the proposed algorithm; one is with multiple images taken by a fixed video camera, the otheris using a single image. The two setups differ in the means for collecting noise
variance information.
4.1

Multiple-images caseIn this experiment, the measurements of noise variances are obtained by capturing multiple shots of a static scene from a fixed viewpoint with fixed camera
parameters. From multiple images, a histogram of output intensities is createdfor each pixel. From the histogram, the noise-free output intensity Õ is determined by taking the mode of the distribution, assuming that the noise-free intensity should correspond to the most frequently observed signal. The pixelwise
histograms are then merged together to form the histogram h(O, Õ) for each

10
1

n
o
i
t
c
n
u
F
e
s
n
o
p
s
e
R
e
s
r
e
v
n
I

1

Ground Truth
method

0.9

tu
pn
i
de
zli
a
rm
o
N
0.8

0.8

0.8

0.7

0.6

0.5

0.6

0.5

0.4

0.5

0.4

0.3

0.4

0.3

Max diff.: 0.052
RMSE : 0.025

0.2

0.1

0
0

0.1

0.2

-3

x 10

Normalized output
0.3

0.4

0.5

0.6

0.7

0.3

Max diff.: 0.037
RMSE : 0.022

0.2

0.1

0
0.8

0.9

1

0

0.1

0.2

-4

x 10

Normalized output
0.3

0.4

0.5

0.6

0.7

Estimation
Observation

1oN

Max diff.: 0.037
RMSE : 0.024

0.2

0.1

0
0.8

0.9

1

0

0.1

0.2

-3

x 10

Normalized output
0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.8

0.9

1

1

1.6

e
c
n
a
i
r
a
V
e
s
i

Our

tu
pn
i
de
zli
a
rm
o
N

0.7

0.6

Ground Truth
method

0.9

Our

tu
pn
i
de
zli
a
rm
o
N

0.7

1

Ground Truth
method0.9

Our

Estimation
Observation

1.4

Estimation
Observation

0.9

0.8

1.2

se
cn
ari
a
V

se
cn
ari
a
V

0.8

se
cn
air
a
V

0.7

1

0.6

0.6

0.5

0.8

0.4

0.6

0.4

0.3
0.4
0.2
0.2
0.2
0.1

0
0

0.1

0.2

Normalized output
0.3

0.4

0.5

0.6

0.7

0
0.8

0.9

(a) Sony DCR-TRV9E (B)

1

0

0.1

0.2

Normalized output
0.3

0.4

0.5

0.6

0.7

0
0.8

0.9

1

(b) Sony DCR-TRV900 NTSC (G)

0

0.1

0.2

Normalized output
0.3

0.4

0.5

0.6

0.7

(c) Sony DSR-PD190P (R)

Fig. 3. Results of our estimation method. Top row: comparison of inverse response
functions. Bottom row: measured noise variance and fitting result.
Table 1. Mean RMSE and disparity ofthe estimated inverse response functions interms of normalized input. Three different
scenes were used for each camera.
(a)

(b)

(c)

Fig. 4. Recorded scenes corresponding to
the results in Figure 3 (a-c).

CameraMean RMSE Disparity
A. DCR-TRV9E
0.026
0.053
B. DCR-TRV900
0.024
0.040
C. DSR-PD190P
0.033
0.055

output intensity level Õ. Finally, the noise distribution p(O|Õ) is computed bynormalizing the histogram h(O, Õ) as
h(O, Õ)p(O|Õ) = P.
O h(O, Õ)(25)

Results We used three different video cameras for this experiment: Sony DCRTRV9E (Camera A), Sony DCR-TRV900 NTSC (Camera B), and Sony DSRPD190P (Camera C). To obtain the ground truth of Camera C, we used Mitsunaga
and Nayar’s method [4], and the Macbeth color checker-based method [1], andcombined these results by taking the mean. For Camera A and B, we used only
the Macbeth color checker-based method [1] to obtain the ground truth because
the exposure setting was not available in these cameras. The results obtained by
the proposed method are compared with the ground truth curves.
Figure 3 shows the results of our algorithm. The top row shows the plot of
the estimated inverse response functions with the corresponding ground truth
curves. The bottom row shows the estimated and measured distributions ofnoise variances; the horizontal axis is the normalized output, and the vertical

11

RMSE vs. Number of samples

Our method
Matsushita's method

0.35
0.3
0.25

E
S 0.2
M
0.15
R
0.1

0.5

Max difference vs. Number of samples
Our method
Matsushita's method

e
c 0.4
n
e
r 0.3
e
f
f
i
d0.2
x
a
0.1

M

0.05
0

0
25

50

100

Number of images

299

25

50

100

Number of images

299

Fig. 5. Comparison between our method and Matsushita and Lin’s method [15]. Our
method uses noise variance, but not profiles of noise distributions. Our method works
well even when the sampling number is relatively small.

axis corresponds to the noise variance. Figure 4 shows the scenes used to obtainthese results.
Figure 3 (a) shows an estimation result using the blue channel of Camera A.
The maximum difference is 0.052 and the RMSE is 0.025 in terms of normalized
input. As shown in the bottom of (a), the noise variances in lower output levels
contain severe measured errors. Our algorithm is robust against such errors because of the use of adaptive weighting factors. Figure 3 (b) shows the result of
Camera B (green channel). The maximum difference is 0.037 and the RMSE is
0.022. Figure 3 (c) shows the estimation result of Camera C (red channel). The
input frames are obtained by setting the camera gain to 12 db which causes highnoise level. The maximum difference is 0.037 and the RMSE is 0.024.
Table 1 summarizes all the experimental results. For each camera, three different scenes are used. The algorithm is applied to RGB-channels independently,
therefore 9 datasets for each camera are used. Disparity represents the meanof maximum differences in normalized input. From these results, the proposed
method performs well even though the algorithm only uses the noise variance as
input.
4.2

Comparison with another noise-based estimation method

Figure 5 shows the comparison between our method and Matsushita and Lin’s
method [15]. Unlike other estimation methods, these two methods take noise as
input. We use Camera B for the comparison.
As shown in the result, the estimation results are equivalent when the numberof images is relatively large. However, Matsushita and Lin’s method breaks downwhen the number of samples becomes small, and our method shows significant
superiority. In statistics, it is known that variance of measured from samples’
variance is inversely proportional to the number of the samples. Therefore, themeasured variance becomes more stable than the profile of noise distribution
does, as the number of samples increases. In addition, Matsushita and Lin’s
symmetry criterion naturally requires large number of samples to make the noise
profiles smooth, while it does not hold in the lower number of samples in Figure 5.
These are why our method works well when the number of samples is relatively
small.

12

4.3

Single-image caseWe describe a single-image case where only one shot of the scene is available. In
this setup, the distribution of noise variances are collected from uniformly colored
image regions. However, the measured noise distribution is often insufficient todetermine the inverse response functions because the limited measurements do
not span the entire range of output levels. To better constrain the problem, we
use a prior model p(g) of the inverse response functions obtained from the DoRF
as done in [12] and [15].
Using the prior model p(g), the MAP (maximum a posteriori) estimation is2
performed by maximizing the cdf p(g|σO(Õ)) which represents the probabilitymof the inverse response function being g when the measured noise variances are2σO(Õ) asm22ĝ = argmax p(g|σO(Õ)) = argmax log p(σO(Õ)|g) + log p(g) .
(26)
m
m
g

g

2
(Õ)|g) is defined as
The likelihood p(σOm2
p(σO(Õ)|g) =m
12
exp −λp E3 (g; σO(Õ)) ,mZ

(27)

where Z is the normalization factor, and λp is a regularization coefficient that
determines the weight on the evaluation function E3 . We empirically set λp to
2 × 104 in the experiments. The prior model p(g) is formed using a multivariate
Gaussian mixture model asp(g) =KX

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

(28)

i=1

where N represents a normal distribution with mean µ and covariance matrix Σ,and αi is a weight factor. The prior model is obtained using the PCA coefficientsof the inverse response functions in the DoRF by applying the cross-entropy
method [25]. The number of normal distributions K is set to 5 in our experiments.
Results We used a Canon EOS-20D camera for the experiment. To obtain the
ground truth, we used Mitsunaga and Nayar’s method [4] using images takenwith different exposures. Since our focus is on estimating the inverse response
functions from the measured noise variances, we photographed a scene composedof relatively flat and uniformly colored surfaces, so that the noise variances canbe easily obtained. The left image in Figure 6 shows one of two scenes used forthe experiment. We photographed them five times each at six different camera
gains (ISO 100 ∼ 3200). We manually selected 21 homogeneous image regions
to obtain the noise variances as input. In total, we ran our estimation algorithm
60 times (= 2 scenes × 5 shots × 6 ISO levels) for each RGB color channel.
Figure 6 summarizes the results of estimation at different ISO levels. The
noise level increases with the ISO gain level, as shown by the cropped images
on the top. The results indicate that the estimation is unaffected by the greater
noise level. The mean RMSE is almost constant across the different ISO levels,
which verifies that our method is not sensitive to the noise level.

13

0.045

0.043

E
S

R

0.041

M
R
n
a0.039
e

G
B

M
0.037

ISO 100

ISO 200

ISO 400

ISO 800

ISO 1600

ISO 3200

Fig. 6. Relationship between the noise level and mean RMSE of the estimates. Left
image shows one of the photographed scenes. Top row shows magnification of a part ofthe image at different ISO levels. Bottom row shows the mean RMSE of RGB channels
at each ISO gain level, and demonstrates that our estimation method is independentof noise levels.

5

Conclusions

In this paper, we have proposed the method for estimating a radiometric response
function using noise variance, not noise distribution, as input. The relationshipbetween the radiometric response function and noise variances in input and output domains is explicitly derived, and this result is used to develop the estimation
algorithm. The experiments are performed for two different scenarios; one is withmultiple shots of the same scene, and the other is only from a single image. Theseexperiments quantitatively demonstrate the effectiveness of the proposed algorithm, especially its robustness against noise. With our method, either special
equipment or images taken with multiple exposures are not necessary.
Limitations It is better for our method that the measured noise variances covera wide range of intensity levels. Wider coverage provides more information to
the algorithm, so the problem becomes more constrained. This becomes an issue, particularly in the single-image case. In the single-image case, we used a
simple method to collect the noise variances, but more sophisticated methods
such like [20] can be used to obtain more accurate measurements that could
potentially cover a wider range of intensity levels.

Acknowledgement
The authors would like to thank Dr. Bennett Wilburn for his useful feedback on
this research.

References
1. Chang, Y.C., Reid, J.F.: Rgb calibration for color image analysis in machine vision.
IEEE Trans. on Image Processing 5 (1996) 1414–1422
2. Nayar, S.K., Mitsunaga, T.: High dynamic range imaging: Spatially varying pixel
exposures. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2000) 472–479

14
3. Mann, S., Picard, R.: Being ‘undigital’ with digital cameras: Extending dynamic
range by combining differently exposed pictures. In: Proc. of IS & T 48th Annual
Conf. (1995) 422 – 428
4. Mitsunaga, T., Nayar, S.K.: Radiometric self-calibration. In: Proc. of Comp. Vis.
and Patt. Recog. (CVPR). (1999) 374–380
5. Grossberg, M.D., Nayar, S.K.: What is the space of camera response functions?
In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2003) 602–609
6. Mann, S.: Comparametric equations with practical applications in quantigraphic
image processing. IEEE Trans. on Image Processing 9 (2000) 1389–1406
7. Debevec, P.E., Malik, J.: Recovering high dynamic range radiance maps fromphotographs. In: Proc. of ACM SIGGRAPH. (1997) 369–378
8. Tsin, Y., Ramesh, V., Kanade, T.: Statistical calibration of ccd imaging process.
In: Proc. of Int’l Conf. on Comp. Vis. (ICCV). (2001) 480–487
9. Pal, C., Szeliski, R., Uyttendale, M., Jojic, N.: Probability models for high dynamic
range imaging. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2004) 173–180
10. Grossberg, M.D., Nayar, S.K.: What can be known about the radiometric response
function from images? In: Proc. of European Conf. on Comp. Vis. (ECCV). (2002)
189–205
11. Kim, S.J., Pollefeys, M.: Radiometric alignment of image sequences. In: Proc. ofComp. Vis. and Patt. Recog. (CVPR). (2004) 645–651
12. Lin, S., Gu, J., Yamazaki, S., Shum, H.Y.: Radiometric calibration from a single
image. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2004) 938–945
13. Lin, S., Zhang, L.: Determining the radiometric response function from a single
grayscale image. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2005) 66–73
14. Wilburn, B., Xu, H., Matsushita, Y.: Radiometric calibration using temporal irradiance mixtures. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2008)
15. Matsushita, Y., Lin, S.: Radiometric calibration from noise distributions. In: Proc.
of Comp. Vis. and Patt. Recog. (CVPR). (2007)
16. Takamatsu, J., Matsushita, Y., Ikeuchi, K.: Estimating camera response functions
using probabilistic intensity similarity. In: Proc. of Comp. Vis. and Patt. Recog.
(CVPR). (2008)
17. Matsushita, Y., Lin, S.: A probabilistic intensity similarity measure based on noise
distributions. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2007)
18. Janesick, J.R.: Photon Transfer. SPIE Press (2007)
19. Schechner, Y.Y., Nayar, S.K., Belhumeur, P.N.: Multiplexing for optimal lighting.
IEEE Trans. on Patt. Anal. and Mach. Intell. 29 (2007) 1339–1354
20. Liu, C., Freeman, W.T., Szeliski, R., Kang, S.B.: Noise estimation from a single
image. In: Proc. of Comp. Vis. and Patt. Recog. (CVPR). (2006) 901–908
21. Healey, G.E., Kondepudy, R.: Radiometric ccd camera calibration and noise estimation. IEEE Trans. on Patt. Anal. and Mach. Intell. 16 (1994) 267 – 276
22. Alter, F., Matsushita, Y., Tang, X.: An intensity similarity measure in low-light
conditions. In: Proc. of European Conf. on Comp. Vis. (ECCV). (2006) 267–280
23. Consul, P.C.: Generalized Poisson Distributions: Properties and Applications. Marcel Dekker Inc, New York (1989)
24. Nelder, J.A., Mead, R.: A simplex method for function minimization. Computer
Journal 7 (1965) 308–312
25. Botev, Z., Kroese, D.: Global likelihood optimization via the cross-entropy methodwith an application to mixture models. In: Proc. of the 36th Conf. on Winter simul.
(2004) 529–535