Estimating Intrinsic Images from Image Sequenceswith Biased Illumination
Yasuyuki Matsushita1 , Stephen Lin1 , Sing Bing Kang2 , andHeung-Yeung Shum1
1

Microsoft Research Asia
3F, Beijing Sigma Center, No.49, Zhichun Road Haidian District,
Beijing 100080, China
2
Microsoft Research
One Microsoft Way, Redmond, Washington 98052-6399, U.S.A.
{yasumat, stevelin, sbkang, hshum}@microsoft.com

Abstract. We present a method for estimating intrinsic images from afixed-viewpoint image sequence captured under changing illumination directions. Previous work on this problem reduces the influence of shadowson reflectance images, but does not address shading effects which can significantly degrade reflectance image estimation under the typically biased
sampling of illumination directions. In this paper, we describe how biased
illumination sampling leads to biased estimates of reflectance image derivatives. To avoid the effects of illumination bias, we propose a solution that
explicitly models spatial and temporal constraints over the image sequence.
With this constraint network, our technique minimizes a regularization function that takes advantage of the biased image derivatives to yield reflectanceimages less influenced by shading.

1

Introduction

Variation in appearance caused by illumination changes has been a challenging
problem for many computer vision algorithms. For example, face recognition is complicated by the wide range of shadow and shading configurations a single face canexhibit, and image segmentation processes can be misled by the presence of shadowsand shading as well. Since image intensity arises from a product of reflectance andillumination, one approach for dealing with variable lighting is to decompose animage into a reflectance component and an illumination component [7], also known
as intrinsic images [1]. The reflectance image, free of illumination effects, can thenbe processed without consideration of shadows and shading.
Decomposition of an image into intrinsic images, however, is an underconstrained
problem, so previous approaches in this area introduced additional constraints tomake the problem tractable. In [6], it is assumed that the illumination componentis spatially smooth while the reflectance component exhibits sharp changes, suchthat low-pass filtering of the input image yields the illumination image. Similarly,[3] assumes smoothness of illumination and piecewise constant reflectance, so thatremoving large derivatives in the input image results in the illumination image. In
addition to illumination smoothness, Kimmel et al. [5] include constraints that thereflectance is smooth and the illumination image is close to the input image.

Instead of relying on smoothness constraints, Tappen et al. [10] proposed alearning-based approach to separate reflectance edges and illumination edges ina derivative image. Although this method successfully separates reflectance andshading for a given illumination direction used in training, it is difficult to createsuch a prior to classify edges under arbitrary lighting. Another edge-based method
was proposed by Finlayson et al. [11] that suppresses color derivatives along theillumination temperature direction to derive a shadow-free image of the scene. In
addition to shadow edges, this approach may remove texture edges that also havea color derivative in the illumination temperature direction.
Rather than using only a single input image, Weiss [9] deals with the simpler
scenario of having an image sequence captured from a fixed viewpoint with changing
illumination conditions. This method employs a ML estimation framework basedon a prior that illumination images yield a Laplacian distribution of derivatives
between adjacent pixels. Experimental results demonstrate that this technique efficiently and robustly removes cast shadows from reflectance images. Shading on
non-planar surfaces, however, can significantly degrade ML estimation of intrinsic
images by altering the distribution of derivatives, especially in the typical case of
a biased illumination distribution that is not centered around the surface normalsof the adjacent pixel pair. More recently, Matsushita et al. [12] extended Weiss’s
method to handle the scenes where the Lambertian assumption does not hold. Using the reflectance image estimated by ML estimation as a scene texture image,their method derives time-varying reflectance images instead of assuming a single
reflectance image.
In our proposed method, we also take as input an image sequence and analyze thederivative distributions. Because of the effects of illumination bias on the derivative
distributions, we present an alternative method for processing image derivatives,
based on explicit modeling of spatial and temporal constraints over the image sequence. With this constraint network, a reflectance image and a set of illumination
images are estimated by minimizing a function based on smoothness of illumination and reflectance. Although the derivative distributions are unsuitable for ML
estimation, our technique nevertheless takes advantage of derivative distribution
information to spatially vary the weight of the smoothness constraints in a manner
unlike previous regularization-based methods.
The goal of this work is closely related to that of photometric stereo with unknown light sources and spatially varying albedo. One strong assumption in mostuncalibrated photometric stereo approaches [19–21] is that there are no cast shadows. However, it is clear that this assumption does not hold in many situations forreal world scenes. Yuille et al. [22] have proposed a method to handle cast shadows
using robust statistics; however, one drawback of the method is that it assumes
a single point source in each image. Photometric stereo yields accurate results,
but generally it is necessary to assume limiting conditions. While the photometricstereo framework relies on the structural smoothness, our method relies more on thesmoothness of reflectance and illumination images. In the context of photometricstereo, Wolff and Angelopoulou [4] acquired multiple stereo pairs of images of thesame scene under different illumination conditions. With two stereo pairs they ob-tain a stereo pair of photometric ratio images, in which the albedo term is removed
in order to extend geometric stereo reconstruction to smooth featureless surfaces.
The remainder of the paper is organized as follows. Sec. 2 details the problemof illumination bias and the resulting effects of shading on derivative distributions.
In Sec. 3, we describe the constraints on the energy minimization process and theinfluence of the derivative distribution. Our algorithm is presented in Sec. 4, followedby experimental results in Sec. 5 and a conclusion in Sec. 6.

2Effect of Illumination Bias

Before describing the effect of illumination bias on derivative distributions, we begin
with a brief review of intrinsic images and the ML estimation technique.
Under the Lambertian assumption, as expressed in the following equation, an
input image I arises from a product of two intrinsic images: the reflectance image ρand the illumination image L. Since the viewpoint of the image sequence is fixed, ρdoes not vary with time t. The illumination is comprised by an ambient term α anda direct illumination term LD , which is the product of the illumination intensity
E, a binary shadow presence function g and the inner product of surface normal n
and illumination direction l:I(x, y, t) = ρ(x, y)L(x, y, t)= ρ(x, y) LD (x, y, t) + α(x, y, t)= ρ(x, y) E(t)g(x, y, t) n(x, y) · l(t) + α(x, y, t)= ρ(x, y)E(t) g(x, y, t) n(x, y) · l(t) + α0 (x, y, t)(1)

where n · l is always non-negative, and α0 indicates the ambient term normalized by
the illumination intensity E. In the ML estimation framework of [9], n derivative
filters fn are first applied to the logarithms of images I(t). For each filter, a filteredreflectance image ρn is then computed as the median value in time of fn ? log I,
where ? represents convolution:log ρ̂n (x, y) = mediant {fn ? log I(x, y, t)}.

(2)

The filtered illumination images log Ln (x, y, t) are then computed using the estimated filtered reflectance images log ρ̂n according tolog L̂n (x, y, t) = fn ? log I(x, y, t) − log ρ̂n (x, y).

(3)

Finally, a reflectance image ρ and illumination images L are recovered from thefiltered reflectance images ρn and illumination images Ln (t) through the following
deconvolution process,X(log ρ̂, log L̂) = h ?
fnr ? (log ρ̂n , log L̂n )(4)nwhere fnr is the reversed filter of fn , and h is the filter which satisfies the following
equation:
h?

X
n


fnr ? fn = δ.

(5)

n1

n2

n1

n2

(a)
(b)
Fig. 1. Illumination conditions. (a) Uniform illumination, (b) biased illumination.

From (2), it can be shown that for two adjacent pixels with intensities I1 (t) andI2 (t),I1 (t)
ρ1 E(t) g1 · n1 · l(t) + α01ρ̂n = mediant.
(6)= mediant ·I2 (t)
ρ2 E(t) g2 · n2 · l(t) + α02

We assume that α0 is constant over an image, i.e., α01 (t) = α02 (t). Cast shadows affectthis equation only when g1 6= g2 . Since this instance seldom occurs, cast shadows do
not affect the median derivative values used in ML estimation. It can furthermore
be seen that when n1 = n2 , shading does not affect ML estimation since n1 ·l = n2 ·land consequently ρ̂n = ρ1 /ρ2 .
When a pair of pixels have different surface normals, ML estimation can also
deal with shading in cases of unbiased illumination samples. For a pair of adjacentpixels with surface normals n1 and n2 , the set Ωl of illumination samples l(t) areunbiased only under the following condition:
(7)
medianl(t)∈Ωl n1 · l(t) − n2 · l(t) = 0.
In other words, the illumination is unbiased for a given pair of pixels when the illumination image value L(x, y) of both pixels is the same for the median derivativevalue. Otherwise, the illumination distribution is biased. Figure 1 shows an illustration of unbiased illumination and biased illumination for a given pair of pixels.
With unbiased illumination as given in (7), it can be seen that (6) results in thecorrect value ρ1 /ρ2 .
When a pair of adjacent pixels have different surface normals, illumination bias
will cause the ML estimation to be incorrect, because n1 · l 6= n2 · l for the median
derivative value. In this case, the illumination ratio in (6) does not equal to one, andconsequently ρ̂n 6= ρ1 /ρ2 . This can be expected since different shading is present in
the two pixels for every observation.
The case of different surface normals with illumination bias is a significant one,because for a pair of adjacent non-planar pixels, unbiased illumination is rare. So
for most pairs of non-planar pixels, ML estimation fails to compute the correctreflectance ratio and the estimated reflectance image contains shading. Figure 2
shows a typical result of ML estimation applied to a synthetic scene with nonplanar surfaces. A ball on a plane is lit from various directions as exemplified inFig. 2. Shading effect remains on reflectance estimate with ML estimation. Left: an input
image sequence, Right: the estimated reflectance image with ML estimation.

the input images on the left side of Figure 2. Although the illumination samples areunbiased for some pairs of pixels, they are biased for most pairs of adjacent pixels.
As a result, shading remains in the estimated reflectance image as shown on theright side of the figure.

3

Solution Constraints

Since ML estimation is generally affected by shading, we propose an alternative
solution method based on the constraints described in this section. Let us denotei, j as labels for illumination directions, p, q for adjacent pixels, and N, M for thenumber of observed illumination conditions and the number of pixels in an image,respectively.
From a sequence of images, we can derive spatial constraints between adjacentpixels (inter-pixel) and temporal constraints between corresponding pixels under
different light directions (inter-frame). Moreover, we employ smoothness constraints
to make the problem tractable.
[Inter-frame constraint]
Assuming that the scene is composed of Lambertian surfaces, the reflectancevalue at each point ρ is constant with respect to time. We can thereby derive a
temporal constraint as follows:
Lp (ti )Ip (ti )=,Ip (tj )Lp (tj )0 ≤ i, j < N ; i 6= j.

(8)This does not determine the absolute values of Ls; however, it fixes the ratios among
Lp s.
[Inter-pixel constraint]
Letting ωp be a set of pixels that are neighbours of p,
Ip (ti )ρp Lp (ti )=·,Iq (ti )ρq Lq (ti )0 ≤ i < N ; q ∈ ωp .

(9)

I1 (t n ) ρ1 L1 (t n )=⋅I 2 (t n ) ρ 2 L2 (t n )xy

I1 (t m ) L1 (t m )=I1 (t n ) L1 (t n )tFig. 3. Inter-frame and inter-pixel constraints. A set of constraints composes a constraint
network.

This constraint can possibly be applied to non-neighboring pixels, however, we
restrict this to be applied only to neighboring pixels because we use the flatness
constraint and smooth reflectance constraint in the energy minimization step.
These constraints can be derived from (1), and they compose a 3-D constrained
network about L and a 2-D constrained network about ρ as illustrated in Fig. 3.
We use them as hard constraints and force ρ and L to always satisfy the following
equation:!2!2
X
X
Ip (ti ) ρp Lp (ti )Ip (ti )Lp (ti )
+
= 0.
(10)
−
−
·
Ip (tj ) Lp (tj )Iq (ti ) ρq Lq (ti )
p,q,i;q∈ω
p,i,j;i6=jp[Smoothness constraints]
In addition to the above constraints, our technique favors smoothness over bothρ and L. Smoothness is a generic assumption underlying a wide range of physicalphenomena, because it characterizes coherence and homogeneity. Based on the factthat retinal images tend to be smooth according to natural image statistics [16], we
assume that both ρ and L are smooth as well. By formulating these two assumptions into an energy function, we derive intrinsic images by an energy minimization
scheme.
The choice of energy function E(ρ, L) is a critical issue. Various kinds of energyfunctions that measure the smoothness of observed data have been proposed. For
example, in regularization-based methods [13, 14], the energy minimization function
makes the observed data smooth everywhere. This generally yields poor results
at object boundaries since discontinuities should be preserved. One discontinuitypreserving function is Geman and Geman’s MRF-based function [15].
Although our method assumes smoothness of L, this condition clearly does nothold when the surface normals of adjacent pixels differ. In such instances, the smoothL constraint should be weakened. To estimate the amount of difference between
neighboring surface normals, we use the information present in the derivative distribution.
If a pair of pixels lie on a flat surface, the values of Ip (t)/Iq (t) are almost always
equal to 1 except when only one of the pixels is shadowed, as discussed in [9]. We
use this strong statistic and define an error function based on the hypothesis that
Ip and Iq share the same surface normal:
 

 
Ip
Ip
(11)
− arctanepq (ti ) = arctan mediant
Iq
Iq
In (11), mediant (Ip /Iq ) corresponds to ML estimation, which gives the ratio ofreflectance if p and q are co-planar. We evaluate the angle between the ratio ofreflectance and the ratio of observed intensity to determine if the observation supports the flatness hypothesis. To determine the amount of support for the flatness
hypothesis, a threshold  is used:
1 epq (ti ) <  : acceptξpq (ti ) =(12)0 epq (ti ) ≥  : reject
Finally, we compute the ratio of the number of acceptances ξ to the number of total
observations N to test the surface-flatness hypothesis. When the surface flatness f
is high, it is likely that L is smooth. On the other hand, when f is not high, thesmoothness assumption for L should be weakened. We define the flatness f as thesquare of the acceptance ratio:
fpq =

P

i ξpq (ti )

N

2

.

(13)

Using the surface flatness evaluated by (13), we define an energy function EΩ
to minimize:XEΩ =Ep (∆ρp , ∆Lp (t))
p

=

X X
p q∈ωp(ρp − ρq )2 + λfpq (ti ) Lp (ti ) − Lq (ti )

2

(14)

where λ is a coefficient that balances the smoothness of ρ and L.
Equation (14) always converges to a unique solution because EΩ is convex withrespect to ρp and Lp . This is confirmed by taking Ep ’s Hessian Hp
#
 "P
 2
1
0
∂ Ep /∂ρ2 ∂ 2 Ep /∂ρ∂L
q∈ωp
P
=
Hp =
,
(15)
0
λ q∈ωp fpq∂ 2 Ep /∂L∂ρ ∂ 2 Ep /∂ L2

the leading principal minors of which are 1 > 0, λfpq > 0 where λ > 0, fpq > 0, sothat the function Ep is strictly convex.
P Since the sum of convex functions is convex,
EΩ is also convex because EΩ = p Ep .

4

Algorithm

With the constraints described in the preceding section, our algorithm can proceed
as follows.
[Step 1 : Initialization] Initialize ρ and L.
[Step 2 : Hard constraints] Apply the inter-frame and inter-pixel constraints
expressed in (10). Since it is difficult to minimize the two terms in (10) simultaneously, we employ an iterative approach for minimization.
1. Inter-frame constraint. Update Lp (ti ).
X Ip (ti )Lp (ti ) ←Lp (tj ) /(N − 1)
Ip (tj )

(16)

j6=i

2. Inter-pixel constraint. Update Lp (ti ) and ρp with ratio error β. Letting Mωp be
the number of p’s neighboring pixels,X
Ip (ti ) ρq Lq (ti )/Mωp .
(17)·
βp (ti ) =Iq (ti ) ρp Lp (ti )q∈ωpSince the error ratio βp (ti ) can be caused by some unknown combination of ρand L, we distribute the error ratio equally to both ρ and L in (18) and (20),respectively.
q(18)
Lp (ti ) ← βp (ti )Lp (ti ),
βp =

X
i

ρp ←


βp (ti ) /N,p

βp ρ p .

(19)
(20)

3. Return to 1. unless Equation (10) is satisfied.
[Step 3 : Energy minimization]
Evaluate the energy function (14), and find ρ and L that lower the total energy.
If the total energy can still decrease, update ρ and L, then go back to Step 2. Otherwise, we stop the iteration. By fixing ρq and Lq in (14), the energy minimization is
performed for each Ep using the conjugate gradient method. The conjugate gradientmethod is an iterative method for solving linear systems of equations which converge faster than the steepest descent method. For further details of the algorithm,
readers may refer to a well presented review [18].

5Experimental Results

To evaluate our method, we carried out experiments over one synthetic image sequence and three real world image sequences. In these experiments, we used 11
different lighting conditions, and set λ = 0.4. We used  = 0.02 which is empiricallyobtained in Equation (12) for all experiments. The determination of  has dependency on minimum signal-to-noise ratio. Starting with constant initial values, ρ and
L are iteratively updated. There is no restriction about the initial values, however,
we used flat images for their initial values because of the smoothness assumption.

Fig. 4. Input image samples from synthetic scene. Illumination samples are chosen to bebiased.

Fig. 5. Estimated reflectance images. Left : our method, Center : the ground truth, Right: ML estimation.

5.1

Synthetic scene

For a synthetic scene, we prepared a Lambertian scene with a sphere and a plane
as shown in Figure 4. The illumination samples are biased in most cases, since mostof them lie on the left-hand side of the images. Figure 5 shows the result of our
method, the ground truth, and the result of ML estimation from left to right. Dueto the scaling ambiguity of the reflectance images, we adjusted the scaling of eachreflectance image for better comparison and display. As we can clearly see in the
figure, our method can successfully derive shading-free a reflectance image which is
close to the ground truth.
5.2

Real world scenes

For real world scenes, we captured two image sequences of toy scenes and used
Yale Face Database B [17]. Figure 6 and Figure 8 show the result of reflectanceestimation from Lambertian scenes. In both figures, the left image shows the resultof our method, while the right image shows the result of ML estimation. As we
can see clearly in them, our method handles shading more correctly and shading
effect is much reduced in our reflectance estimates. Figure 7 and 9 show the estimated illumination images and corresponding input images. In illumination images,
reflectance edges such as texture edges are well removed.
On the other hand, Figure 10 shows a negative result, especially on the hair.
Since the human hair shows high specularity, and it is hard to model it as Lambertian. This non-Lambertian property affects to our method and turns those area
into white. This is because our method is based on Lambertian model, and it implies that our method does not handle specular reflections well. However, as for the

(a)
(b)
Fig. 6. Toy scene 1. Estimated reflectance images. (a) Our method, (b) ML estimation.

(a)(b)
Fig. 7. (a) Estimated illumination image, and (b) the corresponding input image.

non-specular part such as the human face, shading effect is much reduced by our
method.

6

Conclusion

We have presented a method that robustly estimates intrinsic images by energy minimization. Unlike previous methods, our method is not affected by illumination bias
which generally exists. In our framework, we explicitly modeled spatial and temporal constraints over the image sequence to form a constraint network. Using this as a
hard constraint, we minimized an energy function defined from the assumptions thatreflectance and illumination are smooth. By weighting these smoothness constraints
according to a surface flatness measure estimated from derivative distributions, weestimated intrinsic images with improved handling of shading. Evaluation with bothsynthetic and real world image sequences shows that our method can robustly estimate shading-free reflectance image and illumination images. Some of the next
steps of our research will include the acceleration of energy minimization part andextension of our model to correctly handle specularity.

References
1. H.G. Barrow and J.M. Tenenbaum: Recovering intrinsic scene characteristics from images. In A. Hanson and E. Riseman, editors, Computer Vision Systems. Academic

(a)
(b)
Fig. 8. Toy scene 2. Estimated reflectance images. (a) Our method, (b) ML estimation.

(a)(b)
Fig. 9. (a) Estimated illumination image, and (b) the corresponding original image.
Press, 3–26, 1978.
2. E.H. Adelson and A.P. Pentland: The Perception of Shading and Reflectance. In D.
Knill and W. Richards (eds.), Perception as Bayesian Inference, 409–423, 1996.
3. A. Blake: Boundary conditions of lightness computation in mondrian world. In Computer Vision, Graphics and Image Processing , 32, 314–327, 1985.
4. L. B. Wolff and E. Angelopoulou: 3-d stereo using photometric ratios. In Proc. of theThird European Conference on Computer Vision (ECCV), pp 247–258, 1994.
5. R. Kimmel, M. Elad, D. Shaked, R. Keshet and I. Sobel: A Variational Framework forRetinex. In International Journal of Computer Vision, 52(1), 7–23, 2003.
6. E.H. Land: An alternative technique for the computation of the designor in the Retinex
theory of color vision. In Proc. Nat. Acad. Sci. , 83, 3078–3080, Dec. 1986.
7. E.H. Land: The Retinex theory of color vision. In Scientific American , 237(G), No. 6,
108–128, Dec. 1977.
8. E.H. Land, and J.J. McCann: Lightness and retinex theory. In Journal of the OpticalSociety of America, 61(1), 1–11, 1971.
9. Y. Weiss: Deriving intrinsic images from image sequences. In Proc. of 9th IEEE Int’l
Conf. on Computer Vision, 68–75, Jul., 2001.
10. M.F. Tappen, W.T. Freeman, and E.H. Adelson: Recovering Intrinsic Images from aSingle Image. In Advances in Neural Information Processing Systems 15 (NIPS), MIT
Press, 2002.
11. G.D. Finlayson, S.D. Hordley, and M.S. Drew: Removing Shadows from Images. In
Proc. of European Conf. on Computer Vision Vol.4, 823–836, 2002.

(a)
(b)
Fig. 10. Reflectance images estimated from Yale Face Database B. (a) Our method, (b)ML estimation. The result shows the limitation of our method, i.e. high specularity affects
the result.
12. Y. Matsushita, K. Nishino, K. Ikeuchi and M. Sakauchi : Illumination Normalization
with Time-dependent Intrinsic Images for Video Surveillance. In Conf. on Computer
Vision and Pattern Recognition (CVPR), Vol.1, pp. 3–10, 2003.
13. A. Blake: Comparison of the efficiency of deterministic and stochastic algorithms forvisual reconstruction. In IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(1):2–12, 1989.
14. B.K.P. Horn and B. Schunk: Determining optical flow. In Artificial Intelligence,17:185–
203, 1981.
15. S. Geman and D. Geman: Stochastic relaxation, Gibbs distributions, and the Bayesian
restoration of images. In IEEE Transactions on Pattern Analysis and Machine Intelligence, 6:721–741, 1984.
16. B.A. Olshausen and D.J. Field: Emergence of simplecell receptive field properties bylearning a sparse code for natural images. In Nature, 381:607-608, 1996.
17. A.S. Georghiades, P.N. Belhumeur, and D.J. Kriegman: From few to many: Generativemodels for recognition under variable pose and illumination. In IEEE Int. Conf. on
Automatic Face and Gesture Recognition, 277-284, 2000.
18. J.R. Shewchuck: An introduction to the conjugate gradient method without agonizing
pain. Tech. Rep. CMU-CS-94-125, Carnegie Mellon University, 1994.
19. H. Hayakawa: Photometric stereo under a light-source with arbitrary motion. In Journal of Optical Society of America A., 11(11):3079–3089, 1994.
20. R. Basri and D. Jacobs: Photometric stereo with general, unknown lighting. In Proc.
of Computer Vision and Pattern Recognition(CVPR), Vol.2, pp. 374–381, 2001.
21. A.S. Georghiades, P.N. Belhumeur and D.J. Kriegman: Illumination-Based Image Synthesis: Creating Novel Images of Human Faces Under Differing Pose and Lighting. In
Proc. Workshop on Multi-View Modeling and Analysis of Visual Scenes, pp. 47–54,
1999.
22. A.L. Yuille, D. Snow, R. Epstein, P. Belhumeur: Determining Generative Models for
Objects Under Varying Illumination: Shape and Albedo from Multiple Images Using
SVD and Integrability. In International Journal on Computer Vision., 35(3), pp 203–
222. 1999.