IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 26, NO. 10,

OCTOBER 2004

1

Illumination Normalization with
Time-Dependent Intrinsic Images
for Video Surveillance
Yasuyuki Matsushita, Member, IEEE, Ko Nishino, Member, IEEE,
Katsushi Ikeuchi, Fellow, IEEE, and Masao Sakauchi, Member, IEEE
Abstract—Variation in illumination conditions caused by weather, time of day, etc., makes the task difficult when building video
surveillance systems of real world scenes. Especially, cast shadows produce troublesome effects, typically for object tracking from afixed viewpoint, since it yields appearance variations of objects depending on whether they are inside or outside the shadow. In thispaper, we handle such appearance variations by removing shadows in the image sequence. This can be considered as apreprocessing stage which leads to robust video surveillance. To achieve this, we propose a framework based on the idea of intrinsic
images. Unlike previous methods of deriving intrinsic images, we derive time-varying reflectance images and corresponding
illumination images from a sequence of images instead of assuming a single reflectance image. Using obtained illumination images,
we normalize the input image sequence in terms of incident lighting distribution to eliminate shadowing effects. We also propose anillumination normalization scheme which can potentially run in real time, utilizing the illumination eigenspace, which captures theillumination variation due to weather, time of day, etc., and a shadow interpolation method based on shadow hulls. This paper
describes the theory of the framework with simulation results and shows its effectiveness with object tracking results on real scene
data sets.
Index Terms—Intrinsic images, reflectance, shadow removal, illumination normarization, video surveillance, robust tracking.

æ
1

INTRODUCTIONV

IDEO surveillance systems involving object detection
and tracking require robustness against illumination
changes caused by variation of, for instance, weather
conditions. Difficult situations arise not only by the changeof illumination conditions, but also by the large castshadows of surrounding structures, i.e., large buildings
and trees. Since most visual tracking algorithms rely on theappearance of the target object, typically using color,
texture, and feature points as cues, these shadows degrade
the quality of tracking. In urban scenes, where building
robust traffic monitoring systems is of special interest, it isusual to have large shadows cast by tall buildings
surrounding the road, e.g., Fig. 2a. Building a robust video
surveillance system under such an environment is achallenging task. To make the system insensitive todramatic change of illumination conditions and robust
against large static cast shadows, it would be valuable tocancel out those illumination effects from the image. Y. Matsushita is with Microsoft Research Asia, 3/F, Beijing Sigma Center,
No. 49, Zhichun Road, Hai Dian District, Beijing, China 100080.
E-mail: yasumat@microsoft.com.
. K. Nishino is with the Department of Computer Science, Columbia
University, 1214 Amsterdam Avenue, MC 0401 New York, NY 10027.
E-mail: kon@cs.columbia.edu.
. K. Ikeuchi and M. Sakauchi are with the Institute of Industrial Science, theUniversity of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan.
E-mail: ki@cvl.iis.u-tokyo.ac.jp, sakauchi@sak.iis.u-tokyo.ac.jp.
Manuscript received 1 July 2003; revised 22 Mar. 2004; accepted 1 Apr. 2004.
Recommended for acceptance by R. Basri.
For information on obtaining reprints of this article, please send e-mail to:tpami@computer.org, and reference IEEECS Log Number TPAMI-0151-0703.
0162-8828/04/$20.00 ß 2004 IEEE

sequence. Our goal is to “normalize” the input image
sequence in terms of the distribution of incident lighting to
remove illumination effects including shadow effects. We
should note that our method does not consider shadows
cast by moving objects but those cast by static objects suchas buildings and trees. To achieve this goal, we propose anapproach based on intrinsic images. Our method is composed of two parts as shown in Fig. 1.
The first part is the estimation of intrinsic images, which
is an offline process, depicted in A of Fig. 1. In this part,first, the scene background image sequence is estimated toremove moving objects from the input image sequence.
Using this background image sequence, we then derive
intrinsic images using our estimation method which isextended from Weiss’s ML estimation method [21]. Using
the estimated illumination image, which is a part of intrinsic
images, we are able to robustly cancel out the illumination
effects from input images of the same scene, enabling many
vision algorithms such as tracking to run robustly. After thederivation, we construct a database using Principal Component Analysis (PCA), which we refer to as illumination
eigenspace. This linear subspace captures the variation oflighting conditions in the illumination images. The subspace is used for the following direct estimation method.
The second part is a direct estimation of illumination
images, shown in B of Fig. 1. Using the preconstructed
illumination eigenspace, we estimate an illumination image
directly from an input image. To obtain accurate illumination images, shadow interpolation using shadow hulls is
accomplished.
Published by the IEEE Computer Society

2IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

Fig. 1. System diagram for illumination-normalization.

In the remainder of this paper, we first overview related
work in Section 1.1 and, in Section 2, we propose a method
to derive time-varying reflectance images Rðx; y; tÞ andcorresponding illumination images Lðx; y; tÞ. Derivation of
the illumination-invariant images using the illumination
images Lðx; y; tÞ is described in Section 3. In Section 4, we
propose a method to estimate R and L directly from an
input image using PCA. In addition, to obtain moreaccurate illumination images, we use a shadow hull-based
interpolation method, which is described in Section 4.1.
Experimental results are described in Section 5. Finally, we
conclude the paper in Section 6.

1.1 Related Work
Barrow and Tenenbaum proposed to consider every retinal
image as a composition of a set of latent images, which theyrefer to as intrinsic images [7]. One type of the intrinsic
images, R, contains the reflectance values of the scene,
while the other type, L, contains the illumination intensities,
and their relationship can be described by I ¼ R  L. Since
illumination images L represent the distribution of incident
lighting onto the scene while reflectance images R depict

VOL. 26,

NO. 10,OCTOBER 2004

the surface reflectance properties of the scene, thisrepresentation becomes useful to analyze and manipulate
the reflectance/lighting properties of the captured scene.
While decomposing a single image into intrinsic images,
namely, a reflectance image and an illumination image,remains a difficult problem [5], [7], [4], deriving intrinsic
images from image sequences has seen great success.
Recently, Weiss developed an ML estimation framework[21] to estimate a single reflectance image and multiple
illumination images from a series of images captured from afixed view point but under significant lighting condition
variation. Finlayson et al. [6] proposed a similar approach toours independently. They derive the scene texture edges
from the lighting-invariant image, and by subtracting thoseedges from the raw input images, they successfully derive
shadow-free images of the scene. We also take advantage ofthe fact that the reflectance image essentially models thescene texture in a manner invariant to lighting conditions.
We accomplish edge substitution between the reflectanceimage and illumination images, enabling robust derivationof scene-texture-free illumination images. In a single image
framework, Tappen et al. [17] proposed an approach to use
learning-based shading classifier accompanied with chromaticity-based classification method. Categorization was
done by assuming that differences in chromaticity values
indicate reflectance edges and then classifying edges having
similar chromaticity values by using a shading patternclassifier [13] which is prelearnt using synthetic images.
Although this method successfully separates shading andreflectance components, the shading classifier must accountfor the illumination direction, which is generally difficult todetermine at each local point.
Several other works on shadow detection have been
proposed. Deterministic model-based approaches to detectshadow regions are proposed by Kilger [14] and Koller et al.
[3] that exploit gray level, local, and static features. In
statistical approaches, Stauder et al. [9] and Jiang et al. [2]proposed a nonparametric approaches independently thatuse color, global, and dynamic features for enhancing object
detection.

Fig. 2. (a) An input image iðx; y; tÞ, (b) reflectance image derived by ML estimation rw ðx; yÞ, (c) illumination image derived by ML estimation lw ðx; y; tÞ,(d) our time-varying reflectance image rðx; y; tÞ, and (e) our illumination image lðx; y; tÞ.

MATSUSHITA ET AL.: ILLUMINATION NORMALIZATION WITH TIME-DEPENDENT INTRINSIC IMAGES FOR VIDEO SURVEILLANCE

2

INTRINSIC IMAGE ESTIMATION

Our work is inspired by ML estimation method [21] of
estimating intrinsic images from a sequence of images
captured from a fixed viewpoint under varying illumination. In this method, natural image statistics [8] are used aspriors to estimate a single reflectance image and multiple
illumination images. Although the method assumes thescene to be Lambertian, it is still effective to extract thescene texture, not a reflectance image, from a scene whichcontains specular surfaces. In that case, however, theestimated reflectance image does not accurately depict
surface reflectance property of the scene and some portionof the reflectance property will be included in illumination
images. If we consider real world scenes, it is often difficult
to expect the Lambertian assumption to hold. A typicalexample is white lines on a road surface, which show
variable reflection with respect to illumination changes. We
prepared a synthetic scene shown in Fig. 2a. The scene iscomposed of a road surface and pedestrian stripes that havedifferent reflectance property from the road surface. Withan assumption that the view point is fixed and the scene isroughly Lambertian, the ML estimation-based method
derives the scene reflectance image by taking the medianof filtered input image in log domain. Denoting log
reflectance of the road surface and stripes with r1 and r2 ,respectively, the reflectance difference at the boundary of
the road surface and the stripes becomes jr1  r2 j. Since theML estimation method assumes that the scene is Lambertian, i.e., r1 and r2 are constants, the difference jr1  r2 j isalso constant. However, if the Lambertian assumption does
not hold, e.g., r2 is varying, jr1  r2 j is no longer constant. In
that case, the error fjr1  r2 ðtÞj  mediant jr1  r2 ðtÞjg appears in illumination images as a ghost of scene texture asshown in Fig. 2c.
To handle this problem, we propose to assume a set oftime-varying reflectance images Rðx; y; tÞ instead of a timeinvariant reflectance image Rðx; yÞ. To derive the timevarying reflectance images, we start with estimating a scenetexture image using the ML estimation method [21]. We
denote the scene texture image which is the reflectanceimage estimated by the ML estimation method and theillumination image with subscript w, i.e., Rw and Lw , andour reflectance image and illumination image, R and L,
respectively. Applying the ML estimation method, a single
reflectance image Rw ðx; yÞ and a set of illumination images
Lw ðx; y; tÞ are estimated. Our goal is to derive time-varying,i.e., lighting condition dependent, reflectance images
Rðx; y; tÞ, and corresponding illumination images Lðx; y; tÞthat do not contain scene texture.
Iðx; y; tÞ ¼ Rðx; y; tÞ  Lðx; y; tÞð1Þ

We use lower-case letters to denote variables in log domain,e.g., r represents the logarithm of R. With nth derivative
filters fn , a filtered reflectance image rwn is computed bytaking median along the time axis of fn ? iðx; y; tÞ, where ?
represents convolution. We used two derivative filters, i.e.,f0 ¼ ½0 1  1 and f1 ¼ ½0 1  1T . With those filters, input
images are decomposed into intrinsic images by Weiss’s

3

method as described in (2). The method is based on thestatistics of natural images [8].
r^wn ðx; yÞ ¼ mediant ffn ? iðx; y; tÞg:ð2ÞThe filtered illumination images lwn ðx; y; tÞ are then computed by using estimated filtered reflectance image rwn .
l^wn ðx; y; tÞ ¼ fn ? iðx; y; tÞ  r^wn ðx; yÞ:ð3Þ

To be precise, l is computed by l ¼ i  r in the unfiltereddomain in Weiss’s original work while we estimate lin the derivative domain for the following edge-based
manipulation.
We use the output of the ML estimation method as initialvalues of our intrinsic image estimation. As mentioned
above, the goal of our method is to derive time-dependent
reflectance images Rðx; y; tÞ and their corresponding illumination images Lðx; y; tÞ. The basic idea of the method is toestimate time-varying reflectance components by cancellingthe scene texture from initial illumination images. To factor
out the scene textures from the illumination images andassociate them with reflectance images, we use the texture
edges of rw . We take a straightforward approach to remove
texture edges from lw and derive illumination images
lðx; y; tÞ with (4) and (5). These equations describe that if
the magnitude of a pixel in the gradient of the texture
image, i.e., jrwn ðx; yÞj, is larger than a threshold T , then thisis an evidence of a texture edge. In that case, the texture
edges that appear in the gradient of the illumination image
are removed from the illumination image, and the removededge is added to the time-varying reflectance image. Here,
we assume that there is no scene texture edge which has notappeared in rwn . Since the scene texture edges are the edges
which stay at the same position, it is reasonable to assume
that they are captured in rwn .

0 ifjrwn ðx; yÞj > Tð4Þ
ln ðx; y; tÞ ¼
lwn ðx; y; tÞ otherwise;
rn ðx; y; tÞ ¼

rwn ðx; yÞ þ lwn ðx; y; tÞ if jrwn j > T
rwn ðx; yÞ otherwise;

ð5Þ

where T represents a threshold value. While we currently
manually set the threshold value T used to detect texture
edges in rwn , we found the procedure is not so sensitive tothe threshold as long as it covers texture edges well.
Since the operation is linear, the following equation isimmediately confirmed.
fn ? iðx; y; tÞ ¼ rwn ðx; yÞ þ lwn ðx; y; tÞ¼ rn ðx; y; tÞ þ ln ðx; y; tÞ:ð6ÞFinally, time-varying reflectance images rðx; y; tÞ and scenetexture-free illumination images lðx; y; tÞ are recovered fromfiltered reflectance rn and illumination images ln throughthe following deconvolution process as done in [21].
!
X
r
^
^
f ? ð^
rn ; ln Þ ;ð7Þð^r; lÞ ¼ g ?
nn

where fnr is the reversed filter of fn , and g is the filter which
satisfies the following equation:

4IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

g?

X

VOL. 26,

NO. 10,OCTOBER 2004!
fnr ? fn

¼ :

ð8Þ

n

Reconstruction of images from edges has been studied
for a long time. The problem can be formulated as aboundary value problem, and an alternative technique thathas been widely accepted is the multigrid method. Multigrid methods solve the boundary value problem efficiently
by changing the problem into a set of linear algebraic
equations. Multigrid methods can be used to reconstruct
images from edges as an alternative to the method
described above. Readers are referred to a good overview
on multigrid solvers by Press et al. [20].
To demonstrate the effectiveness of our method forderiving time-dependent intrinsic images, we rendered a
CG scene which contains cast shadows and surface patcheswith different reflectance properties, which is analogous toreal road surfaces, e.g., pedestrian stripes. Fig. 2 shows aside-by-side comparison of the results of applying the ML
estimation method and our method. The first row is the CG
scene, where the scene has the property that the histogramof derivative-filtered output is sparse, which is the required
property of the ML estimation-based decomposition method and also is the statistics usually found in natural images.
As can be seen clearly, texture edges are successfully
removed from our illumination image while they obviously
remain in illumination image derived by the ML estimation
method. Considering an illumination image to be an image
which represents the distribution of incident lighting, ourillumination image is much better since incident lighting has
nothing to do with the scene reflectance properties.
Fig. 3. An input image I (a) and the illuminance-invariant image N (b).

3

SHADOW REMOVAL

Using the obtained scene illumination images by our
method, the input image sequence can be normalized interms of illumination.
To estimate the intrinsic images of the scene wherevideo surveillance systems are to be applied, it is necessary
to remove moving objects from the input image sequence
because our method requires the scene to be static.
Therefore, we first create background images in each shorttime range (T ) in the input image sequence, assumingthat the illumination does not vary in that short time
period. We simply use the median image of the short input
sequence as the background image, but of course the more
complicated methods would give the better backgroundimages [11]. The assumption here is that moving objects inthe scene are not observed at the same point longer thanthe background in T . These background images Bðx; y; tÞare used for the estimation of intrinsic images. Using theestimation method described in the former section, each
image in the background image sequence is decomposed
into corresponding reflectance images Rðx; y; tÞ and illumination images Lðx; y; tÞ.
Bðx; y; tÞ ¼ Rðx; y; tÞ  Lðx; y; tÞ:ð9ÞOnce decomposed into intrinsic images, any image whose
illumination condition is captured in the series of Bðx; y; tÞ
can be normalized with regards to its illumination condition
by simply dividing the input image Iðx; y; tÞ by its

corresponding estimated illumination image Lðx; y; tÞ.
Through the normalization, cast shadows are also removed
from the input images.
Since the incident lighting effects are fully captured inillumination images Lðx; y; tÞ, the normalization by dividing
with L corresponds to removing the incident lighting
distribution from the input image sequence. Let us denote
the resulting illuminance-invariant image with Nðx; y; tÞ,that can be derived by the following equation:
Nðx; y; tÞ ¼ Iðx; y; tÞ=Lðx; y; tÞ:ð10ÞFig. 3 shows the result of our normalization method.
Fig. 3a shows the input image I and Fig. 3b represents theilluminance-invariant image N. Note shadows of thebuildings are removed in N.
Since we consider the time-dependent reflectancevalues, accurate normalization of incident lighting can bedone compared to using illumination images that arederived by ML estimation method. Fig. 4 depicts thedifference between using our illumination image L(Figs. 4a and 4c) and using illumination image Lw by
ML estimation (Figs. 4b and 4d). As can be seen in theimages, the pedestrian stripes appear over the trucks inevery right-hand side image. On the other hand, whenusing our illumination images, the ghost of pedestrian
stripes is almost vanished as seen in every left-hand side
image. This is because our method handles reflectanceMATSUSHITA ET AL.: ILLUMINATION NORMALIZATION WITH TIME-DEPENDENT INTRINSIC IMAGES FOR VIDEO SURVEILLANCE

5

widely used in signal processing, statistics, and neuralcomputing. The reason why we employed PCA is that weassume the observed illumination images have Gaussian
distribution. The basic idea in PCA is to find the basic
components ½s1 ; s2 ; . . . ; sn  that explain the maximum
amount of variance possible by n linearly transformed
components. Fig. 6 shows the hyperplane constructed by
mapping illumination images onto the eigenspace visualized with the first three eigenvectors.
In our case, we mapped Lw ðx; y; tÞ into the illumination
eigenspace, instead of mapping Lðx; y; tÞ. This is because,
when given an input image, the reflectance image
Rw ðx; yÞ is useful to eliminate the scene texture bycomputing Iðx; y; tÞ=Rw ðx; yÞ, and the resulting image
becomes Lw ðx; y; tÞ. We keep the mapping between
Lw ðx; y; tÞ and Lðx; y; tÞ to derive final Lðx; y; tÞ estimates.
First, an illumination space matrix is constructed bysubtractingLw , which is the average of all Lw , i.e.,PLw ¼ n1 n Lw , from each Lw and stacked column-wise.
P ¼ fLw1  Lw ; Lw2  Lw ; ; Lwn  Lw g:

ð11Þ

P is an N  M matrix, where N is the number of pixels in
the illumination image and M is the number of illumination
images Lw . We made the covariance matrix Q of P as
follows:
Q ¼ PPT :Fig. 4. The difference of the normalization results between using ourillumination images L and illumination image Lw by ML estimation. The
left images of each pair, (a) and (b), show the results using ourillumination image L, and the right ones, (c) and (d), are the results using
illumination image Lw derived by ML estimation.

variation properly while the reflectance values are fixed in
the ML estimation method.

4

ILLUMINATION EIGENSPACE FOR DIRECT
ESTIMATION OF ILLUMINATION IMAGES

The intrinsic image estimation method described in theformer section is fully offline since the method basicallyrequires estimation of the background images and taking
the median over the accumulated images to obtain the final
intrinsic images. Computation of the intrinsic images using
the method described in Section 2 is fast for typical framesizes of current videos. However, it is necessary to firstestimate background images to adopt the method. Therefore, for our use the computation of the intrinsic images isdone offline. However, realtime processing is required forpractical use. In this section, we describe our approach torealtime derivation of illumination images for shadow
removal. Our method first stores a lot of illumination
images captured under different illumination conditions.
Using stored illumination images, realtime estimation ofillumination image from an input image is accomplished.
We propose illumination eigenspace to model variations ofillumination images of the scene. The illumination eigenspace is an eigenspace into which only illumination effects
are transformed. We use principal component analysis(PCA) to construct the illumination eigenspace of a target
scene, in our case, the crossroad shown in Fig. 5. PCA isð12Þ

Finally, the eigenvectors ei and the corresponding eigenvalues i of Q are determined by solving,i ei ¼ Qei :ð13Þ

To solve (13), we utilized Turk and Pentland’s method [16],which is useful to compute the eigenvectors when thedimension of Q is high. Fig. 6 shows the illumination
eigenspace on which all the illumination image, Lw ðx; y; tÞ,of 2,048 images from 120 days (7:00 - 15:00) are mapped.
Using the illumination eigenspace, direct estimation ofan illumination image can be done given an input image
which contains moving objects. We consider that the global
similarity of the illumination image is measured by the
distance weighted by the contribution ratio of eigenvectors(eigenvalues) in the illumination eigenspace. Thus, we firstdivide the input image by a reflectance image to get a
pseudoillumination image L which includes dynamic
objects: L ¼ Iðx; y; tÞ=Rw ðx; yÞ. For the division, we use
the scene texture image Rw as a reflectance image. Then, byusing this pseudoillumination image as a query, best
approximation of the corresponding illumination image L^
is estimated from the illumination eigenspace based on thefollowing distance metric,ﬃ
X qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
L^w ¼ arg min8Lwi
wj ðF ðL ; jÞ  F ðLwi ; jÞÞ2 ;ð14Þ
jwhere F is a projection function which maps an illuminaPtion image onto the jth eigenvector, and wj ¼ j =  i
where i s are the eigenvalues. The term F ðLwi Þ ispremapped onto the illumination eigenspace, each Lwi is apoint in the eigenspace. To achieve the minimization of (14),
we used nearest neighbor (NN) search in the illumination

6

IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

VOL. 26,

NO. 10,OCTOBER 2004

Fig. 5. Direct estimation of intrinsic images. Each row shows different weather conditions. (a) An input image I, (b) the pseudoillumination image Lw ,(c) the estimated illumination image L^w by the NN search in the illumination eigenspace, and (d) the corresponding background image B to (c).

eigenspace with the assumption that the distribution of thepseudoillumination image is also Gaussian. We also assume
that the moving objects are small enough compared to theimage size and the error of Rw is small. For NN search, theSR-tree method [18] is used which is known to be fast

Fig. 6. Illumination eigenspace constructed using 2,048 images from
120 days data of a crossroad.

especially for high-dimensional and nonuniform data
structures such as natural images. Finally, the illumination
image Lðx; y; tÞ is derived using the mapping table from L^w
to L. The number of stored images for this experiment was
2,048 and the contribution ratio was 84.5 percent at 13
dimensions, 90.0 percent at 23 dimensions, and 99.0 percent
at 120 dimensions. We chose to use 99.0 percent of
eigenratio for this experiment. The compression ratio wasapproximately 17:1, and the disk space needed to store thesubspace was about 32 MBytes when the image size is
320  243.
Results of illumination image search is shown in Fig. 5.
In this figure, starting with the left-hand side column, thefirst column shows input images I, the second column
shows pseudoillumination images L , and the third column
corresponds to estimated illumination images L^w . The rightend column shows the background images which correspond to the estimated illumination images. The NN searchin PCA is reasonably robust to estimate the most similar
illumination image Lw from the pseudoillumination image
L . However, since the sampling of the illumination images
is sparse, there are slight differences in the shadow shapes.
It is possible to acquire the exact illumination image L when
the database is dense enough, but it is not easy to preparesuch a database. To solve this problem, we propose anapproach to compute intermediate illumination images byinterpolating shadow regions using geometry estimatedMATSUSHITA ET AL.: ILLUMINATION NORMALIZATION WITH TIME-DEPENDENT INTRINSIC IMAGES FOR VIDEO SURVEILLANCETABLE 1
Dimension of the Illumination Eigenspace, Contribution Ratio,
and NN Search Cost

from sampled cast shadow regions and sunlight angles. The
details are described in Section 4.1.
As for the computational cost, the average time of the
NN search is shown in Table 1 with MIPS R12000 300MHz,
when the number of stored illumination images is 2,048, theimage size is 360  243 and the number of output search
results is 5. Since the input image is obtained at the intervalof 33ms (at 30 frames/sec), the estimation time is fastenough for realtime processing.

4.1 Shadow Interpolation Using Shadow Hulls
Unfortunately, it is difficult to store all illumination images
under every possible illumination condition. NN search in
the illumination eigenspace gives good results, however, it
is often the case that they are slightly different from the true
illumination image due to the limitation of the number ofstored illumination images. Therefore, we propose tointerpolate NN search results to generate the final estimateof the illumination image. We assume global intensitychanges are linear as long as they are densely sampled, but
the motion of cast shadows cannot be represented by linear
intensity interpolation. To achieve the interpolation of the
shadow motion, we propose to use Shadow Hull s. The ideaof shadow hull is the same as Visual Hull (VH) except using
cast shadow silhouette instead of object’s silhouette andillumination source as a view point in VH. VH, also known
as Shape-From-Silhouette (SFS) [15], [1], is a popular
3D reconstruction method which estimates the shape ofthe object from multiple silhouette images. These traditional
approaches basically assume that the object is static.
Recently, extending the traditional SFS formulation tohandle the shape of a rigidly moving object over time anda dynamic articulated object are well investigated [12], [10].
Our purpose is to reconstruct the rough shape of an
object (such as a building) which gives enough information
to compute the cast shadow regions given illumination
directions. Even though parts of many objects cannot be
represented by an SFS-based technique, such as concavities,
it does not matter for our method because our purpose is
not reconstructing the 3D scene geometry but generatingintermediate shadow regions. To accomplish shadow hullreconstruction, we use estimated shadow region maps S,
sunlight angles which are immediately computed from thetime-stamp, and calibration parameters. Assuming that thescene is roughly planar, we pitch the shadow volume using
the shadow map S and illumination direction in the world
coordinate system. Taking the intersection of different
shadow volumes pitched using different shadow maps,
the shadow hull can be estimated. The resulting hull is notnecessarily precise, but it gives enough scene geometry tocompute cast shadow region between sampled illumination

7

conditions. Also, we do not expect this method to work well
if the sampling of the illumination variation is not
sufficiently dense.
To compute the intermediate illumination images, wefirst compute the shadow region map S for each stored
illumination image. Since an illumination image depicts theintensity distribution of reflected lighting from the scene,
we can assume that low intensity pixels in the illumination
image represent shadowed area. By thresholding, shadow
regions S are derived from the illumination images. To
automatically determine the threshold for shadowed area
determination, we adopt the clustering technique of Otsu
[19]. The threshold is computed from maximizing the
between-class scatter by minimizing within-class variances.
In our case, we classify pixels into two classes, shadowed
and lit, assuming shaded pixels can be categorized intoeither of them. We first create an intensity histogram foreach illumination image to obtain the probability density
function pðiÞ where i indicates an intensity value. Assumingthat shadowed pixels have relatively lower intensity than lit
pixels, we define the cumulative probability functions P for
shadowed (Ps ) and lit (Pl ) area using a threshold value T .
Ps ðT Þ ¼

T
X

pðiÞ;

Pl ðT Þ ¼

i¼imin

imax
X

pðiÞ:

ð15Þ

i¼T

In the same manner, we define the mean of the shadowed(s ) and lit (l ) area as functions of the threshold T as
follows:
s ðT Þ ¼

T
X
i¼iminipðiÞ;l ðT Þ ¼

imax
XipðiÞ:

ð16Þ

i¼T

Finally, the optimum threshold value Topt can be obtainedby the following equation:
Topt ¼ arg maxfPs ðT Þ  Pl ðT Þ  ðs ðT Þ  l ðT ÞÞ2 g:
8T

ð17Þ

Using Topt , we compute shadowed area for each
illumination image and store the shadowed area in shadow
map S. The shadow map S is a Boolean map which stores 1
if the pixel is shadowed, otherwise 0.
Shadow regions are then mapped onto the world
coordinate, and shadow volumes are computed using
shadow regions associated with sunlight angles. By taking
the intersection of shadow volumes in the 3D space, we get
the rough geometry of objects casting shadows, which hasenough information for computing intermediate castshadows (Fig. 7a). Fig. 7b shows the results of shadowinterpolation in a CG scene using an estimated shadow hull.
The dark regions show the interpolated shadow regions,while the lighter regions represent the sampled shadow
regions.
Shadow interpolation using shadow hulls is useful toestimate the intermediate shadow shapes between sampled
lighting conditions. Fig. 8 shows the interpolated result of
the real world scene. The left-hand column represents theestimated shadow regions, while the right-hand column
shows the corresponding illumination images. The top andbottom row represent the images under sampled illumination conditions, and the middle row depicts the interpolated

8IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

Fig. 7. Interpolation of cast shadow using a shadow hull. (a) Computing
shadow hulls using shadow regions associated with sunlight angles.
(b) Result of shadow interpolation.

result. With the assumption that illumination angle ismonotonously changing, intermediate shadow regions canbe approximately computed using a shadow hull pitchedby only two neighboring shadow map samples.
For computational efficiency, we adopt a multiple planarslice representation of shadow hull as shown in Fig. 9. For
each z along the vertical axis, slices parallel to the x-y plane
are assumed. Cast shadows on the world plane are thenprojected onto these slices using the illumination direction
to obtain the region where a shadow volume penetrates
each slice. Finally, the intersections of the projected shadow
regions are taken as a shadow hull region on the slice. These
intersections can be computed with simple Boolean operations. The reason why a set of 2D representation is adopted
instead of 3D representation is that our objective is togenerate intermediate shadow regions and, for that purpose, dense recovery of the shadow hull is not required. The
sampling rate along the z-axis is currently manually set,
however, we have confirmed that as long as it is not too
sparse virtually the same result is obtained. As describedabove, intermediate shadow shapes can be computed withsome projection and intersections computations that arevery efficient.
Once we obtain the intermediate shadow map Sint ðx; yÞ,
we can compute the intermediate illumination image Lint
by interpolating the first k illumination images Lk that are

VOL. 26,

NO. 10,OCTOBER 2004

Fig. 8. Shadow hull-based shadow interpolation. Figures in the top andbottom rows are shadow regions and sampled illumination images. The
middle row shows the interpolated shadow region. Note that theinterpolated result describes only the area of the shadow and is used
to compute the intermediate illumination image. The grid is overlaid forbetter visualization.

computed by NN search. With the kth illumination image
obtained by NN search Lk and wk which is the weighting
factor, the distance from Lw to Lwk (See Section 4) in theillumination eigenspace, the intermediate illumination
image Lint can be computed as follows:

Fig. 9. Illustration of shadow hull construction. To efficiently represent
the shadow hull, a set of 2D slices of the shadow hull along the verticalaxis is used instead of a volumetric representation.

MATSUSHITA ET AL.: ILLUMINATION NORMALIZATION WITH TIME-DEPENDENT INTRINSIC IMAGES FOR VIDEO SURVEILLANCE

9

Fig. 10. Comparison with the ground truth. (a) Interpolated result using our method, (b) the ground truth, (c) result of image differencing between (a)and (b).

P

wk Sk ðx; yÞLk ðx; yÞLint ðx; yÞ ¼ Sint ðx; yÞ k P
k wk Sk ðx; yÞP
wk ð1  Sk ðx; yÞÞLk ðx; yÞ;þ ð1  S int ðx; yÞÞ k P
k wk ð1  S k ðx; yÞÞ

ð18Þ

where Sk ðx; yÞ is a map which holds 1 if the pixel ðx; yÞ in Lk is
inside the shadow region, otherwise 0. As shown in (18),pixels in shadow and pixels outside the shadow are
separately interpolated. Practically, sometimes it is possiblePP
to observe that k wk ð1S k ðx; yÞÞ ¼ 0 or k wk ðS k ðx; yÞÞ ¼ 0.
In that case, the term whose denominator equals zero istreated as zero because the term does not contain any
information.
The resulting intermediate illumination image is shownin middle right in Fig. 8. Fig. 10 shows the comparison
between the result of our method and the ground truth. We
can notice the slight difference between them from Fig. 10c;
however, it shows a globally correct shadow shape which isuseful to remove shadow effects from the input image.
To evaluate the effectiveness of our shadow interpolation
method, we compared the results with simple image
interpolation as shown in Fig. 11. Both our result and image
interpolation result are computed from two neighboringillumination images. As we can clearly see in the figure, ourmethod generates quite similar cast shadow shape to theground truth while the simple interpolation yields a big
ghost at the shadow transition area. To quantify thedifference from the ground truth, we evaluated the absoluteintensity difference per pixel. Table 2 shows the average ofthe absolute difference in 12 sets of interpolation results. In
our shadow interpolation method, the error is much smaller
compared to image interpolation.

5

EXPERIMENTAL RESULTS

We evaluated our shadow elimination method by object
tracking based on block matching using two different
matching functions. One is using sum of squared differences (SSD) as a matching function, and the other is using
normalized correlation function (NCF) instead of SSD. NCF
is considered to be robust to illumination effects because thematching score is largely independent of linear variations ofshading on objects. We chose the block matching algorithm
for actual tracking to show even the simplest and widelyused tracking method can achieve good results after
utilizing our illumination normalization preprocess. The
block matching method based on SSD is accomplished by

Fig. 11. Comparison of the ground truth, our interpolation result, and image interpolation result. (a) The ground truth, (b) interpolation result using our
method, and (c) the result using simple image interpolation. Both (b) and (c) are computed from two illumination images.

10IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

TABLE 2Absolute Differences of Pixel Intensity from the Ground Truth

The top column shows the mean difference across the whole image andthe second column shows the mean difference only at the shadowtransition regions.

pursuing the most similar window in the neighboring frame
evaluated by (19).
(
M1
X N1
X
ft ðx þ m; y þ nÞSSDðx; yÞ ¼ min
8i;jm¼0 n¼0

 ft1 ðx þ m þ i; y þ n þ jÞ

2

)

ð19Þ
:

For the case of using NCF, the most similar window isfound by evaluating (20).
NCF ðx; yÞ ¼ max
8i;jPM1 PN1
m¼0
n¼0 ft ðx þ m; y þ nÞft1 ðx þ m þ i; y þ n þ jÞ
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ:
PM1 PN12
m¼0
n¼0 ft1 ðx þ m þ i; y þ n þ jÞ
ð20Þ
For the experiments, we chose image sequences containing frames of vehicles crossing boundaries of cast shadowssince we focus especially on the advantage of our shadow
elimination. The parameter set of vehicle tracking for eachsequence is totally even, i.e., the same initial window
position, 10  10 pixels of window, 10 pixels for maximal
search distance, and 10fps of frame-rate.

VOL. 26,

NO. 10,OCTOBER 2004

The result is shown in Fig. 12. In Fig. 12, the time axis isrepresented from top to bottom. The first column of eachpair and the second column of that represent results of theblock matching-based tracking applied to the original image
sequence and image sequence with our preprocessing of
illuminance normalization, respectively. In the original
image sequence in Figs. 12a, 12b, and 12c, we get typicalerroneous results where the block matching fails at shadowboundaries, because of the large intensity variation between
inside and outside the shadow. On the other hand, after
proper illumination normalization using our method, we
get successful results.
We accomplished the tracking experiments over 502
vehicles in 11 sequences under different lighting conditions.
Since we cannot have the ground truth, we carefully
evaluated the tracking results. The results are shown in
Table 3 and Table 4. In the tables, Ocorrect , Oerror , Ncorrect , andNerror mean:
Ocorrect : Count of the correct tracking results on theoriginal input image sequences.
. Oerror : Error count of the tracking results on theoriginal input image sequences.
. Ncorrect : Count of the correct tracking results on the
preprocessed image sequences using our method.
. Nerror : Error count of the tracking results on the
preprocessed image sequences using our method.
The success rate using the block matching method
based on SSD over the original input image sequences was
55.6 percent, while with normalized input it improved to69.3 percent as shown in Table 3. The effectiveness of our
method is clearly confirmed by the result that 45.3 percent
of originally failed results were rescued by our method.
On the other hand, 11.5 percent got worse after applying
our method. This poor effect happens typically when the.

Fig. 12. Results of tracking based on block matching. Along row from top to bottom it shows the frame sequence. The first column of each pair, (a),(b), and (c), shows the tracking result over the original image sequence, and the second column of each pair, (d), (e), and (f), shows thecorresponding result after our preprocessing.

MATSUSHITA ET AL.: ILLUMINATION NORMALIZATION WITH TIME-DEPENDENT INTRINSIC IMAGES FOR VIDEO SURVEILLANCE

11

TABLE 3
Tracking Result over 502 Sequences Using a Block Matching
Method Based on SSD

TABLE 4
Tracking Result over 535 Sequences Using a Block Matching
Method Based on NCF

shadow-edge cast on the vehicle surface largely differsfrom the shadow-edge in the illumination image. It
happens because our method currently can handle only
2D shadows on the image plane, but the actual shadow iscast three-dimensionally on the scene. When the gap of the
shadow-edge position is large, the error of the normalization gets large, as a result, the block matching fails. Our
system currently does not handle this problem since theerror rate is small compared to the improved correct rate,
but we are investigating on handling cast shadows threedimensionally.
Table 4 shows the tracking result using NCF. NCF
basically is not sensitive to overall illumination variations,however, it is vulnerable to the case where a shadowboundary is observed inside the window. Suppressing suchcases, the tracking error rate reduces from 32 percent to25.2 percent using our illumination normalization technique. The improvement in the tracking accuracy is not verysignificant using NCF, however, it is confirmed that ourillumination normalization method raises the accuracy ofSSD-based tracking to the same level of that of NCF-based
tracking. Also, it is important to note that our preprocess
does not affect the NCF-based method.
From the experimental results, we confirmed that our
method significantly improves the accuracy of tracking
using SSD-based block matching. The SSD-based method is
the simplest method in the class of tracking methods that do
not handle the appearance variation caused by illumination
effects. Actually, many rich matching algorithms use
SSD-based matching as their core component due to its
simplicity for implementation and low computational cost.
Therefore, we consider that our preprocess method can beused with this class of methods to improve the accuracy ofthe matching. As for the methods that are naturally robust
against the appearance variations caused by illumination
effects such as NCF-based methods, our method may notcontribute much to improve the accuracy. However, wehave confirmed that our method does not affect the NCFbased method at all and even slightly improves theaccuracy.

extended the method to properly handle surfaces with
nonrigid reflectance properties. This is accomplished bymodifying pseudoderivative illumination images withregards to the scene texture edges that can be derived
from the pseudoreflectance image estimated throughML estimation algorithm. The proposed method is remarkably robust and does not require the information of scene
geometry, camera parameter, and lighting condition at all,but requires the camera to be fixed and several lighting
conditions to be observed. As a key component of ourframework, we proposed to utilize illumination eigenspace, apreconstructed database which captures the illumination
variation of the target scene, to directly estimate illumination images for elimination of lighting effects of the scene,including elimination of cast shadows. As for the intermediate illumination images that cannot be represented bylinear combinations of sampled illumination images, we
proposed to use shadow hulls for interpolating cast shadow
regions using sunlight angles and estimated camera
parameters.
The effectiveness of the proposed method was shown bycomparing the tracking results between the original image
sequence and the image sequence with our method used asa preprocess. Since our method is used as a preprocessing
stage, we believe this method can be applied to many video
surveillance systems to increase the robustness against
lighting variations. The SSD-based matching method is
widely used as the core component of many tracking
algorithms and have been already integrated in manysystems. We have shown that our method can be used toincrease the reliability of such SSD-based matching methods. We also have shown that our preprocessing method
does not contribute much to the matching algorithms such
as NCF-based methods that are naturally robust against
appearance variation caused by illumination. However, ourmethod does not negatively affect NCF-based tracking
results and it even slightly improves the accuracy. Therefore, we conclude that our preprocessing method iseffective for the class of matching methods that do nothandle the appearance variation caused by illumination bythemselves which we strongly believe has significant use.
In addition to the illumination estimation using illumination eigenspace, we have investigated direct estimation ofillumination images corresponding to real scene images
using the illumination eigenspace and shadow interpolation
based on shadow hulls. Though our current implementation of the shadow interpolation in research code is not fastenough for realtime processing, we believe the frameworkhas the potential to be processed in realtime.

6

CONCLUSIONS

We have described a framework for normalizing illumination effects of real world scenes, which can be effectively
used as a preprocess for robust video surveillance. We
believe it provides a firm basis to improve existing
monitoring systems. We started with a previous method
to derive intrinsic images from image sequences, and12IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE,

REFERENCES
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]

[11]
[12]

[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]

A. Laurentini, “The Visual Hull Concept for Silhouette-Based
Image Understanding,” IEEE Trans. Pattern Analysis and MachineIntelligence, vol. 16, no. 2, pp. 150-162, Feb. 1994.
C. Jiang and M.O. Ward, “Shadow Identification,” Proc. IEEE Int’l
Conf. Computer Vision and Pattern Recognition, pp. 606-612, 1992.
D. Koller, K. Daniilidis, and H.H. Nagel, “Model-Based Object
Tracking in Monocular Image Sequences of Road Traffic Scenes,”
Int’l J. Computer Vision, vol. 10, pp. 257-281, 1993.
E.H. Adelson and A.P. Pentland, “The Perception of Shading andReflectance,” Perception as Bayesian Inference, pp. 409-423, 1996.
E.H. Land, “The Retinex Theory of Color Vision,” Scientific Am.,
vol. 237G, no. 6, pp. 108-128, Dec. 1977.
G.D. Finlayson, S.D. Hordley, and M.S. Drew, “Removing
Shadows from Images,” Proc. European Conf. Computer Vision,
vol. 4, pp. 823-836, 2002.
H.G. Barrow and J.M. Tenenbaum, “Recovering Intrinsic Scene
Characteristics from Images,” Computer Vision Systems, pp. 3-26,
1978.
J. Huang and D. Mumford, “Statistics of Natural Images andModels,” Proc. Conf. Computer Vision and Pattern Recognition,
pp. 541-547, 1999.
J. Stauder, R. Mech, and J. Ostermann, “Detection of Moving Cast
Shadows for Object Segmentation,” IEEE Trans. Multimedia, vol. 1,
no. 1, pp. 65-76, Mar. 1999.
K. Grauman, G. Shakhnarovich, and T. Darrell, “A BayesianApproach to Image-Based Visual Hull Reconstruction,” Proc. IEEE
Conf. Computer Vision and Pattern Recognition, vol. 1, pp. 187-194,
2003.
K. Toyama, J. Krumm, B. Brumitt, and B. Meyers, “Wallflower:
Principles and Practice of Background Maintenance,” Proc. Int’l
Conf. Computer Vision, pp. 255-261, 1999.
K.M. Cheung, S. Baker, and T. Kanade, “Shape-from-Silhouette of
Articulated Objects and Its Use for Human Body Kinematics
Estimation and Motion Capture,” Proc. IEEE Conf. Computer Vision
and Pattern Recognition, vol. 1, pp. 77-84, 2003.
M. Bell and W.T. Freeman, “Learning Local Evidence for Shading
and Reflection,” Proc. Int’l Conf. Computer Vision, vol. 1, pp. 670677, 2001.
M. Kilger, “A Shadow Handler in a Video-Based Real-Time Traffic
Monitoring System,” Proc. IEEE Workshop Applications of Computer
Vision, pp. 11-18, 1992.
M. Potmesil, “Generating Octree Models of 3D Objects from Their
Silhouettes in a Sequence of Images,” Computer Vision Graphics andImage Processing, vol. 40, pp. 1-20, 1987.
M. Turk and A. Pentland, “Eigenfaces for Recognition,” The J.
Cognitive Neuroscience, vol. 3, no. 1, pp. 71-86, 1991.
M.F. Tappen, W.T. Freeman, and E.H. Adelson, “Recovering
Intrinsic Images from a Single Image,” Advances in Neural
Information Processing Systems 15, MIT Press, 2002.
N. Katayama and S. Satoh, “The SR-Tree: An Index Structure forHigh-Dimensional Nearest Neighbor Queries,” Proc. 1997 ACM
SIGMOD Int’l Conf. Management of Data, pp. 369-380, 1997.
N. Otsu, “A Threshold Selection Method from Gray Level
Histograms,” IEEE Trans. Systems, Man, and Cybernetics, vol. 9,
pp. 62-66, 1979.
W.H. Press, S.A. Tuekolsky, W.T. Vetterling, and B.P. Flannery,
Numerical Recipes in C: The Art of Scientific Computing, second ed.
Cambridge Univ. Press, 1992.
Y. Weiss, “Deriving Intrinsic Images from Image Sequences,” Proc.
Ninth IEEE Int’l Conf. Computer Vision, pp. 68-75, July 2001

VOL. 26,

NO. 10,OCTOBER 2004

Yasuyuki Matsushita received the BEng,
MEng, and PhD degrees in electrical engineering from the University of Tokyo in 1998, 2000,
and 2003, respectively. Currently, he is anassociate researcher at Microsoft Research
Asia. His research interests include photometricmethods in computer vision, image-based rendering, and modeling from images. He is amember of the IEEE.

Ko Nishino received the BE and ME degrees
from The University of Tokyo in 1997 and 1999,
respectively. He received the PhD degree in
information science from The University ofTokyo in 2002. Since 2002, he has been a
postdoctoral research scientist at Columbia
University. His research interests span computer
vision and computer graphics. He has publishedseveral papers on photometric and geometricproblems in scene and object modeling that
include physics-based vision, image-based modeling, and rendering. He
is a member of the IEEE and ACM.
Katsushi Ikeuchi received the BE degree from
Kyoto University in 1973 and the PhD degreefrom the University of Tokyo in 1978. After
working at the MIT AI Laboratory for three years,
ETL for five years, and CMU Robotics Institutefor 10 years, he joined the University of Tokyo in
1996, and is currently a full professor. His
research interest spans computer vision, robotics, and computer graphics. In these research fields, he has received several awards,
including the David Marr Prize in computational vision for the paper“Shape from Interreflection,” and IEEE R&A K-S Fu memorial besttransaction paper award for the paper “Toward Automatic Robot
Instruction from Perception—Mapping Human Grasps to Manipulator
Grasps.” In addition, in 1992, his paper, “Numerical Shape from Shading
and Occluding Boundaries,” was selected as one of the most influential
papers to have appeared in the Artificial Intelligence Journal within thepast 10 years. His IEEE activities include general chair, IROS95,
ITSC00, IV01; program chair, CVPR96, ICCV03; Associate Editor, IEEE
TRA, IEEE TPAMI; distinguished lecture SPS (2000-2002) , RAS (20042006). Dr. Ikeuchi was elected as an IEEE fellow in 1998. He is the EICof the International Journal of Computer Vision.
Masao Sakauchi received the BSci degree inelectrical engineering from the University ofTokyo in 1969 and the MS and PhD degrees inelectronics engineering from the University ofTokyo in 1971 and 1975, respectively. He is a
deputy director general of the National Instituteof Informatics, Japan, and a professor in theInstitute of Industrial Science at the University ofTokyo. He has acted as the general chairman
and program chairman of 10 international conferences and workshops, including the IEEE International Workshop on
Machine Vision and Machine Intelligence (1987), IAPR and IEEE
International Conference on Document Analysis and Recognition(ICDAR ’93) (1993), IEEE International Conference of Multimedia
Processing and Systems (IEEE Multimedia 96) (1996), and Intelligent
Transportation Systems Conference ITSC ’99 (1999). He has authoredmore than 340 refereed papers in the research fields of multimedia
databases, multimedia systems, image processing and understanding,spatial data structures, geographical information systems, and fault
tolerant computing.

. For more information on this or any other computing topic,
please visit our Digital Library at www.computer.org/publications/dlib.