Discrete camera calibration from pixel streams
Etienne Grossmann1 , José António Gaspar and Francesco Orabona

Abstract
We consider the problem of estimating the relative orientation of a number of individual photocells -or pixels- that hold fixed relative positions. The photocells measure
the intensity of light traveling on a pencil of lines. We assume that the light-field thus
sampled is changing, e.g. as the result of motion of the sensors and use the obtained
measurements to estimate the orientations of the photocells.
Our approach is based on correlation and information-theory dissimilarity measures. Experiments with real-world data show that the dissimilarity measures arestrongly related to the angular separation between the photocells, and the relation canbe modeled quantitatively. In particular we show that this model allows to estimate
the angular separation from the dissimilarity. Although the resulting estimators are
not very accurate, they maintain their performance throughout different visual environments, suggesting that the model encodes a very general property of our visual world.
Finally, leveraging this method to estimate angles from signal pairs, we show how
distance geometry techniques allow to recover the complete sensor geometry.

1. Introduction
Early experiments in psychological research, such as the famous experiments withthe prism and the inverted glasses [16, 32, 31], revealed that much of the geometry of
the human’s vision is known beforehand and is saved in the brain as imaging experiences. In particular, those psychological experiments showed that a person wearing
distorting glasses for a few days, after a very confusing and disturbing period, could
learn the necessary image correction to restart interacting effectively with the environment.
This amazing learning capability of the human’s vision system clearly contrasts
with current calibration methodologies of artificial vision systems, which are still strongly
grounded to the a priori knowledge of the projection models. Calibration is typicallybootstrapped by using some basic properties of the imaging system. In particular, in all
the cases that we are aware of, the local topology is known: pixel (i, j) is a neighbor of(i + 1, j); the triplet (i, j), (i + 1, j), (i + 2, j) is approximately aligned. It is possible
to detect and localize edges, extrema, corners and features in uncalibrated images, just
1 This work was partially supported by TYZX, Inc, by the Portuguese FCT POS_C program that includes
FEDER funds, and by the EU-project URUS FP6-EU-IST-045 062

Preprint submitted to Elsevier

April 1, 2009

Figure 1: A discrete camera consists of a number of photocells (pixels) that measure the light traveling along
pencil of lines. More generally, one could consider that each photocell is characterized by a point spreadfunction defined on the 3D sphere.

as well as on calibrated images. In short, it is undeniable that locally, and for many
practical purposes, an uncalibrated image looks just like a calibrated image.
This paper departs from traditional computer vision on this point: we do not assume
that pixels (i, j) and (i + 1, j) are neighbors. In fact, we are given pixels indexed bya single index i. Since such a camera is peculiar by the absence of a-priori known
topological information, it is only natural to call it a discontinuous, or discrete camera.
The problem that we address is to determine the geometry of a discrete camera that
produced a stream of values x (i, t), where the integer index i denotes the photocell
or pixel, and t, the time. This problem may seem too general to pertain to computer
vision. Indeed, digital cameras are unlikely to ever produce discontinuous images2 .
Yet we must also consider the computation that occur in biological systems, where it isunclear how much of the geometry of the sensor is known beforehand, and how much isdetermined by processing visual stimuli. We must also consider robots equipped witharrays of photocells, and new visual sensors [28] that may have variable geometry.
Our goal here is thus to determine the 3D direction pointed by each pixel, under
the assumption that pixels sample light that travels along rays that intersect in a uniquecenter of projection (see Fig. 1). In the impossibility of performing local image processing, traditional calibration techniques [34, 12] are out of the question.
Less traditional non-parametric methods that assume a smooth image mapping andsmooth motion [20, 10] can obviously not be applied either. By using controlled-lightstimuli or known calibration, matches could be obtained, allowing to use match-based
non-parametric techniques [25]. In the technical report [6], a non-central projection
sensor samples the light-field in a discontinuous way: the response at each pixel isthe convolution of the lightfield with a possibly multimodal function (this point is not
mentioned explicitly); this “point-spread function” is only sufficiently estimated to reconstruct images on a fixed plane, and calibration is done with controlled stimuli. In
this study however, we wish to exclude known calibration objects and other controlled
2 Except in rare cases such as a camera with a stereo adapter or split mirror.

2

stimuli.
Given the problem, a mixture statistical signal processing and geometry seems themost natural approach and this is indeed how we proceed.
1.1. Related work
Our approach to calibration is closely related to the work of Pierce and Kuipers [23],who measure the dissimilarity, or distance, between sensor elements that are not necessarily light sensors. The elements are then embedded in a metric space using metricscaling [17], which also determines the dimension of the space. A relaxation method
then improves this embedding, so that the Euclidean distance between sensor elements
better matches the dissimilarity between the sensor inputs. Getting very close to the
problem addressed in the present paper, the authors use this method to reconstitute thegeometry of a rectangular array of visual sensors that scans a fronto-parallel image.
Going further, Olsson et al. [21] use the information metric of [4] as a more appropriate method to measure the distance between visual or other sensor elements.
They also show how visual sensors -the pixels of the camera of a mobile robot- can bemapped to a plane, either using the method of [23], or their own, that embeds sensor
elements specifically in a square grid.
The works of Olsson et al. and of Pierce and Kuipers are very interesting to computer vision researchers, but they cannot calibrate an arbitrary discrete camera, since
the geometry of the embedding space is either abstract or fixed to a grid. In either case,it lacks an explicit connection to the geometry of the sensor. Filling this gap is one of
the motivations for the present work.
Since the present work exploits statistical properties of the light-field of the world
surrounding a light sensor, this article is also related to research on the statistical properties of real-world images [13]. In this field a theoretical framework which is gaining
more and more support considers the properties of the visual system to be reflections of
the statistical structure of natural images, because of evolutionary adaptation processes.
Our approach is similar in spirit, trying to replicate the plasticity of the visual systemof biological systems. That research has put in evidence fundamental properties, interms of local, global and spectral statistics, of real-world images, that has been also
exploited for computer vision tasks, such as classification [33], image restoration [7]and 3D inference [24].
Although these results are of great interest, they are not directly applicable in ourcase, mainly because we lack images.
Moreover, these statistics are about images formed on a planar image plane, which
is a hindrance in our case: first, we do not want to exclude the case of visual sensor
elements that are separated by more than 180 degrees, such as the increasingly popularomnidirectional cameras. Also, the local statistical properties of perspective images
depend of the orientation of the image plane with respect to the scene, except in special
constrained cases such as the fronto-parallel “leaf world” of Wu et al. [35]. Defining
images on the unit sphere thus appears as a natural way to render image statistics independent of the sensor orientation, at least with proper assumptions on the surrounding
world and/or the motion of the sensor.

3

Figure 2: The process of estimating the geometry of an unknown discrete camera.

Algorithm 1 Procedure for estimating the geometry of a discrete camera.
Input: Signals produced by N rigidly connected light sensors at T time instants:
x (i, t) , 1 ≤ i ≤ N , 1 ≤ t ≤ T . The light sensors are rigidly connected
and point towards unknown 3D directions X1 , . . . , XN , Xi ∈ R3 , X ⊤ X = 1.
Output: Estimates of the directions X1 , . . . , XN .
Algorithm:
1. Estimate the distance dij between the temporal signals x (i, .) and x (j, .),for each 1 ≤ i, j ≤ N .

2. Estimate the angular separations θij = cos Xi⊤ Xj .
3. Embed the angular separation in the 3D sphere: find X1 , . . . , XN s.t.
Xi⊤ Xi = 1, Xi⊤ Xj = cos (θij ), 1 ≤ i, j ≤ N .

1.2. Proposed approach
In the present work, we assume that the statistical properties of the data streams
produced by pairs of sensor elements depends only on the angular separation between
the elements. This assumption, when one considers the images as observations of arandom field defined on the sphere, is equivalent [26] to assuming that the random
field is homogeneous 3 - we could say “isotropic” to use the common image processing
vocabulary.
This assumption of homogeneity clearly does not hold in an anisotropic world,
unless the orientation of the sensor is uniformly distributed amongst all unitary transformations of the sphere. We thus adopt this last assumption, in order to ensure thehomogeneity of the observations. In terms of computer vision and robotics, our assumption amounts to saying that the sensor is randomly oriented, so that each photocell
is just as likely to sample the light-field in any direction.
The great practical utility of this assumption is that, as a consequence, the statistical properties of a pair of data streams generated by two photocells depend only on theangle separating them. In this situation, one can envision estimating the angular separa3 The definition of a homogeneous random field defined on the sphere is that the covariance between two
random variables sampled at two locations of the sphere depends only on the angular separation between thetwo points and that the expectation does not depend on the sampling location

4

Figure 3: Left: The camera used to sample omnidirectional images (image mirrored). Right: A calibrated
omnidirectional image mapped to a sphere.

tion from a measure of dissimilarity (e.g. correlation or information distance) between
the streams, and this is precisely what we do (Sec. 3). In order to estimate angles from
dissimilarity measures, we empirically observe the relation between these two quantities. We may then (Sec. 4.1) embed these angular estimates using techniques fromdistance geometry [5]. The whole process is outlined in Figure 2 and Algorithm 1. Finally Sec. 5 presents some conclusions and a list of questions to be addressed by future
research.
This work is the continuation of [11] and also includes some of the improvements
brought by our previous work [9]. In the latter, it is shown that a metric between
signals based on the correlation can just as well be used, instead of a metric based on
information distance. The embedding method used in [9] is described more thoroughly
in the present paper, allowing others to reproduce our results more easily.
2. Discrete camera model
Before entering into the details of our methodology for estimating the sensor geometry, we define the discrete camera and explain how to simulate it using an omnidirectional image sensor.
We define a discrete camera as a set of N photocells indexed by i ∈ {1, . . . , N },pointing in directions Xi ∈ R3 . T, on the number of pixels N , acquiring along time t,brightness measurements x (i, t) in the range {0, . . . , 255}. The directions of the lightrays, contrarily to conventional cameras, are not necessarily organized in a regular grid.
Many examples of cameras can be found under these definitions. One example is thelinear camera, where all the Xi are co-planar. Another example is the conventional perspective camera which comprises a rectangular grid of photocells that are enumerated
in our model by a single index i,nh i%W io⌋
Xi | Xi ∼ K −1 ⌊i/W
,
0
≤
i
<
HW
1

where W , H are the image width and height, K is the intrinsic parameters matrix, %represents the integer modulo operation and ⌊.⌋ is the lower-rounding operation. Cam5

eras equipped with fisheye lenses, or having log- polar sensors, can also be modeledagain by setting Xi to represent the directions of the light-rays associated to the image
pixels. In the same vein, omnidirectional cameras having a single projection center, as
the ones represented by the unified projection model [8], also fit in the proposed model.
In this paper we use a calibrated omnidirectional camera to simulate various discrete
cameras.
2.1. Simulated image sensor
We simulate a discrete camera with known Euclidean geometry by sampling a calibrated panoramic image with unique projection center at fixed locations. Since thecamera is calibrated, it is straightforward to locate the position (u, v) in the panoramic
image corresponding to the 3D direction Xi of a photocell that is part of the simulated discrete camera. In the present work, we use bilinear interpolation to measure thegraylevel value at non-integer coordinates (u, v).
Images are acquired by a VStone catadiopric camera consisting of a perspective
camera fitted to a hyperbolic mirror, shown in Figure 3, left. This system is modeled as
single projection center camera [8] with a 360◦ × 210◦ field of view and a ∼ 45◦ blindspot at the south pole (Fig. 3, right). The mirror occupies a 453 × 453 pixel region of
the image. The angular separation between neighboring pixels in the panoramic image
is usually slightly smaller than 0.5◦ . Also, some mild vignetting occurs, that could becorrected. Apart for these minor inconveniences, simulating a discrete camera by anomnidirectional camera presents many advantages: no other specialized hardware isneeded and each omnidirectional image can be used to simulate many discrete camera“images”, as in Fig. 4 (a). With respect to perspective cameras, the available field ofview allows to study very-wide-angle discrete cameras.
2.2. Data acquisition
For the purpose of studying the effect of angular separation on joint photocell signal
statistics, we use the 31-pixel planar camera (or “probe”) shown in Fig. 4 (b). Thisprobe design allows to study the effect of angular separations ranging from 0.5 to 180
degrees and each sample provides 465=31(31-1)/2 pixel pairs. In the “tighter” part ofthe discrete camera layout, there exists a slight linear dependence between the valuesof consecutive pixels due to aliasing.
The camera is hand-held and undergoes “random” general rotation and translation,according to the author´s whim, while remaining near the middle of the room, at 1.0
to 1.8 meters from the ground. We acquired three sequences consecutively, in very
similar conditions and joined them in a single sequence totaling 1359 images, i.e. approximately 5 minutes of video at 4.5 frames per second.
To simulate the discrete camera, we randomly choose an orientation (i.e. half agreat circle) such that all pixels of the discrete camera fall in the field of view of thepanoramic camera. Figure 4 (b) shows two such choices of orientations. For eachchoice of orientation, we produce a sequence of 31 samples x (i, t), 1 ≤ i ≤ 31, 1 ≤t ≤ 1359, where each x (i, t) ∈ {0, . . . , 255}. Choosing 100 different orientations, we
obtain 100 discrete sensors and 100 arrays of data. Appending these arrays we obtain
31 signals x (i, t) of length to 135900.

6

Sampled Directions
1
0.8
0.6
0.4
0.2
0
-1

(a) Two sensor instances

-0.5

0

0.5

1

(b) Linear discrete sensor(c) Four images from a sequence of 1359 for calibration learning.
Figure 4: (a)Two instances of the linear discrete camera, inserted in an omnidirectional image. Pixels locations are indicated by small crosses connected by white lines. (b) Geometry of a discrete camera consistingof a planar array of thirty one (31) pixels, spanning 180◦ in the plane. The first two pixels are separated by
0.5◦ , the separation between consecutive photocells increases geometrically (ratio ≃ 1.14), so that the 31st
photocell is antipodal with respect to the first. (c) Four images of a set to be used for calibration learning.

Considering any pair of pixel (indices) 1 ≤ i, j ≤ 31, the angular separation θij isthe main object to estimate by our method, and it is known in our test setup. We definein the next section two possible metrics d(i, j) between signals x(i, t) and x(j, t), 1 ≤t ≤ 135900 and, as a result, get a dataset, D of angle-distance pairs:
D = {(θij , d (i, j)) | 1 ≤ i, j ≤ 31} .

(1)

We will then use this dataset to build a functional relation that allows to estimate theangular separation θij from the information distance d (i, j).
3. Distances between pairs of signals vs angular separations
In this section, we define the measures of distance between signals, namely correlation and information distance, and use them to build a functional relation between theangle separations of photocells and the distances of pixel streams they acquire.
3.1. Correlation distance
We call correlation distance between signals x (t) and y (t), 1 ≤ t ≤ T , the quantitydc (x, y) =1
(1 − C (x, y)) ,27where C (x, y) is the correlation between the signals. It is easy to verify that dc (., .) isa distance.
For the task considered in this paper, it is natural to prefer the correlation distanceover the covariance or the Euclidean distance kx − yk, because both vary with signal
amplitude (and offset, for the latter), whereas dc (., .) is offset- and scale-invariant.
3.2. Information distance
Given two random variables x and y (in our case, the values produced by individual
pixels of a discrete camera) taking values in a discrete set {1, . . . , Q}, the information
distance between x and y is [4]:d (x, y)= H (x|y) + H (y|x) = 2H (x, y) − H (y) − H (x) ,(2)

where H (x, y) is the Shannon entropy of the paired random variable (x, y), and H (x)and H (y) are the entropies of x and y, respectively. It is easy to show that Eq. (2)defines a distance over random variables. This distance is bounded by H (x, y) ≤log2 Q, and can be conveniently replaced by the normalized information distance:d¯(x, y) = d (x, y) /H (x, y) .

(3)

The advantage of (3) over (2) is that its properties are more invariant to Q: Eq. (2) isbounded by log2 Q and is not necessarily equal to its upper bound when x and y are
independent. In contrast, Eq. (3) is bounded by 1 independently of Q and is exactly 1
if and only if x and y are independent [4].
3.3. Estimating the information distanceCaution should be taken when estimating the information distance (3) from finite
samples x (t), y (t), 1 ≤ t ≤ T : it is relatively common knowledge that replacing
unknown probabilities px (q) by sample frequencies pˆx (q) = |{t|x (t) = q}| /T 4 inthe expressions of the entropy results in the biased plug-in estimator Ĥ (x), with expectancy
P 1
 
n o
Q − 1 1 − q px (q)1
E Ĥ = H −
.
(4)
+
O
+
2T
12T 2
T3

This bias in turn causes a bias in the information distance estimates in Eq. (3). Whilecorrecting for the first bias term (Q − 1) /2T , is easy and results in the Miller-Madow
estimator, correcting for the other terms is more delicate [30]. In addition to the plug-in
and Miller-Madow estimators, we also consider here the estimator of Paninski [22].
A first consideration in the choice of estimators is the number of bins Q. Since thebias decreases slowly and is to a large extent proportional to Q, it is advantageous tore-quantize the signals from 256 graylevels to a more parsimonious representation withQ < 256 bins. When re-quantizing signals for estimating their information distance,
we found a slight advantage in choosing bins that maximize the entropy, i.e. bins that
contain equal numbers of values. All signal re-quantization in this paper is done using
this scheme, also used in [21].
4 |.| denotes the set cardinal.

8Mean Information Distance, Q=16

0.98
0.96
0.94
0.92

Mean Information Distance, Q=64

Plug-In
Miller-Madow
Paninski

0.98

Information DistancePlug-In
Miller-Madow
Paninski
Information Distance

Information Distance

Mean Information Distance, Q=4

0.96
0.94
0.92

2500

5000

10000

Signal Length

20000

40000

Plug-In
Miller-Madow
Paninski

0.98
0.96
0.94
0.92

2500

5000

10000

Signal Length

20000

40000

2500

5000

10000

20000

40000

Signal Length

Figure 5: Each of these plots shows how bias of various estimators is affected by signal length. Each curve isthe average of 150 information distance estimates, taken between signals generated by photocells separated
by approximately 5 degrees. The signal length is varied from 1359 to 43488. A good estimator is one
for which the mean value changes little with the signal length. The three considered information distanceestimators are: the plug-in estimator, the Miller-Madow corrected estimator and Paninski´s estimator.

In order to choose a suitable estimator of the information distance, we first consider
the stability of each estimator with respect to signal length. In Fig. 5, each curve is theaverage of 150 information distance estimates, taken between signals generated by two
photocells separated by approximately 5 degrees, when the signal length varies approximately from 1350 to 435005. The leftmost curve shows that, for Q = 4, all estimators
behave pretty well, although the asymptotic value appears smaller than with Q > 4(two plots at right). This last effect results from the information lost in quantization.
For Q = 16, the middle curve shows that Paninski´s estimator is stablest, but not much
more than the Miller-Madow estimator, while the plug-in estimator is clearly far fromits asymptotic value when the signals length is 1350. Finally, for Q = 64, the rightmost curve shows the trend of the previous curve greatly increased: as Q increases,
the term Q2 in the bias expression for H (X, Y ) or H (X|Y ) in Eq. (2) becomes more
important and more samples are needed to reduce it. Not shown in these curves, thestandard deviations of all three estimators are comparable.
From the previous result, the Miller-Madow estimator appears like a good choice,
due to its relatively good performance, and to the fact that it is is very easy to compute.
In contrast, the Paninski estimator requires the computation of a table of same length
as the signal, which renders it impractical with signals of variable length.
In order to further assess the advantage of the Miller-Madow estimator over theplug-in estimator, we show the results of other experiments performed by varying Q
and the angular separation between the photocells.
Figure 6 shows the average information distance Eq. (3) between the outputs ofsensors with known angular separation. Each plot holds four curves, corresponding to
Q = 2, 4, 16 and 64. The left plot uses bias correction, whereas the right plot doesnot. The curves in the left plot are better grouped, showing that, when bias reduction isapplied, the information distance estimates depend less on Q. This is important if one
is interested in estimating properties of the information distance that depend on angularseparation and on the ambient light field, rather than on Q and T .
5 In order to obtain signals of length greater than 1359, the original length of our signal, we concatenate
many signals corresponding to different sensor orientations.

9Mean Information Distance vs Angle

Mean Information Distance vs Angle
1

Info. Distance (Bias-Corrected)

Info. Distance (Uncorrected)

1
0.95
0.9
0.85
0.8
0.75
Q=2
Q=4
0.5

Q=16
Q=64

0.25
0
0

10

20 30

60

90

120 150 180

0.95
0.9
0.85
0.8
0.75
Q=2
Q=4
0.5

Q=16
Q=64

0.25
0
0

Angular Separation (deg)

10

20 30

60

90

120 150 180

Angular Separation (deg)Figure 6: Dependency of the information distance d (i, j) on the angular separation θij . The height of eachplotted point is the average of 100 estimates of the information distance between signals x (i, 1 ≤ t ≤ 1359)
and x (j, 1 ≤ t ≤ 1359). The abscissa is the angular separations between pixels i and j of the discrete
camera. Each curve corresponds to a different number Q of quantization bins. The quantization bins are all
taken to have weight 1/Q. Right: Information distances computed with the first-order bias correction termof Eq. (4) removed. Left: Distances computed without correction.

The conclusion of these experiments is that adding the first correction term in
Eq. (4) to entropy estimates in Eq. (3) is justified and beneficial, especially when only
short sequences are available. For the above reasons, we will use the Miller-Madow
estimator in the sequel, and use Q = 4 bins, which is appropriate for short sequences.
One should note that correcting for bias in the entropy estimates is less importantin the situation of Olsson et al. [21], where the distance measurements are embedded
directly in an abstract space. Obtaining accurate information distance estimates is important in our approach because our goal is to accurately map information distance toangles.
3.4. Estimating angular separations from inter-signal distances
As explained earlier, our a-priori knowledge of the world will be encoded in a functional relation mapping a measure of discrepancy between two signals, to the angularseparation between the photocells that generated the signals. We now build this functional relation and assess its effectiveness at estimating angles.
In the previous section we detailed the computation of the distances between pixel
streams, and therefore completed the construction of the angles and distances dataset
D (Eq. 1). From this dataset, we can finally build a constant by parts model of the
expectancy of the distance knowing the angle. After verifying and, if needed enforcing,the monotonicity of this model, we invert it, obtaining a graph of angles, θ as a functionof (correlation or information) distances, d:θ = F (d)(5)

Strict monotonicity has to be enforced for the correlation-based data, owing to therelatively small number of data points used for each quantized angle.
Figure 7 shows the resulting graphs. This figure shows one of the major issues
that appear when estimating the angular separation between pixels from the correlation
or information distance: the graphs become very steep for large values of the distance,
10

Information-Distance-to-Angle Model
120

100

100
Predicted Angle

Predicted Angle

Correlation-Distance-to-Angle Model
120

80
60
40
20

80
60
40
20

0

0
0

0.1

0.2
0.3
Correlation Distance

0.4

0.5

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Information Distance1

Figure 7: Models relating correlation (left) or information distance (right) to angular separation between
photocells. These models were build from simulated signals produced by the linear probe of Fig. 4(a).
Signals of length T = 135900, acquired indoors were used.

indicating that small changes of the distance result in large changes in the estimated angle. Hence, a small uncertainty in the estimated distance implies a large uncertainty inthe estimated angle. This problem is somewhat mitigated, in the case of the correlation
distance, by limiting the abscissa to values [0, 1/2]. For small distance values, on theother hand, the curves are much flatter, suggesting that small angles can be determinedwith greater accuracy. Both trends are particularly true for the information distance.
3.5. Experimental validation
We now assess how well angles can be estimated from the graphs obtained in theprevious section. As described in Sec. 2.2, we collect in each omnidirectional image the samples necessary to simulate 100 discrete linear sensors, each one having 31
photocells (100 times 31-tuples per image). For each 31-tuple of signals, we estimate
the 31 · 30/2 correlation and information distances between pairs of signals and map
these distance to angles using the models of Fig. 7. We assess the quality of the angle
estimators by comparing these estimated angles with the known true angles.
Figure 8 (a) and (b) show the boxplots the 100 estimated angles versus the true
angles for the correlation- and information-based estimators, respectively. These plotsshow that the true angles are well estimated for small angles, while the spread andapparent bias of the estimators increases sharply for angles greater than five or ten
degrees.
Figure 8 (c) and (d) better illustrate this trend: Figure 8 (c) plots the mean estimated angle versus the true angle. This figure shows that the estimators are essentiallyunbiased until five degrees and that they underestimate angles after that. Figure 8 (d)plots the mean absolute error in the estimated angles.
Figure 8 shows how the proposed angle estimators perform on signals producedin the same conditions of environment and motion as the signals used to build thecorrelation and information distance-to-angle models. Since we will want to use theseestimators in other situations, we must also evaluate their performance in more generalconditions.

11Angles from Correlation Distance
180

150

150
Estimated Angular Separation

Estimated Angular Separation

Angles from Information Distance
180

120
90
60
30
20
10
5
1
0

120
90
60
30
20
10
5
1
0

0

1

5

10

20

30

60

90

120

150 180

0

1

5

10

True Angular Separation

20

30

60

90

120

150 180

True Angular Separation(a) Mean estimated angle v.s. true angle(b) Mean absolute error in estimated angles

Mean Angle Estimate from Signal Distance180

Abs. Error in Estimated Angles
Mean Est. from Info. Distance
Mean Est. from Correlation

160

90
Mean Absolute Error (deg)

Estimated Angle

Info. DistanceCorrelation

80

140
120
100
80
60
40
20

70
60
50
40
30
20
10

0

0
0

20

40

60

80 100 120 140 160 180
True Angle0(c) Mean estimated angle v.s. true angle

2040

60 80 100 120 140 160 180
True Angle (deg)(d) Mean absolute error in estimated angles

Figure 8: Validation of the distance-to-angle models of Figure 7. Boxplots (a) and (b) show the 5th percentile,
first quartile, median, third quartile and 95th percentile of 100 angles estimated from signals of length T =1359 simulated using the indoor sequence of Fig. 4 (c). The angles are estimated from information distancein (a) and correlation in (b). The curves in (c) show the mean estimated angle v.s. the true angle, and (d)shows the mean absolute error in the estimated angles.

Figure 9 shows how the proposed angle estimators perform on signals captured
outdoors and indoors by the same camera setup as above. This time, the camera isheld mostly horizontally, but the virtual discrete cameras still have pretty random orientations, owing to the very wide field of view. Four of the 2349 images are shown in
Fig. 9 (a). Again, we sampled 100 time 31-tuples of pixel-values per image and compared the true angles with the angles estimated using the correlation and information
distance based estimators.
Figure 9 (b) shows the mean of the 100 estimated angles, while Fig. 9 (c) shows themean absolute error. This figure shows that, globally, the estimators actually performbetter on this dataset than on the dataset acquired in the same conditions used to build
the estimators.
Using the same outdoors and indoors sequence, but sampling horizontal virtual
cameras allows to evaluate the performance of the angle estimators when the assumption of uniform sensor orientation is violated. Since the camera was mostly horizontal
during the capture of the outdoors and indoors sequence of Fig. 9, and the virtual sensors sampled in the present experiment remained horizontal, these virtual sensors in
effect have a very restricted set of orientations.
12(a) Four images of a sequence for testing the calibration methodology.
Mean Angle Estimate from Signal Distance180

Abs. Error in Estimated Angles
Mean Est. from Info. Distance
Mean Est. from Correlation

160

90
Mean Absolute Error (deg)

Estimated Angle

Info. DistanceCorrelation

80

140
120
100
80
60
40
20

70
60
50
40
30
20
10

0

0
0

20

40

60

80 100 120 140 160 180
True Angle0(b) Mean estimated angle v.s. true angle

2040

60 80 100 120 140 160 180
True Angle (deg)(c) Mean absolute error in estimated angles

Figure 9: (a) Four out of 2349 images of an outdoors and indoors sequence used to validate the angle
estimators. (b) Mean of 100 sets of angles estimated from this sequence, plotted versus the true angle. (c)Mean absolute error in the estimated angles.

Nevertheless, Fig. 10 shows that the signals captured in these conditions allow toestimate angles with an accuracy that is less good than that of Fig. 9, but not much lessgood. One should note that the mean estimated angle is always smaller than the true
angle. This seems to indicate that the correlation is greater between signals obtainedby sampling along the horizon.
Finally we go a step further in testing the generalization ability of the proposed
angle estimators. We now use different sensor, point-and-shoot Olympus Stylus 300
camera that can acquire 160 × 120 videos at ∼ 15 frames per second. Multiple sequences of typically 1050 images (70 seconds) were taken indoors and outdoors, while
the camera was moving, predominantly horizontally and forward, but also with different orientations. The horizontal and vertical fields of view cover 50 and 35 degrees,
respectively. These sequences were concatenated, forming a single sequence of 22822
images. Figure 11 (a) show four of the images in the sequence. Obviously, this sequence does not allow to simulate a sensor with 180 degrees of field of view. In thisexperiment, we use a linear sensor consisting of 31 photocells spanning 35 degrees.
The first two pixels are separated by 0.5◦ , the separation between consecutive photocells increases geometrically, so that the 31st photocell is 35 degrees away from thefirst.
Figure 11 (b) and (c) show the precision and accuracy of the estimated angles. This
figure shows that both estimators slightly overestimate angles, possibly due to the hightexture contents that is often present in the images. The estimator based on information

13

Mean Angle Estimate from Signal Distance180

Abs. Error in Estimated Angles
Mean Est. from Info. Distance
Mean Est. from Correlation

160

90
Mean Absolute Error (deg)

Estimated Angle

Info. DistanceCorrelation

80

140
120
100
80
60
40
20

70
60
50
40
30
20
10

0

0
0

20

40

60

80 100 120 140 160 180
True Angle0(a) Mean estimated angle v.s. true angle

2040

60 80 100 120 140 160 180
True Angle (deg)(b) Mean absolute error in estimated angles

Figure 10: Results of the angle estimators estimated from the sequence of Fig. 9, but with the virtual linear
sensor kept mostly horizontal. (a) Mean of 100 sets of angles, plotted versus the true angle. (b) Meanabsolute error in the estimated angles.

distance performs somewhat better than that based on correlation.
4. Calibrating a discrete camera
Having seen how to estimate the angle between two pixels from the correlation or
the information distance between their respective data streams, we now estimate thewhole sensor geometry from angle estimates. In this section, we produce sequencesof pixel signals in the same conditions as in Sec. 2.2, except that the sensor shape isdifferent. The distance between all pairs of signals are then estimated, and the angular separation between the pixels are estimated using Sec. 3.4. Finally, these angle
estimates are embedded in the sphere using the algorithms defined in this section.
4.1. Embedding points in the sphere
The last step in order to solve the discrete camera calibration problem will be taken
by solving the problem:
Problem 1) Spherical embedding problem: Given angle estimates θij , 1 ≤ i, j ≤N , find points Xi on the unit sphere, separated by angles approximately equal toθij , i.e. Xi⊤ Xj ≃ cos θij , for all i, j.
This problem can be reduced to the classical problem of distance geometry [5]:Problem 2) Euclidean embedding problem: Given distance estimates Dij , 1 ≤ i, j ≤N , find points Yi in a metric vector space, such that, for all i, j, kYi − Yj k ≃ Dijp
Indeed, by defining an extra point Y0 = (0, 0, 0), and distances Dij = 2 − 2 cos θij ,the mapping of the first problem to the second is immediate.

14

(a) Four images from a sequence of 22822 for testing the calibration methodology.
Mean Angle Estimate from Signal Distance

35
Mean Absolute Error (deg)

Estimated Angle

Abs. Error in Estimated Angles

Mean Est. from Info. Distance
Mean Est. from Correlation

50
40
30
20
10
0

Info. DistanceCorrelation

30
25
20
15
10
5
0

0

5

10

15
20
True Angle

25

30

35

0

(b) Mean estimated angle v.s. true angle

5

10

15
20
25
True Angle (deg)3035(c) Mean absolute error in estimated angles

Figure 11: (a) Images 300, 8400, 11300 and 13100, in a set of 22822 images acquired with an off-the-shelfdigital camera. (b) Mean of 100 sets of angles, plotted versus the true angle. (c) Mean absolute error in theestimated angles.

4.2. Factorization method
Solutions to both problems (with exact equality, rather than approximate) were
published by 1935 [29]6 . Schoenberg´s Theorem 2 [29] states that if the matrix C withterms Cij = cos θij is positive semidefinite with rank r ≥ 1, then there exist points on
the unit (r − 1) −dimensional sphere that verify Xi⊤ Xj = Cij for all i, j. This result
directly suggests the following method for embedding points in the 2-sphere:
1. Build the matrix C with terms Cij = cos θij , 1 ≤ i, j ≤ N .
2. Compute, using the SVD decomposition, the rank-3 approximation C̃ = U U ⊤of C, where U is N × 3. In essence, this involves square-rooting the three largest
singular values of C and using them to scale the first three components of theleft singular values 7 .
3. Define Xi = (Ui1 , Ui2 , Ui3 ) / k(Ui1 , Ui2 , Ui3 )k.
In the absence of noise, this very simple algorithm gives the exact solution, up to a unitary transformation. When there is noise, however, it is not optimal in many ways. In
particular, it does not take into account the structure of the approximately multiplicative
noise in the angles θij .
6 Schoenberg cites [15] and [19] as, respectively, previous solutions to Problems 1 and 2. These are theearliest references we are aware of.
7 Using Matlab/Octave notation:[u,s,v]=svd(C); U=u(:,1:3)*sqrt(s(1:3,1:3))

15

Sensor diameter = 5◦

Sensor diameter = 10◦

Sensor diameter = 20◦

Figure 12: True and estimated pixel layouts of a discrete camera consisting of photocells lying on a triangulargrid. To estimate the layout, each pixel was sampled 1359 times. The information distance between pairs of
pixels was estimated from these values, and converted to angular estimates using our non-parametric model.
The angular distances were then embedded in the sphere. For visualization, they are aligned by the usual
procrustes method, mapped to the plane by projective mapping with unit focal length, and line segments
indicate the original pixel neighborhood relations.

In practice, this algorithm is very sensitive to noise: Fig. 12 shows the results of
the reconstruction, when the pixels lie on a triangular grid of the image plane, scaled
so that the angular diameter of the whole discrete camera is 5◦ (left), 10◦ (middle) or20◦ (right). Please note that the lines connecting the photocell positions serve only tobetter illustrate the correspondence.
The reconstructed grids shown here are typical of what researchers reproducing
our method may encounter. It is easy to see that for small-diameter discrete cameras,
the estimated shape is pretty well determined by the information distance between the
pixel streams, whereas, when the angular diameter increases, the noise in the estimatedangles gradually overwhelms the embedding algorithm. It is common for the factorization method to break down with sensors of 20◦ diameter or less [11].
A different method is thus needed if we wish to obtain better calibration results. It iseasy to verify that this situation is not directly tractable by variable-error factorization
methods used in computer vision, such as [14]. In the next section, we will define andvalidate a more performant algorithm.
4.3. Robust nonlinear embedding method
Noting that the error in the estimated angles is approximately proportional to theactual angle suggests that the embedding method should weigh less heavily large angular estimates.
One such method is Sammon´s algorithm [27], which we adapt to spherical -rather
than Euclidean- embedding. In this paper, we minimize the weighted sum of the differences of the angle estimates with the internal products of directions:X̂ = arg min
XX

wi,j Xi⊤ Xj − Cij

i,j

16

2

,

(6)

Algorithm 2 Spherical embedding
Input: Distance matrix Θ with terms θij , 1 ≤ i, j ≤ N, θij ∈ [0, π]
Output: Embedded points Xi , 1 ≤ i ≤ N on the unit sphere.
1. Compute the minimum spanning tree T of Θ.
2. Select the triplet (i, j, k) of points such that (i, j) ∈ T , (i, k) ∈ T and theˆ is maximized.
spherical angle jik
3. Choose Xi , Xj and Xk s.t. Xi⊤ Xj = cos (θij ), Xi⊤ Xk = cos (θik ) andXj⊤ Xk = cos (θjk ).
4. Define the set of already embedded points B = {i, j, k}.
5. Define the set N of neighbors of B that are not in B: N ={i | ∃j ∈ B, (i, j) ∈ T } .
6. For each element i in N , taken from the closest to B to the furthest:
7.
Find j ∈ B s.t. (i, j) ∈ T and k ∈ B that maximizes the sphericalˆ
angle jik.
Triangulate Xi from θij and θik (disambiguate using other points of

8.
B).

9.
Update N := N \ {i}, B := B ∪ {i}.
10. End while
11. If B =6 {1 . . . N }, then goto 5.

where X = {X1 , X2 , · · · , XN } and
n
o
(
1
1
max 0, 1−C
− 1−C
ij
o
wij =
1
η

if Cij 6= 1
otherwise.

To reflect the fact that big angles are less well estimated, we set C0 = 0.9, so that
estimates greater than acos (0.9) ≃ 25◦ be ignored. The other parameter, η is set
to 1, allowing the points Xi to stray a little bit away from the unit sphere. Our implementation is inspired by the second-order iterative method of Cawley and Talbot(http://theoval.sys.uea.ac.uk/~gcc/matlab/default.html). For initialization, we modify the algorithm by Lee, Slagle and Blum [18], which operates on theplane, so that it operates on the sphere.
Combining Sammon´s algorithm to that of Lee, Slagle and Blum is not an unusual
practice for embedding in Euclidean space, e.g. [2].
The algorithm by Lee, Slagle and Blum works by steps, adding one point to aset of already embedded points at each step, starting from a triangle. The embedded
points verify 2N − 3 of the original distances θij exactly, i.e. the minimum numberof distances that uniquely define the positions of all N points up to an isometry. The
algorithm has quadratic complexity in N .
In our case, the algorithm is adapted to embed in the unit sphere an N × N distancematrix with terms θij , 1 ≤ i, j ≤ N , θi,j ∈ [0, π]. The algorithm is detailed in Algorithm 4.3. It is easy to see that this algorithm favors using small angle estimates and17

Figure 13: Calibrations of two different sensors covering more than one hemisphere. On the left, a band-likesensor consisting of 85 photocells, calibrated from correlations (estimated: smaller, true: bigger). On theright, a discrete camera covering more than 180×360◦ , of 168 photocells, calibrated from the information
distance (estimated: smaller, true: bigger). Each ball represents a photocell except the big black balls,
representing the optical center.

ignoring larger ones: first, in step 1, the edges that are part of the minimum spanning
tree are by nature short (for example, the nearest-neighbor graph is a subgraph of theminimum spanning tree). Then, when adding a point to the set B of already embedded
points (step 7), this is done by triangulating with two points of B that are necessarily
ˆ in step 7 cannot be greater thanclose to the added point. Also note that the angle jik
π/3, by definition of the minimum spanning tree, and the triangle used to localize the
newly added point is far from a degenerate configuration. In addition to the usefulproperties just mentioned, this algorithm is simple to implement.
4.4. Testing the generality of the method
We now evaluate the results of the complete calibration methodology, and in particular the robust embedding algorithm on data produced by the angle-estimating methodof Sec. 3.4. In this section we test two main cases: (i) calibrating non planar discrete
cameras having enlarged fields-of-view, using the same image acquisition device described in Sec. 2.1, and (ii) calibrating a rectangular discrete camera, but based on adifferent image acquisition device, namely a (consumer) digital camera.
Calibration from an indoors-and-outdoors sequence. The main objective here is testing the calibration of large fields-of-view while using calibration images that have completely different illuminations and textures (as compared to indoors image datasets).
For this purpose, we produce sequences of pixel signals in the same conditions as previously, using the outdoors and indoors sequence shown in Fig. 9(a), except that thediscrete-sensor shape is different. The information and correlation distances between
pixels is then estimated from these signals, the angular separation between the pixels
is estimated using Sec. 3.4, and the embedding method of Sec. 4.1 is applied to theseangle estimates.
Figure 13 shows the results of our calibration method on sensors covering morethan a hemisphere, which thus cannot be embedded in a plane without significant distortion. It should be noted that, although the true sensor is each time more than hemispheric, the estimated calibration is in both cases smaller. This shrinkage is a known
effect of some embedding algorithms, which we could attempt to correct.

18

(a) Correlation-Based(b) Information Distance-Based

Figure 14: Reconstructed and true pixel layouts of a discrete camera consisting of photocells lying on arectangular grid (a,b). The sensor used differs from that with which the models of Fig 6 were built. The
reconstructions are obtained by first estimating the pairwise angular distances, then embedding the angles
in the sphere (see text). For visualization, the reconstructions are aligned by the usual procrustes method,
mapped to the plane by projective mapping with unit focal length. Added line segments show the true pixel
neighborhood relations. Plot (a) is obtained from the correlation distance, and plot (b) from the information
distance.

Calibration using a different image sensor. Figure 14 shows how our method applies
to signals produced by a different sensor from the one used to build the distance-toangle models, namely an Olympus Stylus 300 camera. An 8-by-8 square grid pixels
spanning 34 degrees was sampled along a 22822 image sequence taken indoors andoutdoors. From this sequence, the estimated angles were generally greater than thetrue angles, which explains the absence of shrinkage. The higher angle estimates werepossibly due to higher texture contents of the sequence. The estimated angles were also
fairly noisy, possibly due to the sequence length, and we surmise that longer sequences
would yield better results.
These results represent typical results that researchers reproducing our method may
encounter.
5. Discussion
This paper addressed the problem of determining the geometry of a set of photocells
in a very general setting. We have shown that a discrete camera can be calibrated toa large extent, using just two pieces of data: a table relating information distances
to angles; and a long enough signal produced by the camera. The main assumptions
are that the camera motion directs each pixel uniformly in all directions, and that theenvironments in which the table is built and the one in which the calibrated sensor
evolves are statistically similar.
Our algorithm proceeds in three main computational steps. We first estimate the information distance between each pair of signals. Then, using a pre-built table relating
distances to angles, we transform the information distances into estimated angles. In
the final step we embed the angles in a sphere, with a method that takes into account thenoise introduced by the estimation process. It should be noted that the first and third
steps are in themselves open problems that are the object of much research [1, 3, 5].
In particular we adapted two embedding methods [18, 27] and obtained an algorithm

19adapted for the approximately multiplicative noise in the estimated angles. Experiments on simulated sensors assess the validity of the proposed approach.
As said above, we assumed an isotropic light field, which can be easily produced
manually moving the camera in any environment with enough visual stimuli. Dropping this assumption might allow to calibrate, for example, a sensor that is moving
mostly in one direction e.g. the eyes of a fly that mostly flies forward. It could also beinteresting to generalize the proposed method to calibrate a discrete camera with a noncentral calibration model, e.g. an earthworm with eyes all over the body. This poses
an immediate challenge, because describing the geometric relation between two nonintersecting rays requires four parameters, while a single angle describes the relation
between two intersecting rays. Future work will address this open problems.
References[1] J. Beirlant, E. Dudewicz, L. Gyorfi, and E. van der Meulen. Nonparametric entropy estimation: An overview. International J. of Mathematical and Statistical Sciences, 6:17–39,
1997.
[2] G. Biswas, A.K. Jain, and R.C. Dubes. Evaluation of projection algorithms. IEEE Transactions on Pattern Analysis and Machine Intelligence, 3(6):701–708, 1981.
[3] P. Biswas, T.-C. Liang, K.-C. Toh, Y. Ye, and T.-C. Wang. Semidefinite programming approaches for sensor network localization with noisy distance measurements. IEEE Transactions on Automation Science and Engineering, 3(4):360–371, 2006.
[4] J. P. Crutchfield. Information and its metric. In L. Lam and H. C. Morris, editors, Nonlinear Structures in Physical Systems–Pattern Formation, Chaos and Waves, pages 119–130.
Springer-Verlag, 1990.
[5] Jon Dattorro. Convex Optimization & Euclidean Distance Geometry. Meboo Publishing,
2005.
[6] R. Fergus, A. Torralba, and W. T. Freeman. Random lens imaging. Technical Report MIT
CSAIL TR 2006-058, Massachusetts Institute of Technology, 2006.
[7] W. T. Freeman, E. C. Pasztor, and O. T. Carmichael. Learning low-level vision. International Journal of Computer Vision, 40(1):25–47, 2000.
[8] C. Geyer and K. Daniilidis. A unifying theory for central panoramic systems and practicalapplications. In proc. ECCV, volume II, pages 445–461. Springer-Verlag, 2000.
[9] E. Grossmann, J. A. Gaspar, and F. Orabona. Calibration from statistical properties of thevisual world. In European Conf. on Computer Vision, 2008, 2008.
[10] E. Grossmann, E-J Lee, P. Hislop, D. Nistér, and H. Stewénius. Are two rotational flows
sufficient to calibrate a smooth non-parametric sensor? In proc. IEEE CVPR, 2006.
[11] E. Grossmann, F. Orabona, and J. A. Gaspar. Discrete camera calibration from the information distance between pixel streams. In Proc. Workshop on Omnidirectional Vision,
Camera Networks and Non-classical Cameras, OMNIVIS, 2007.

20

[12] R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision. Cambridge
University Press, 2000.
[13] A. Hyv arinen, J. Hurri, and P. O. Hoyer. Natural Image Statistics — A probabilistic approach to early computational vision. to be published by Springer-Verlag, 2008. Preprint
available at http://www.naturalimagestatistics.net/.
[14] M. Irani and P. Anandan. Factorization with uncertainty. In European Conf. on Computer
Vision, 2000.
[15] L. Klanfer. Metrische charakterisierung der kugel. In Ergebnisse eines mathematischen
kolloquiums, volume 4, pages 43–45, 1933.
[16] I. Kohler. Experiments with goggles. Scientific American, 206:62–72, 1962.
[17] W. J. Krzanowski. Principles of Multivariate Analysis: A User’s Perspective. Clarendon
Press, Statistical Science Series, 1988.
[18] R. C. T. Lee, J. R. Slagle, and H. Blum. A triangulation method for the sequential mappingof points from n-space to two-space. IEEE Trans. Computers, 26(3):288–292, 1977.
[19] K. Menger. New foundations of euclidean geometry. American Journal of Mathematics,
53:721–745, 1931.
[20] D. Nistér, H. Stewenius, and E. Grossmann. Non-parametric self-calibration. In proc.
ICCV, 2005.
[21] L. Olsson, C. L. Nehaniv, and D. Polani. Sensory channel grouping and structure fromuninterpreted sensor data. In NASA/NoD Conference on Evolvable Hardware, 2004.
[22] L. Paninski. Estimation of entropy and mutual information. Neural Computation, 15:1191–
1254, 2003.
[23] D. Pierce and B. Kuipers. Map learning with uninterpreted sensors and effectors. Artificial
Intelligence Journal, 92(169–229), 1997.
[24] B. Potetz and T. S. Lee. Scaling laws in natural scenes and the inference of 3d shape. In
NIPS – Advances in Neural Information Processing Systems, pages 1089–1096. MIT Press,
2006.
[25] S. Ramalingam, P. Sturm, and S. Lodha. Towards complete generic camera calibration. In
Proc. CVPR, volume 1, pages 1093–1098, 2005.
[26] R. Roy. Spectral analysis for a random process on the sphere. Annals of the institute ofstatistical mathematics, 28(1), 1976.
[27] J. W. Jr. Sammon. A nonlinear mapping for data structure analysis. IEEE Transactions on
Computers, C-18:401–409, 1969.
[28] G. Schmidt and D. T. Moore. Tapered gradient index microlenses for compound lens arrays.
In Proc. SPIE Vol. 6342, International Optical Design Conference, 2006.
[29] I. J. Schoenberg. Remarks to Maurice Fréchet’s article “Sur la définition axiomatique d’une
classe d’espaces distanciés vectoriellement applicable sur l’espace de Hilbert”. Annals of
Mathematics, 36(3):724–732, 1935.

21

[30] T. Schuermann. Bias analysis in entropy estimation. J. Phys. A: Math. Gen, 37:L295–L301,
2004. arXiv:cond-mat/0403192v3.
[31] F. W. Snyder and N. H. Pronko. Vision with spatial inversion. University of Wichita Press,
1952.
[32] G. M. Stratton. Some preliminary experiments on vision without inversion of the retinal
image. Psychological Review, 3(6):611–617, Nov 1896.
[33] A. Torralba and A. Oliva. Statistics of natural image categories. Network: Computation in
Neural Systems, 14:391–412, 2003.
[34] R. Tsai. An efficient and accurate camera calibration technique for 3D machine vision. In
IEEE Conf. on Computer Vision and Pattern Recognition, 1986.
[35] Y. N. Wu, S.-C. Zhu, and C.-E. Guo. From information scaling of natural images toregimes of statistical models. Technical Report 2004010111, Department of Statistics,
UCLA, 2004.

22