∗†A BSTRACT

targets [20] or laser pointers [18] that can be shown to the camera orcameras. However, these approaches require adjusting the shutter
speed or camera aperture, which prevents simultaneous geometricand photometric calibration.
Our approach extends earlier methods [17, 21], whose implementation is now available in OpenCV [15]. They rely on a printed
grid that is shown to the camera. However, the grid detection isless reliable than using markers [5]. By contrast, our method ismore ﬂexible since it works with any planar textured surface. Furthermore, unlike ours, these methods neither handle extrinsic calibration of non-overlapping ﬁelds of view nor any photometric information. A method that calibrates multiple camera is presented
in [19]. It also recovers the relative poses of multiple cameras but
requires each camera to see all calibration target positions, which,in practice, is a very severe limitation.
In our system, the user only has to move the calibration pattern infront of the cameras. We use a powerful computer vision technique
to detect it in the images [12] both robustly and in real-time. This
provides homographies between the pattern and the images, which
are then used to compute the intrinsic camera parameters as in [21,17]. However, unlike these earlier approaches, we also use themto compute the poses of the cameras with respect to each other.
Finally, we rely on the intensity variations within the pattern andacross the images to achieve photometric calibration. Our all-inone approach includes the following contributions that both reduce
the user’s workload and give our system its additional capabilities:

We propose a fully automated approach to calibrating multiple cameras whose ﬁelds of view may not all overlap. Our technique only
requires waving an arbitrary textured planar pattern in front of thecameras, which is the only manual intervention that is required. The
pattern is then automatically detected in the frames where it is visible and used to simultaneously recover geometric and photometriccamera calibration parameters.
In other words, even a novice user can use our system to extract
all the information required to add virtual 3D objects into the scene
and light them convincingly. This makes it ideal for Augmented
Reality applications and we distribute the code under a GPL license.
1

I NTRODUCTION

We propose a simple-to-use camera calibration system that can handle multiple cameras whose ﬁelds of view do not necessarily overlap. It estimates the geometry of the cameras, their photometricresponses, and an environmental lighting map. The only manualintervention required involves waving an arbitrarily textured planarpattern in front of the cameras. In other words, in one single operation, our system yields all the information required by sophisticatedAugmented Reality applications to draw virtual 3–D objects at theright locations and then light them convincingly, as shown in Fig. 1.
In the remainder of the paper we will refer to the recovery ofgeometric properties as geometric calibration and the estimationof photometric responses and lighting parameters as photometriccalibration. Since geometric camera calibration is now very wellunderstood from a theoretical point of view [9], there are of course
a number of practical techniques, including some commercial systems that are powerful but expensive. Most rely on retro-reﬂective• Fully automated detection of the calibration pattern. Our
approach to detecting the calibration pattern is both real-timeand robust. Furthermore, it automatically selects the framesthat yield the most informative homographies and rejects ambiguous ones. As a result, it becomes easy to collect large
amounts of calibration data, which in turn yields excellent numerical stability and accuracy without any extra manual intervention.

∗ This work was supported in part by the Swiss National Science Foundation† Contact author. E-mail:julien.pilet@epﬂ.ch1-4244-0651-X/06/$20.00 ©2006 IEEE• Handling one or more cameras with non-overlapping
fields of view. Our system selects the position of the patternthat is seen by the largest number of cameras to deﬁne thecommon referential. As depicted by Fig. 2, it then expresses

69

p0

scene. However, they are not well adapted to computing the relative
positions of multiple cameras. Furthermore, they lack robustness.
We therefore focus here on methods that rely on a calibration
object or pattern because they are much more reliable. Commercially available systems rely on tri-dimensional calibration objects
with retro-reﬂective markers, spheres, or disks on them, which areoften complex to build and cumbersome. A notable exception areapproaches such as [13] that only involve a moving wand with twomarkers on it. However, all methods that depend on retro-reﬂective
markers require changing the shutter speed to reliably extract themarkers from the images, which prevents simultaneous geometricand photometric calibration. Similarly, using light patterns emitted
by a laser pointer as in [18] is a very practical approach to geometriccalibration but requires modifying the cameras settings.
Our system falls into the category of recent approaches that rely
on a planar target moved in front of the camera [17, 21]. This is
attractive because such a target can be built by simply printing apattern on a sheet of paper. These earlier techniques, however, were
only designed to recover the internal parameters of a single camera.
Here we are also interested in the positions and orientations of thecameras with respect to each other. The method closest to ours that
we are aware of is presented in [19]. As ours, it provides both the
intrinsic parameters and relative poses of a multi-camera systemusing a planar target. However the parameters are estimated via
factorization of a matrix built from homographies between image
pairs. This requires that all the positions of the planar target must be
seen from all the cameras simultaneously. This is a major limitationthat our method does not have.
The papers discussed above focus on the geometric aspects ofthe problem, and underplay the pattern detection and the frame selection issues. In practice, the target position is often provided byhand to perform the experiments, which is inconvenient. [15] uses
a black and white checkerboard as a calibration object to facilitate
the point extraction procedure. However the repetitive nature of thepattern can easily result in ambiguous poses, and its detection is
known to be error-prone. Here we solve this problem using a stateof the art feature point-based method [11], which can deal with arbitrary textures and is reliable enough not to require manual ﬁlteringof misdetections.

c3

p1

c0

c2
c1
c3

c0 c1

p1

c0
c2

c1

c3

c2

p0

the pose of the cameras that do not see this position by composing relative displacements between camera pairs. Theseare easy to estimate from the homographies and give initial
estimates that are then reﬁned by bundle adjustment.
• Photometric calibration for free. We use the very same
set of images both for geometric and photometric calibration,
which means that no additional user intervention is required toobtain the latter. We take advantage of the fact that our patterneffectively samples the space of surface normals to sample anirradiance map which we deconvolve to render virtual objects
with correct specular reﬂections and cast shadows, such as
those of Fig. 1(c) and 1(d).

2.2 Photometric Calibration

Our system is embedded in a distributed framework that allows
real-time processing of the output of several cameras plugged into
separate but networked computers. To demonstrate its effectiveness, we distribute the code under a GPL license [1].
In the remainder of the paper, we ﬁrst compare our approach toexisting calibration techniques. We then introduce the techniques
we use to recover ﬁrst the geometric camera pose parameters andthen the photometric and illumination ones. Finally, we evaluate
the accuracy of our system and present Augmented Reality results.
2

Real-world illumination can be captured using calibration objects
such as reﬂective spheres. This can be done as a preprocessing
step [3] if the illumination remains unchanged. It can also be donein real-time [10] using a calibration object built by adding a mirror
ball to a 2–D square marker. While effective, this approach is much
more difﬁcult to deploy than ours since constructing such an object is much more involved than simply printing a textured pattern.
Similarly, [16] relies on omni-directional stereo cameras, which requires specialized hardware instead of the ordinary cameras we use.
The photometric calibration we get may not be as accurate as thoseobtained with such lighting probes but we will show that it is amply sufﬁcient to synthesize realistic augmented images while being
much more light-weight.
Another class of approaches [4, 6] focuses on interactive rendering of illumination changes caused by virtual objects in the realscene. However, this typically involves a 3–D scene model, which
is cumbersome to acquire.
Our approach to creating an environmental lighting map is related to the inverse lighting framework of [14] that involves an object of constant albedo and known shape. If the object surface contains a sufﬁciently large number of differently oriented normals, it
is possible to relate a lighting contribution to each direction of adiscretized lighting sphere. Regularized deconvolution then allows
the estimation of the light sources position. In our case, as the grid

R ELATED W ORK

We review brieﬂy existing approaches to geometric and photometric calibration, which are usually performed separately. Our approach is inspired by several of them but performs both simultaneously, which makes it easier to use for Augmented Reality purposes.
2.1 Geometric Calibration
Calibration algorithms that do not require an object known a priori are sometimes called auto-calibration algorithms. They usually
involve a moving camera that is calibrated by simultaneously recovering pose parameters and reconstructing the 3–D structure of the70

that maps each point 1 [X,Y ] of the calibration pattern plane to its
corresponding 2D point [u, v] in camera c for a given pose p. The
scalar λ accounts for the fact that homographic transformation isexpressed in homogeneous coordinates.

moves in front of the cameras, it samples the spaces of normals.
Since we precisely compute those normals, we can directly exploit
the observed changes in pixel intensities to model the illumination
using a very similar approach. This removes the need for a constantalbedo 3–D object and its model.
3

G EOMETRIC C ALIBRATION

H2

H1

Our geometric calibration procedure includes several stages. Our
system ﬁrst computes homographies between the geometric planeof the target object and its image projections. It then retains themost reliable ones and the corresponding frames to estimate the intrinsic camera parameters and the relative pose of the calibration
object with respect to the cameras in the corresponding frames. Inturn, these poses are used to select a common referential and tocompute the positions and orientations of the cameras with respect
to each other. Finally, our system performs global non-linear minimization to reﬁne these estimates. We outline the individual stepsof this process below.

H4

H5

H3

Fronto−parallel viewof the calibration
objectH3
H5

3.1 Computing and Selecting the best Homographies
To compute the homographies between the geometric plane of thetarget object and its image projections, we must ﬁrst detect it in as
many frames as possible. To this end, we use an interest point-based
method that we developed in earlier work [11]. It ﬁrst performs
real-time matching of interest points extracted from a reference image such as the ones of Fig. 3 against those extracted from input
images acquired at run-time under potentially large perspective andscale variations. It then uses standard robust estimation techniques
to compute homographies from these matches and selects the best.

3.1.3 Selecting the Best Homographies
Unlike checkerboard-based methods, our approach to estimating
the homographies is sufﬁciently robust not to generate erroneous
matches that would have to be removed by hand. What happens inpractice is that when the pattern is either occluded or too slanted, it
is simply not detected.
However, this is still not quite enough because some of these
homographies are inherently ambiguous or singular from a calibration point of view. As shown in Fig. 4, these unreliable homographies come in two ﬂavors. First, the ones that insufﬁciently distortthe pattern are ambiguous and, depending on the noise, may yieldwildly different pose estimates. Second, those that distort the pattern too much also produce unreliable pose estimates because theinterest points become difﬁcult to locate precisely enough.
To overcome this problem, we deﬁne a square in the plane attached to the pattern and measure the angle at each corner after
warping it by the homography. If one of the angles is too large or
too close to π2 , the homography is rejected. More formally, a homography H is rejected if at least one of the angles α of the warped
square veriﬁes

3.1.1 Wide-Baseline Matching
We achieve fast and robust wide-baseline matching by formulating
it as a classiﬁcation problem, which lets us shift much of the computational burden to a training phase. We treat the set of all possible appearances of each individual interest point as a class, which
we refer to as its view-set. During training, given a fronto-parallelview of the target object, interest points are extracted and numerous
synthetic views of their possible appearance under perspective distortion are generated and used to train a set of randomized trees [2].
These trees are then used at run-time to recognize the distorted interest points by deciding to which view-set, if any, their appearance
belongs. In theory, we could have used other kinds of classiﬁers forour purpose but we chose randomized trees because they are robust
and fast, while naturally handling multi-class problems and being
reasonably easy to train.
In this context, a simple and fast interest point detector such as
the Harris corner detector [8] sufﬁces for operation even under large
perspective and scale variations, which results in a further performance increase. While earlier methods require detectors that canbe depended upon to produce very repeatable results, which can bevery time-consuming, the most repeatable object keypoints for the
speciﬁc target object are simply found during the training phase.

cos(α ) > cos(α0 ) orcos(α ) < cos(π − α0 ) orcos( π2 + α1 ) < cos(α ) < cos( π2 − α1 ).
Good results have been achieved using α0 = 0.01 and α1 = 0.005.
The ﬁltered homographies are then sorted by the likelihood assigned by the robust estimator of Section 3.1.2 and only the best
ones are kept. In practice, we retain around ﬁfty for each camera.
As will be shown in Section 5.1, this is enough to guarantee accuracy without imposing an unnecessary computational burden.

3.1.2 Estimating the Homographies
Given enough correspondences, the homography between the target image and the input one is easily estimated from the correspondences by RANSAC followed by non-linear least-squares estimation. For each successful such computation, this gives us a homography Hc←p such that:
⎤
⎤
⎡
⎡
u
X
(1)
λ ⎣ v ⎦ = Hc←p ⎣ Y ⎦11

3.2 Initial Estimation of the Internal Parameters
The internal parameters are ﬁrst estimated for each camera individually from the homographies Hc←p using a method similar tothe ones of [17, 21]. The computation is quite standard and is described in appendix for the sake of completeness. This yields for
1 The symbol  denotes the transpose operator throughout the paper.

71

(a)

each camera c a matrix of internal parameters⎤⎡
τc fc 0 u0c
f c v0c ⎦ ,Kc = ⎣ 0
0
0
1

(b)

(c)

p0

p1
[R,t] p1

[R,t] c(2)2where f c stands for the focal length, τc for the aspect ratio, and(u0c , v0c ) for the principal point. These internal parameters of
all the cameras will then be reﬁned together with the external ones
during the non-linear global minimization of Section 3.4.

c

2

p

0

[R,t] c3c0

3.3 Initial Estimation of the Poses

c1

1

c2

p0

p1

We recover the external parameters of each camera in a common
referential in two steps. First, our algorithm computes external parameters in a coordinate system attached to the calibration pattern
for each frame independently by making use of the internal parameters as estimated previously and the homography related to theframe. It then selects a common referential and computes camera
poses in this referential by composing rotations and translations between pairs of frames.

pc3[R,t]c3

c3[R,t]c3 = [R,t]c3 ←p1 [R,t] p1 ←c2 [R,t]c2 ←p0

the pose of p0 with respect to camera c2 with the displacement between cameras c2 and c3 , which can itself be estimated from theposes of p1 with respect to these two cameras.
Recovering such chains amounts to compute paths between
nodes of a connected graph, which is well understood from an algorithmic viewpoint [7]. More speciﬁcally, we deﬁne a graph whose
nodes correspond to the cameras and the poses of the calibration
object, such as the one depicted by Fig. 5. An edge links a camera
node and a pose node when the camera sees the pose and is labeledwith the corresponding displacement.
Let p0 denote the pose that deﬁnes the common referential, and[R, t]c the pose parameters of camera c in this common referential.
Since we have by deﬁnition [R, t]c = [R, t]c←p0 , we ﬁrst look for apath between the p0 and c nodes, which we write as

3.3.1 Displacements between the Calibration Object and the Individual Cameras
Given a homography and the intrinsic parameters, we estimate thedisplacement of the calibration object with respect to the camera as
described in the appendix. This step gives us the relative rotation R
and translation t of the calibration object with respect to the camera.
Together, they represent the rigid displacement corresponding topose p as seen from camera c. We write this displacement as the
4 × 4 matrixRc←p tc←p[R, t]c←p =
.
(3)
0
1

p0 → cσ (1) → pσ (2) → cσ (3) → . . . pσ (n) → c .
where σ (.) is a mapping function on the indices that deﬁnes thepath. Note that the path alternates object pose nodes and camera
nodes. It gives us a way to compute [R, t]c from the displacements
we just computed since

The reverse displacement is computed by inverting the matrix
which we denote as [R, t] p←c .

[R, t]c=[R, t]c←p0 =[R, t]c←p3.3.2 Handling Non-Overlapping Cameras

σ (n). . . [R, t]c

σ (3) ←pσ (2)[R, t] p

σ (2) ←cσ (1)[R, t]cσ (1) ←p0.

The [R, t] p pose parameters of the calibration object can be estimated similarly.
Obviously, this will only work if the graph has one single connected component. In practice, assuming that no camera has a ﬁeldof view that does not overlap at all any of the others, this is always
the case if we move the calibration pattern sufﬁciently.

In practice, the calibration object may never be seen by all cameras
simultaneously. Our system therefore selects as a common referential the one attached to the pose of the calibration object seen by thelargest number of cameras. It then expresses all the external camera parameters in this referential by composing the displacementsof Eq. 3 and their inverses.
Fig. 2 illustrates this behavior. In this case, p0 provides the common referential because it is seen by three cameras whereas p1 isseen by only two. Even though camera c3 does not see p0 , its external parameters in this referential can be estimated by composing

3.4 Refining the Estimation
Computing displacements by composing pairwise motions is effective but not particularly accurate. Therefore, to reﬁne not only the72

only on its normal. In other words, each image in which the pattern is detected provides a number of samples corresponding to one
individual surface orientation. Because we can easily and automatically collect many such samples, we can simultaneously recover thegain and bias of each camera and create an environmental lighting
model that can be used for relighting purposes.
In this section, we ﬁrst introduce the illumination model we use
for complex environments such as the one depicted by Fig. 1(a). We
then present two complementary approaches to instantiating it. The
ﬁrst one involves solving a linear system of equations derived fromthe calibration sequence to express the illumination as a functionof surface normals. Once the gains and biases have been estimated,
this model can be incrementally updated to reﬂect lighting changes.
The algorithm can therefore be embedded into a real-time application that handles non-constant lighting. However, it is not designed
to synthesize either shadows or specular reﬂections. We have therefore developed a second approach based on deconvolution. It explicitly computes a light distribution that could have produced theobserved pixel intensities. It is more computationally intensive andmakes no provisions for time-varying lighting but allows added virtual objects to cast shadows and produce realistic specularities.

pose parameters but also the intrinsic ones, we minimize with respect to all cameras simultaneously the sum of the reprojection errors for the point correspondences used during the detection step ofSection 3.1.2. This bundle adjustment is expressed asC
P M(c,p) ∑ ∑ ∑ c=1 p=1 k=1 uc,p,kvc,p,k



⎞⎞2
⎛
⎛

Xc,p,k

⎠⎠ ,
⎝
,
Y
− proj⎝Kc , [R, t]c [R, t]−1
c,p,k
p


0

(4)

where
• C is the number of cameras, P the number of poses of thecalibration object, and M(c, p) the number of matches foundby the detection stage for camera c and pose p. M(c, p) = 0 if
pose p of the calibration object is not seen by camera c;• [uc,p,k , vc,p,k ] and [Xc,p,k ,Yc,p,k , 0] are respectively are a 2–
D point and a 3–D point matched by the detection step;• proj(K, [R, t] , M) returns the projection of 3–D point M underpose [R, t] and internal parameters K.
To increase the robustness of our algorithm, we introduce a simple
robust estimator in Eq. 4 to eliminate potentially incorrect correspondences.

4.1 Illumination Model
In a natural environment with extended light sources such as a room
lit by rectangular windows, we cannot assume that we only deal
with a small number of point light sources. Instead, we consider
a potentially inﬁnite number of sources that we assume to be both
directional and outside the capture volume. For our purposes, thisprovides a satisfactory approximation of extended light sources.
By printing the calibration pattern on matte paper, we ensure thatit reﬂects light equally in all directions. As a result, the difference
in intensity values between images taken at the same time by two
different cameras are due to shutter speed or aperture, but not tocamera pose. If we further assume that the same sources are visible from every point of the calibrated volume, the amount of lightreﬂected by a point on the calibration pattern depends only on thesurface normal nt at time t. To be robust to small localization errors,
we do not consider individual points, but small patches π that average the local property around points on the calibration pattern. The
irradiance x at surface patch π and time t can therefore be writtenas

3.5 Online vs. Offline Processing
All the computations up to the initial estimation of pose parameters of Section 3.3 are performed in real-time on an ordinary PC.
The reﬁnement stage of Section 3.4 is much slower and requires a
few minutes for long calibration sequences. In practice, this meansthat we can run the system in real-time mode until the graph of
Section 3.3.2 has indeed one single connected component and we
can obtain initial estimates, which will then be reﬁned ofﬂine. This
makes the system extremely easy to use.
In real-time mode, the processing time varies from one frame tothe other. It is therefore important to drop frames in a synchronous
manner: If a camera drops a frame, all other cameras should also
drop it. Since in our implementation the cameras are connected
to several computers, we developed a networked solution to thisproblem: A process is attached to each camera and is considered
as a network client that connects to a server. Every time a newframe arrives, the client informs the server and puts the image ina queue. The server either ask the client to accept or to drop theframe, depending on the other cameras and on the current computational load. If a frame is accepted, every client is informed andputs the image in a waiting-for-computation queue. At this point,
the computation thread of each client, if ready, can start to detect
the calibration pattern in the image. Upon termination, it sends theresult to the server. Once the server has collected homographiesand matches from every client, it can either accumulate data forcalibration or, if the calibration has already been done, augment theimages.
4Lxπ ,t = ∑ max (nt dl , 0) Ωl ,

(5)

l=0

where Ωl represents the radiosity, or power, of source l and dl its
direction. We assume a linear response for the cameras and write
the pixel intensity I(c, π ,t) of patch π in the image acquired bycamera c at time t asI(c, π ,t) = gc aπ xt + bc ,(6)

where aπ is the average surface albedo over π , gc the camera gain,
and bc its bias. In practice these quantities can only be known up toa scale factor. We use for aπ the mean intensity over π in an image
acquired under uniform diffuse lighting such as the ones of Fig. 3,
which simply amounts to selecting a particular scale factor.

P HOTOMETRIC C ALIBRATION

Convincingly adding virtual objects into a scene not only requires
proper registration, which is what the technique of Section 3 provides, but also photometric modeling so that they can be relighted
correctly and properly blend into the environment.
As discussed in Section 2, most existing techniques perform geometric and photometric calibration as independent steps that require different camera settings. By contrast, our approach relies on
the very same set of images to perform both kinds of calibration.
For each camera, as long as the lighting and camera settings do notchange, the pixel intensities within the calibration pattern depend

4.2 Online Lighting Calibration
Instead of explicitly evaluating the Ω radiances of Eq. 5, we directlycompute the xπ ,t irradiances as a function of the orientation of the nnormals. More speciﬁcally, we simultaneously compute the gains,biases, and xπ ,t irradiance values for the normals we have observed.
We then interpolate this set to estimate the unobserved values.
This yields a light map M(n) such as the one depicted by
Fig. 6(a). When augmenting an image c with a virtual surface734.3 Offline Estimation of Light Distribution for Rendering
Specular Effects and Casting Shadows(a)While the previous method yields satisfying results for fast online
pre-visualization, its accuracy suffers from two numerical problems. First, it does not minimize a physical error since the problem
is intrinsincally non-linear. Second, it solves simultaneously forcamera-related parameters, the gains and the biases, and lightingrelated parameters, the irradiances. As is the case for geometriccamera calibration, performing non-linear optimization and dissociating the estimation of the internal parameters, here the gains andbiases, from that of the external ones, here the irradiances, yields
more reliable results.
Therefore, we developed an ofﬂine approach to doing so. It is
computationally more expensive than the online technique of Section 4.2 but produces a more sophisticated illumination model thatallows for specular highlights, cast shadows, and changing the materials of the virtual objects.
We model the lighting environment as a regularly sampled domeof lights of varying power, such as the one depicted by Fig. 6(b).
We begin by estimating the gain and bias of each camera in a waythat is independent from lighting effects and normalizing the pixel
intensities. We then apply a regularized deconvolution algorithm on
the observed pixel intensities within the calibration object to assign
to each individual light the power that best explains the observations.

(b)whose normal is n, the augmented pixel value is set to I(v, c,t) =gc av Mt (nv ) + bc , where av is the virtual albedo. This lightmap rendering is easily done by the GPU and requires only a short OpenGL
shading language program.
4.2.1 Linear Estimation of the Light Map
The geometric calibration process provides many surface normals
nt and pixel values I(c, π ,t). To solve for the unknown gains, biases, and radiances, we linearize the problem by replacing somevariables in Eq. 6. Let
gc

=

bc

=

4.3.1 Relative Cameras Responses
To estimate the gains and biases, we consider pairs of views taken
at the same time by two different cameras. As we will see, thisremoves the need from separately computing the surface albedo aπ
and the irradiance xπ , thus avoiding an additional source of noise.
Let D be the set of triplets (c, d,t) such that both cameras c and d see
the calibration object at time t. Assuming a Lambertian behavior of
the calibration pattern, we can compute our gains and biases as

1
,
gc
bc
.
gc

For each patch in each detected frame, Eq. 6 can be rewritten as⎡⎤gc−I(c, π ,t)gc + aπ xt + bc = [−I(c, π ,t) 1 aπ ] ⎣ bc ⎦ = 0 . (7)xtargming1 ...gC ,b1 ...bC (c,d,t)∈D π

Putting all these equations together yields a large but sparse linear
system and we ﬁnd our solution as the eigenvector associated to thesmallest eigenvalue.

f (n) =

exp
0

. (8)

4.3.2 Deconvolution
From Eqs. 5,6, and 9, the normalized intensity at time t of a patch
π with normal n can be written asL¯ π ,t) = aπ ∑ max(nt · dl , 0)Ωl ,
I(c,

(10)

l=0

where L is the number of sources used to represent the lighting environment and dl the direction of source l with power Ωl . Our objective is to recover the powers Ωl of each source so as to minimize∀n : Mt+1 (n) = (1 − f (n)) Mt (n) + f (n)xt+1 ,− cos−1 (n · nt+1 ) 2s12

2and express the deconvolution in terms of these normalized images.

The above computation assumes that the lighting does not change
during calibration. However, once the gains and biases are known,
computing a new surface intensity xu from a new observed pixel
cπ ,u is trivial. Thus, if the illumination changes, new frames canupdate the light map at the same time it shades a virtual scene.
Each frame can sample only one normal at a time. Thus, if the lightchanges suddenly, it is not possible to update the whole irradiance
map at once. Instead, we locally update the irradiance around themeasured normal and keep the old values for other normals.
Let Mt (n) be the light map at time t To update its value for asurface of normal n with the recently computed sample xt+1 corresponding to the observed real surface orientation nt+1 , we writeI(c, π ,t) − bc I(d, π ,t) − bd−
gc
gd

We can now deﬁne a normalized image I¯ from the image taken bycamera c as
¯ π ,t) = I(c, π ,t) − bc ,
(9)
I(c,
gc

4.2.2 Incrementally Updating the Light Map

where

∑ ∑

εobs (Ωl ) = ∑if cos (n · nt+1 ) > 0
otherwisec,π ,tand s a constant blurring factor. In practice, our system initializes
the light map using the solution of the linear system of equations
and then updates it for each new frame.

L¯ π ,t) − aπ ∑ max(nt · dl , 0)Ωl
I(c,

2
,

(11)

l=0

which represents the sum of the squares of the differences between
the intensities predicted by the lighting model and the observed
ones. Unfortunately, since the Lambertian model acts as a low-pass

74

Focal stability (5 cameras, 1 sequence)

1900

Camera 1
Camera 2
Camera 3
Camera 4
Camera 5

1850

40

60

80

100
120
140
Number of homographies

160

250
200
150
100
50
0

180

200

1900

Focal length

Optimization time [seconds]

Focal length

1950

300

1950

1800

Focal stability

Optimization time evolution

2000

1850

1800

1750
Camera 2, sequence A
Camera 2, sequence B

20

40

60

80

100 120 140 160 180 200

1700

Number of homographies

100
120
140
Number of homographies

(b)

(c)

(a)

ﬁlter on the lighting, this is an ill-posed deconvolution problem thathas many potential solutions. To obtain a realistic one, we must ﬁrst
constrain the Ωl to be positive as not doing so would allow the existence of lights with negative power.
√ To this end, we optimize withrespect to the square roots Ωl = Ωl , instead of the Ωl themselves.
We then deﬁne two regularization terms to be added to εobs (Ωl ).
The ﬁrst one is a smoothness term similar to the one introduced in[14] that we write as


2
2 2
Ωi − Ωj
εad j Ωl = ∑
.
(12)

τf
f
u0
v0

Cam 1
-0.16%-0.08%
5.57%
0.11%

40

60

Cam 2
0.32%
0.33%
6.49%
0.76%

80

Cam 3
0.27%
0.19%
9.81%
3.15%

Cam 4
1.12%
1.22%
5.85%
-11.7%

160

180

200

Cam 5
0.56%
0.53%
6.00%
-3.56%
Kc

i, j adjacent

Minimizing it forces adjacent sources to be of relatively similar
powers. The second is designed to encourage the good localizationof narrow light sources and is written as
2
2
L 


 2
εstick Ωl = ∑ Ω0,l − Ωl.
(13)photometric ones returned by the procedure of Section 4 lets us
augment the scene with objects that are correctly registered andlighted.

l

In practice, we therefore solve
argmin εobs (Ωl ) + αεad j (Ωl ) + β εstick (Ωl ) ,Ωl

5.1 Accuracy of the Geometric Calibration Parameters(14)

To test the accuracy of the camera parameters that our system recovers, we ﬁrst trained a classiﬁer to recognize the interest points ofthe pattern of Fig. 3(a), as discussed in Section 3.1.1. We then useda 5 cameras setup to record two different sets of video sequences bywaving the pattern in front of the cameras. Finally, we performed
the calibration independently for each set.
Fig. 7(a) depicts the focal lengths recovered using the ﬁrst setof sequences. They are shown as a function of the total numberof homographies retained to perform the computation. When thisnumber climbs above 140, the estimates become quite stable. It
therefore does not make sense to use many more since the computational cost increases almost linearly with this number, as shownin Fig. 7(b).
In Fig. 7(c) and Table 1 we compare the results obtained independently using the two sets of video sequences. In Fig. 7(c), we
superpose the focal length estimates, again drawn as a function of
the total number of homographies retained. As before, once we use
more than 140, the two estimates become very close. This is a very
good indication that they are accurate since they were computed independently. As shown by Table 1, this is true not only for the focal
lengths but also for the principal point locations.
We have not obtained ground truth for the external camera parameters. However, as evidenced by the supplementary videos, thevirtual object appears to be very stable with respect to the calibration pattern, which would not be the case if they were poorly recovered.

where α and β are regularization constants. To this end, we use
a standard Levenberg-Marquardt least-squares solver but deﬁne a
speciﬁc optimization schedule that initially favors a relatively small
subset of light sources and then progressively allows more sources
to contribute to the lighting. At the ﬁrst iteration, the initial powerof the sources Ω0,l are set to zero and the parameter β is given ahigh value. This results in an initially small set of sources withsigniﬁcant amounts of power while the majority of the remaining
ones are still turned off. The next iterations start from the previous
power estimation and the β parameter is progressively decreased toallow the recovery of more wide-spread sources. This results is alighting sphere with smooth but well-deﬁned clusters instead of acompletely smoothed out distribution.
The knowledge of the source positions allows for realisticallyshaded virtual objects. Realistic cast shadows and specularities canbe computed using a Phong shading model. As shown in the next
section, the conjunction of tracking and lighting distribution estimation allows us to render convincing augmented images.
5

R ESULTS

In this section we ﬁrst evaluate the accuracy and reliability of thegeometric camera calibration of Section 3. We then show that using both the corresponding geometric camera parameters and the75

5.2 Adding and Relighting Virtual Objects

We associate to each camera a projection matrix
P = K [R | t] ,We used a three camera set-up similar to the one depicted by
Fig. 1(a) to produce the augmented images of Fig. 8. In the ﬁrst
row, we use only the geometric calibration parameters to draw thevirtual teapot at the right place and use a randomly selected point
light source to relight it. Even though the teapot is correctly registered, the result is unconvincing because the shading patterns do
not match those of the other real objects present in the scene. In the
middle row, we use the output of the online photometric calibration
procedure of Section 4.2 to relight the object. The result is much
improved but highlights and shadows are still missing. As shownin the bottom row of the ﬁgure, using the output of the more sophisticated ofﬂine calibration procedure of Section 4.3 solves bothproblems. Both highlights and shadows now appear at the rightplaces, thus signiﬁcantly increasing the realism. We supply the corresponding videos as supplementary material.
Fig. 9 showcases the ﬂexibility that our multi-camera systemprovides. In two of the three images, the calibration pattern is not
visible in the view we are trying to augment and a monocular approach that relies on detecting it would fail. However, because it
is seen by another camera and because the relative positions of thecamera with respect to each other have been computed, we can nevertheless draw it at the right location, as evidenced by the fact thatthe real box occludes it correctly.
Our system can use any sufﬁciently textured pattern for calibration purposes. In Fig. 10, we demonstrate this by using the pieceof cardboard of Fig. 3(b) to calibrate and augment the scene. Notethat the virtual teapot casts a realistic shadow on the real cardboard,
thereby indicating that both the geometric and photometric calibration are correct. We have also successfully used the book of
Fig. 3(c) for this purpose.
6(15)

where R is a 3×3 rotation matrix and t a translation vector. Without
loss of generality, we can choose Z = 0 as the plane of our calibration pattern. The relation between K, R, and t and the homography
H that maps a point M = [X,Y, 0] of this plane to its corresponding
2D point m = [u, v] under perspective projection can be written as
⎤
⎡
⎤
X
X
Y
⎥
⎢
= K [r1 r2 t] ⎣ Y ⎦ ,∝P= K [r1 r2 r3 t] ⎣
0 ⎦
1
1
(16)
where r1 , r2 and r3 respectively are the ﬁrst, second and third column of the rotation matrix R, and the symbol ∝ means “equal up
to a scale factor”. Identifying the terms of Eq. 1 and Eq. 15 yields
H ∝ K [r1 r2 t]. For convenience, we introduce a 3 × 3 matrix Twith⎡⎤
1 0 t 1
T=⎣ 0 1 t2 ⎦,0 0 t 3




m
1



M
1



⎡

where t = R−1 t = R t that lets us write
H ∝ KRT .

(17)

To ﬁnd the relations between the coefﬁcients of H and the internal parameters, let us now compute the product H ω H where ωis the matrix (KK )−1 , also known as the image of the absolute
conic. We haveH ω H = H K− K−1 H ∝ (KRT) K− K−1 (KRT)C ONCLUSION
andWe have presented a system for simultaneous geometric and photometric calibration of multiple cameras that is extremely easy touse. At run-time, the only manual intervention required involves
waving a planar calibration pattern in front of the cameras. Furthermore, it imposes very few constraints on the relative positions
and orientations of the cameras and can handle complex real-world
illuminations.
A natural extension would be to correct for distorion [21] and vignetting. However, even without it, our system yields both accurate
camera parameters that let us add virtual 3D objects at the appropriate locations and a lighting model that allows us to relight themconvincingly. Since these are the two ingredients that are requiredby Augmented Reality applications, we believe that our approachto calibration is ideally suited for this purpose. We urge interested
readers to download the source code to try it for themselves [1].

(KRT) K− K−1 (KRT)= T R K K−⎤K−1 KRT =
⎡
1
0
−t 1


⎣
0
1
−t
T T=
2 ⎦.



−t 1 −t 2 t 2
By considering the expressions of the elements of the 2 × 2 submatrix on the top-left of the T T matrix, we obtain the two following equations: 
h1 ω h1 − h
2 ω h2 = 0 ,
(18)
ω
h
h
2=0
1
where hi represents the column i of H. It follows that:
⎡
⎤
1
0
−u0
⎦.
τ2
−τ 2 v0ω ∝⎣ 0−u0 −τ 2 v0 τ 2 f 2 + u20 + τ 2 v20

A PPENDIX(19)

By rearranging the terms of Eq. 18 and Eq. 19, we obtain the following linear system in some of the coefﬁcients of ω :Internal Parameters from a Set of Homographies
The computation of internal parameters from a set of homographies
is quite standard [17, 21], but we give it here for the sake of completeness. First, we write the matrix of internal parameters as⎡⎤
τ f 0 u0
f v0 ⎦ ,
K=⎣ 0
0 0 1

AW = h

(20)

where


where f represents the focal length, τ the aspect ratio, and (u0 , v0 )
the principal point. We want to estimate K from a set of homographies that map points on a planar object to points on captured
images for different object positions.

76

A=

2(h11 h31 − h12 h32 ) h221 − h222 2(h21 h31 − h22 h32 ) h231 − h232h11 h32 + h12 h31
h22 h21h32 h21 + h22 h31
h32 h31
⎤
⎡
ω13


⎢ ω ⎥
h211 − h212.
W = ⎣ 22 ⎦ and h =ω23h11 h12ω33

First row
Second row
Third row

Each homography yields such a pair of equations. For each camera,
this produces an over-constrained system that we solve in the leastsquares sense. The internal parameters u0 , v0 , f , and τ can thenbe estimated from ω13 , ω22 , ω23 , and ω33 . In our implementation,
they are reﬁned by the ﬁnal non-linear optimization.

follows.
Recall from Eq. 17 that H ∝ KRT, where R and t are expressedin a coordinate system attached to the planar object. Equivalently,
we write:RT ∝ K−1 H .

Displacement between a Camera and a Planar Object from a
Set of Homographies and the Internal Parameters

Since the columns of R should have a norm equal to 1, the scale
factor can be retrieved, and a ﬁrst estimation of these columns is
obtained as:−1

K h1
r1 = K−1 h  ,Once the internal camera parameters have been recovered from thewhole sequencem we can recover the rotation R and the translationt between a particular pose of the planar object and the camera as

1

−1

K h2
r2 = K−1 h  ,2

r3 = r1 × r2 ,where × denotes the cross-product of two vectors. Because the ho-77

mographies are noisy, the resulting rotation matrix is not orthonormal and we correct it using the procedure given in the appendix
of [21] which seeks the closest orthonormal matrix in the Frobenius sense using a Singular Value Decomposition. Similarly thetranslation vector t can be approximated as:t= [10] M. Kanbara and N. Yokoya. Real-Time Estimation of Light Source
Environment for Photorealistic Augmented Reality. In International
Conference on Pattern Recognition, 2004.
[11] V. Lepetit and P. Fua. Keypoint recognition using randomized trees.
IEEE Transactions on Pattern Analysis and Machine Intelligence,
2006. In press.
[12] V. Lepetit, P. Lagger, and P. Fua. Randomized Trees for Real-Time
Keypoint Recognition. In Conference on Computer Vision and Pattern
Recognition, San Diego, CA, June 2005.
[13] H.G. Maas. Image sequence based automatic multi-camera systemcalibration techniques. ISPRS Journal of Photogrammetry and Remote Sensing, 54(5-6):352–359, December 1999.
[14] S. R. Marschner and D. P. Greenberg. Inverse lighting for photography. In Proceedings of the Fifth Color Imaging Conference, Society
for Imaging Science and Technology, 1997.
[15] Open
SourceComputer
VisionLibrary.
Intel.
http://www.intel.com/technology/computing/opencv/.
[16] I. Sato, Y. Sato, and K. Ikeuchi. Acquiring a Radiance Distribution to
Superimpose Virtual Objects onto a Real Scene. IEEE Transactions
on Visualization and Computer Graphics, 1999.
[17] P. Sturm and S. Maybank. On Plane-Based Camera Calibration: A
General Algorithm, Singularities, Applications. In Conference on
Computer Vision and Pattern Recognition, pages 432–437, June 1999.
[18] T. Svoboda, D. Martinec, and T. Pajdla. A convenient multi-camera
self-calibration for virtual environments. PRESENCE: Teleoperators
and Virtual Environments, 14(4):407–422, August 2005.
[19] T. Ueshiba and F. Tomita. Plane-based calibration algorithm for multicamera systems via factorization of homography matrices. In International Conference on Computer Vision, pages 966–973, 2003.
[20] Viconpeak. http://www.vicon.com/.
[21] Z. Zhang. A ﬂexible new technique for camera calibration. IEEE
Transactions on Pattern Analysis and Machine Intelligence, 22:1330–
1334, 2000.

2K−1 h3
 .
−1
K h1  + K−1 h2 

R EFERENCES[1] CVLab Software.
.
[2] Y. Amit and D. Geman. Shape Quantization and Recognition with
Randomized Trees. Neural Computation, 9(7):1545–1588, 1997.
[3] P. Debevec. Rendering synthetic objects into real scenes: Bridging traditional and image-based graphics with global illumination and highdynamic range photography. In ACM SIGGRAPH, July 1998.
[4] G. Drettakis, L. Robert, and S. Bougnoux. Interactive common illumination for computer augmented reality. In Eurographics Rendering
Workshop, pages 45–56, June 1997.
[5] M. Fiala and C. Shu. Fully Automatic Camera Calibration Using SelfIdentifying Calibration Targets. Technical report, National ResearchCouncil Canada, 2005.
[6] S. Gibson and A. Murta. Interactive Rendering with Real-World Illumination. In Eurographics Workshop on Rendering, June 2000.
[7] J. Gross and J. Yellen. Graph theory and its applications. CRC Press,
Inc., Boca Raton, FL, USA, 1999.
[8] C.G. Harris and M.J. Stephens. A Combined Corner and Edge Detector. In Fourth Alvey Vision Conference, Manchester, 1988.
[9] R. Hartley and A. Zisserman. Multiple View Geometry in Computer
Vision. Cambridge University Press, 2000.

78