Antoine Deleforge

Radu Horaud

2D SOUND-SOURCE LOCALIZATION ON THE BINAURAL MANIFOLD

The problem of 2D sound-source localization based on a robotic binaural setup and audio-motor learning is addressed. We first introduce a methodology to experimentally verify the existence of a locally-linear bijective mapping between sound-source positions and high-dimensional interaural data, using manifold learning. Based on this local linearity assumption, we propose an novel method, namely probabilistic piecewise affine regression, that learns the localization-tointeraural mapping and its inverse. We show that our method outperforms two state-of-the art mapping methods, and allows to achieve accurate 2D localization of natural sounds from real world binaural recordings.

INTRODUCTION

The human extraordinary ability of localizing one or several sound sources from the perceived acoustic signals has been intensively studied in cognition [START_REF] Blauert | Spatial Hearing: The Psychophysics of Human Sound Localization[END_REF], in computational audition [START_REF] Wang | Computational Auditory Scene Analysis: Principles, Algorithms and Applications[END_REF], and more recently in the emerging field of robot hearing [START_REF] Deleforge | The cocktail party robot: Sound source separation and localisation with an active binaural head[END_REF]. The most commonly used cues in binaural computational sound-source localization (SSL) are interaural cues, namely the ITD (interaural time difference), the IPD (interaural phase difference) and the ILD (interaural level difference). A lot of techniques exist to evaluate their values, either in the time domain using cross-correlation [START_REF] Liu | Azimuthal source localization using interaural coherence in a robotic dog: modeling and application[END_REF], or in the time-frequency domain using Fourier analysis [START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF][START_REF] Deleforge | A latently constrained mixture model for audio source separation and localization[END_REF] or gammatone filters [START_REF] Woodruff | Binaural localization of multiple sources in reverberant and noisy environments[END_REF]. Once such cues are computed they need to be mapped to a source position. The vast majority of current SSL approaches mainly focus on frontal azimuth estimation, i.e., 1D localization [START_REF] Liu | Azimuthal source localization using interaural coherence in a robotic dog: modeling and application[END_REF][START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF][START_REF] Woodruff | Binaural localization of multiple sources in reverberant and noisy environments[END_REF][START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF], while very few perform 2D localization [START_REF] Kullaib | 2d binaural sound localization: for urban search and rescue robotics[END_REF]. All these approaches rely on a simplifying geometric or parametric model of sound propagation in the binaural system. Simple models assume a direct path propagation from source to microphones (single attenuation coefficient and delay) while more advanced ones use a sphericalhead model for ITD propagation (Woodworths formula), approximate ILD data from a human head related transfer function (HRTF) dataset with a sine function [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF] or rely on a spiral ear model [START_REF] Kullaib | 2d binaural sound localization: for urban search and rescue robotics[END_REF]. In all these cases, extra parameters are needed such as the distance between microphones, the radius of the head, the sound speed, the shape of the ear or the sine coefficients. Following this view, the number of parameters required to accurately model a real world binaural system may become prohibitively high including the exact shape of the recording device, of the room, and all their acoustic properties. This is often unaccessible in practice, which is an important limitation for real-world SSL. On the other hand, in the particular case of a hearing robot, other sensory data may be available, such as the motor states (proprioception) or source position based on vision.

The task of learning a mapping between two spaces can be summarized as follows: if we are given a set of training couples {(x n , y n )} N n=1 ⊂ R L × R D , how can we obtain a relationship between the latent space R L and the observation space R D such that given a new observation in one space, its associated point in the other space is deduced? This problem has been extensively studied in machine learning, and offers a broad range of applications. In audio, it was notably used in text-to-speech synthesis [START_REF] Kain | Spectral voice conversion for textto-speech synthesis[END_REF], voice conversion [START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF] or articulatory-acoustic mapping systems [START_REF] Toda | Statistical mapping between articulatory movements and acoustic spectrum using a gaussian mixture model[END_REF], but never thoroughly examined for SSL. In this study, we focus on a particular class of mapping techniques where the relationship between the two spaces is approximately locally linear. This approximation is particularly relevant for data lying on a smooth Riemanian manifold, which by definition is locally homeomorphic to an Euclidean space. Given a set of K unknown linear transformations, locally linear mapping can be split in (i) assigning training couples to transformations and (ii) using linear regression to estimate the parameters of each transformation. These two tasks can be achieved within a probabilistic framework, where observations are seen as the realization of random variables X and Y , while assignments to transformations are modeled as hidden variables Z. Although the mixture of linear regression model (MLR) [START_REF] De Veaux | Mixtures of linear regressions[END_REF] represents data by a set of linear transformations, the transformations are not local since Z is independent of X. This is unsuitable for manifold data. Mapping methods based on Gaussian mixture models (GMM) [START_REF] Kain | Spectral voice conversion for textto-speech synthesis[END_REF][START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF][START_REF] Toda | Statistical mapping between articulatory movements and acoustic spectrum using a gaussian mixture model[END_REF], widely used in audio applications, were recently unified by the mixture of probabilistic linear regression model (MPLR) [START_REF] Qiao | Mixture of probabilistic linear regressions: A unified view of GMM-based mapping techiques[END_REF]. In MPLR, posteriors p(Z|x) or p(Z|x, y) are first obtained by estimating a GMM with the EM algorithm on latent data or joint data. Then, transformations parameters are estimated with weighted lin-ear regression. As the two tasks are achieved sequentially, the partitioning of the latent space is not optimal in the maximum likelihood sense. Conversely, the mixture of experts model (MoE) and its probabilistic view [START_REF] Xu | An alternative model for mixtures of experts[END_REF] allows to jointly optimize the partitioning of the latent space and local transformation parameters using EM. However, MoE is a generic model inspired by gating networks, and was not specifically designed for piecewise linear data.

In this paper we show that high-dimensional binaural observations, parameterized by low-dimensional motor-state parameters, enjoy an intrinsic manifold structure, and hence are locally linear. With this property in mind, we propose a novel SSL method based on a probabilistic piecewise affine regression model (PPAR) specifically designed to deal with high-dimensional acoustic data that have an intrinsic lowdimensional manifold structure. PPAR is based on a generative model that is used to learn the mapping between a set of motor states and associated binaural cues. We show that this mapping can be inverted in closed-form to obtain the full posterior density function p(X|y) in the latent space given a new y. We further generalize this inversion to the case of missing and/or redundant observations in order to solve for two-dimensional SSL of natural, sparse sounds. The proposed method uses an audio-motor learning stage with white noise and does not require any a priori knowledge on the system's geometry and parameters.

The remainder of this paper is organized as follows. Section 2 describes the binaural model and the audio-motor paradigm. Section 3 validates the manifoldness of the binaural data. Section 4 describes the PPAR method and its associated Bayesian inverse mapping extended to missing and redundant observations. Section 5, compares PPAR with MPLR [START_REF] Qiao | Mixture of probabilistic linear regressions: A unified view of GMM-based mapping techiques[END_REF], and shows SSL results obtained with PPAR. Section 6 concludes and provides directions for future work.

THE AUDIO-MOTOR MODEL AND DATA

Any SSL method needs a content-independent sound representation containing as much spatial information as possible. Rather than single ILD and ITD coefficients, we use frequency-dependent ILD and IPD values calculated from a short-term Fourier transform (STFT) analysis, as done in e.g. [START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF][START_REF] Deleforge | A latently constrained mixture model for audio source separation and localization[END_REF]. First, the complex-valued spectrograms associated with the two microphones are computed with a 64ms time-window and 8ms overlap, yielding T = 126 frames for a 1s signal. Since sounds are recorded at 16,000Hz, each time window contains 1,024 samples which are transformed into F = 512 complex Fourier coefficients associated to frequency channels between 0 and 8,000Hz. For a binaural recording made in the presence of a single sound source located at x in a listener-centered coordinate frame, we denote with {s (S)

f t } F,T
f,t=1 the complex-valued spectrogram emitted by the source, and with {s

(L) f t } F,T f,t=1 and {s (R) f t } F,T f,t=1
the left (and right) perceived spectrograms. The HRTF model provides a relationships between the emitted and the perceived spectrograms points:

s (L) f t = h (L) f (x) s (S) f t and s (R) f t = h (R) f (x) s (S) f t (1)
where h (L) and h (R) denote the left and right non-linear HRTFs. The interaural transfer function (ITF) is defined by the ratio between the two HRTFs, i.e.,

I f (x) = h (R) f /h (L) f ∈ C. The interaural spectrogram is defined by Îft = s (R) f t /s (L) f t , so that Îft ≈ I f (x).
This way, Îft does not depend on the emitted spectrogram value s (S) f t but only on the emitting source position x. However, this approximation holds only if the source is emitting at (f, t) (i.e. s (S) f t = 0). Therefore, interaural spectrograms of natural sounds are sparse. Missing interaural spectrograms values will be characterized by χ f t = 0 and χ f t = 1 otherwise. They can be determined using a threshold on left and right spectral powers |s

(L) f t | 2 and |s (R) f t | 2 .
We define the ILD spectrogram α and the IPD spectrogram φ as the log-amplitude and phase of the interaural spectrogram Îf,t :

α f t = 20 log | Îft | ∈ R, φ f t = exp(j arg( Îft )) ∈ C (2)
The phase difference is expressed in the complex domain, or equivalently R 2 , to avoid problems due to phase circularity. In the particular case of a sound source emitting white noise from x, we have χ t = 1 for all t, i.e., the sound source is emitting at all (f, t) points. One can thus compute the temporal means ᾱ(x) ∈ R F and φ(x) ∈ R 2F of ILD and IPD spectrograms. These mean vectors will be referred to as the mean ILD and the mean IPD vectors associated to x. The well established duplex theory suggests that ILD cues are mostly used at high frequencies (above 1.5 kHz) while ITD (or IPD) cues are mostly used at low frequencies (below 1.5kHz) in humans. Indeed, ILD values are similar because the HRTF can be neglected at low frequencies, and the phase difference becomes very unstable with respect to the source position at high frequencies. To account for these phenomena, the initial binaural cues are split into two distinct vectors, namely the mean low-ILD and high-ILD and the mean low-IPD and high-IPD vectors, where low corresponds to 96 frequency channels between 0 and 1.5kHz and high corresponds to 416 frequency channels between 1.5kHz and 8kHz.

To automatically gather a large number of such vectors associated with different source positions, we used the same technique and robot setup as in [START_REF] Deleforge | The cocktail party robot: Sound source separation and localisation with an active binaural head[END_REF]. A binaural acoustic dummy-head is mounted onto a pan-tilt (ψ, θ) motor system. The emitter (a loud speaker) is fixed in a reference position at 2.5 meters in front of the robot, while the robot is placed in 160 pan angles ψ ∈ [-160 i.e., N = 9, 600 uniformly distributed motor states. Hence, the source location spans a 320 • azimuth range and a 120 • elevation range in the robot's frame. Notice that there is a one-to-one association between motor states and source locations. They will be denoted by {x n } N n=1 . For each x n ∈ R 2 , two binaural recordings are made: (i) white noise which can be used to estimate ᾱ(x n ) and φ(x n ), and (ii) a randomly picked utterance amongst 362 samples from the TIMIT dataset [START_REF] Garofolo | The DARPA TIMIT acoustic-phonetic continuous speech corpus[END_REF]. These are 50% female, 50% male and they last 1-5s. Sparse interaural spectrograms can be computed from these recordings and are used to test our sound localization algorithm on natural sounds in section 5. All the experiments were carried out in real-world conditions, i.e., a room with natural reverberations and background noise.

THE MANIFOLD OF INTERAURAL CUES

In this section, we analyze the intrinsic structure of the mean high-and low-ILD and IPD vectors previously described. While these vectors live in a high-dimensional space, they should be parameterized by motor states and hence, lie on a lower L-dimensional manifold with L = 2. We propose to experimentally verify the existence of this manifold structure using non-linear dimensionality reduction, and examine whether obtained representations are homeomorphic to the motor-state space. Such a homeomorphism would allow us to confirm (or invalidate) the existence of a locally linear bijective mapping between motor states (or equivalently, source positions) and the interaural data gathered with our setup.

If the interaural data lie in a linear subspace, a standard dimensionality reduction method such as PCA could be used. However, in the case of a non-linear subspace one should use a manifold learning technique, e.g., diffusion kernels [START_REF] Talmon | Supervised source localization using diffusion kernels[END_REF]. Alternatively, we chose to use local tangent-space alignment (LTSA) [START_REF] Zhang | Principal manifolds and nonlinear dimensionality reduction via tangent space alignment[END_REF] because it essentially relies on the assumption that the data are locally linear, which is our central hypothesis. LTSA starts by building a local neighborhood around each high-dimensional observation. Under the key assumption that each such neighborhood spans a linear space of low dimension corresponding to the dimensionality of the tangent space, i.e., a Riemannian manifold, PCA can be applied to each one of these neighborhoods thus yielding as many lowdimensional data representations as points in the data set. Finally a global low-dimensional map is built by optimal alignment of these local representations (see [START_REF] Zhang | Principal manifolds and nonlinear dimensionality reduction via tangent space alignment[END_REF] for details). Twodimensional1 maps obtained using LTSA are shown in Fig. 1. Mean low-ILD, low-IPD, and high-ILD maps are smooth and homeomorphic to the motor-state space, thus confirming that these cues can be used for 2D binaural SSL based on locally linear mapping. However, this is not the case for the mean high-IPD map which features several distortions, elevation ambiguities, and crossings. While these computational experiments confirm the duplex theory for IPD cues, they surprisingly suggest that ILD cues at low frequencies still contain rich enough 2D sound-source position information. We can therefore concatenate full-spectrum ILD and low-frequency IPD vectors to form an observation space in R D (with D = 704), which will be referred to as the ILPD space. Similarly, we can define sparse ILPD spectrograms for general sounds. Fig. 2 shows the result of applying PCA globally to the mean high-ILD data with L = 2. One may observe that the resulting map is extremely distorted, due to the nonlinear nature of the high-ILD manifold. This rules out the use of a linear regression method to estimate the interauralto-localization mapping and justifies the development of an appropriate piecewise-linear mapping method.

PROBABILISTIC PIECEWISE AFFINE MAPPING

Let X ⊂ R L be a subset of the low-dimensional latent space, i.e., the space of sound-source positions (or motor states), and R D be the high-dimensional observation space, i.e., the space of ILPD cues. The computational experiments of section 3 suggest that there exists a smooth, locally linear bijection g : X → Y ⊂ R D such that the set Y = {g(x), x ∈ X } forms an L-dimensional manifold embedded in R D . Based on this assumption, the proposed idea is to compute a piecewiseaffine probabilistic approximation of g from a training data set {(x n , y n )} N n=1 ⊂ X × Y and to estimate the inverse of g using a Bayesian formulation. The local linearity of g suggests that each point y n is the image of a point x n ∈ R k ⊂ X by an affine transformation t k , plus an error term. Assuming that there is a finite number K of such affine transformations t k and an equal number of associated regions R k we obtain a piecewise-affine approximation of g. With each training-couple (x n , y n ) we associate an assignment variable z n = (z 1n . . . z Kn ) ⊤ such that z kn = 1 if y n is the image of x n ∈ R k by t k and 0 otherwise. This allows us to write:

y n = K k=1 z kn (A k x n + b k ) + e n (3)
where the D × L matrix A k and the vector b k ∈ R D define the transformation t k , and e n ∈ R D is an error term capturing both the observation noise and the reconstruction error of affine transformations. If we make the assumption that the error terms e n do not depend on x n , y n or z n , and are iid realizations of a Gaussian variable with 0 mean and diagonal covariance matrix Σ = diag(σ 2 1:D ) we obtain:

p(y n |x n , z kn = 1; Θ) = N (y n ; A k x n + b k , Σ) (4)
where Θ designates all the model parameters [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF]. To make the affine transformations local, we model z n by the realization of a hidden multinomial random variable conditioned by x n :

p(z kn = 1|x n ; Θ) = π k N (x n ; c k , Γ k ) K k=1 π k N (x n ; c k , Γ k ) , (5)
where

c k ∈ R L , Γ k ∈ R L×L and k π k = 1.
We can give a geometrical interpretation of this distribution by adding the following volume equality constraints to the model:

|Γ 1 | = • • • = |Γ K | and π 1 = • • • = π K = 1/K (6)
One can verify that under these constraints, the set of K regions of X maximizing (5) for each k defines a Voronoi diagram of centroids {c k } K k=1 , where the Mahalanobis distance ||.|| Γ k is used instead of the Euclidean one. This corresponds to a compact probabilistic way of representing a general partitioning of the latent space into convex regions of equal volume. Extensive tests on simulated and audio data showed that these constraints yielded lower reconstruction errors, on top of providing a meaningful interpretation of (5). To make our generative model complete, we define the following Gaussian mixture prior on the latent variables:

p(x n ; Θ) = K k=1 π k N (x n ; c k , Γ k ) (7)
This model yields a closed-form and efficient EM algorithm maximizing the observed-data log-likelihood log p(X, Y; Θ) with respect to the model's parameters:

Θ = {Γ k , c k , A k , b k , π k } K k=1 , Σ (8)
Posterior probabilities r

(i) kn = p(z kn = 1|x n , y n ; Θ (i-1)
) are computed in the E-step from ( 4), ( 5) and Bayes inversion. The M-step maximizes the expected complete-data log-likelihood

E (Z|X ,Y ,Θ (i) ) [log p(X, Y , Z|Θ)].
We obtain the following closed-form expressions for the parameters updates under the volume equality constraints (6):

c (i) k = N n=1 r (i) kn r(i) k x n , Γ (i) k = K k=1 r(i) k N |S (i) k | 1 L S (i) k |S (i) k | 1 L (9) A (i) k = Y (i) k X (i) † k , b (i) k = N n=1 r (i) kn r(i) k (y n -A (i)
k x n ), (10)
σ 2 (i) d = 1 K K k=1 N n=1 r (i) kn r(i) k (y dn -a (i)⊤ dk x n -b (i) dk ) 2 , ( 11
)
where † is the Moore-Penrose pseudo inverse operator and:

S (i) k = N n=1 r (i) kn /r (i) k (x n -c (i) k )(x n -c (i) k ) ⊤ r(i) k = K k=1 r (i) kn , A (i) k = (a (i) 1k . . . a (i) Dk ) ⊤ , a (i) dk ∈ R L X (i) k = (r (i) 1 2 k1 (x 1 - x(i) k ) . . . r (i) 1 2 kN (x N - x(i) k )) Y (i) k = (r (i) 1 2 k1 (y 1 - ȳ(i) k ) . . . r (i) 1 2 kN (y N - ȳ(i) k )) x(i) k = N n=1 r (i) kn /r (i) k x n and ȳ(i) k = N n=1 r (i) kn /r (i) k y n Initial posteriors r (0)
kn can be obtained either by estimating a K-GMM solely on X (GMM) or jointly on (X, Y ) (GMM-J). One may see that the MPLR-GMM and MPLR-GMM-J methods proposed in [START_REF] Qiao | Mixture of probabilistic linear regressions: A unified view of GMM-based mapping techiques[END_REF] are strictly equivalent to these initializations followed by a single parameter estimation using [START_REF] Kain | Spectral voice conversion for textto-speech synthesis[END_REF]. Unlike [START_REF] Qiao | Mixture of probabilistic linear regressions: A unified view of GMM-based mapping techiques[END_REF], we use a consistent probabilistic model enforcing partitioning within parameter estimation, a component-dependent reconstruction error Σ, and an Estep allowing to iterate until convergence to a maximum of the log-likelihood. This significantly reduces the global reconstruction error (section 5). Fig. 3 shows a partitioning example obtained with our method using a toy data set.

Let Θ denote the parameters estimated with our algorithm; we describe now a method based on Bayesian inversion to estimate the unknown position x of a sound t,d=1 as defined in section 2 and 3. If we suppose that the observations are assigned to the same position x and transformation z, it follows from our model ( 4), ( 5), ( 7) that the posterior distribution p(x|Y χ ; Θ) is a GMM

K k=1 α k N (x; µ k , V k ) in R L with parameters: µ k = V k Γ -1 k c k + D,T d,t=1 χ dt σ 2 d a dk (y dt -b dk ) , (12)
V k = Γ -1 k + D,T d,t=1 χ dt σ 2 d a dk a ⊤ dk -1 and (13)
α k ∝ (|V k |/| Γ k |) 1 2 exp -1 2 ( D,T d,t=1 χ dt e σ 2 d (y dt -b dk ) 2 + c ⊤ k Γ -1 k c k + µ ⊤ k V -1 k µ k ) (14)
where the weights {α k } K k=1 are normalized to sum to 1. This formulation is more general than the unique, complete observation case (T = 1, χ 1 = 1). Several estimates can be inferred from this posterior distribution. We used the posterior expectation E[x|Y χ ] = K k=1 α k µ k which yields the lowest average localization errors. However, the posterior distribution may have several high modes, in which case the posterior expectation becomes a misleading estimator. This notably happens if the sound spectrum is extremely sparse. In this case, one may preserve the full posterior distribution and, for instance, combine it with other external probabilistic knowledge.

EXPERIMENTS AND RESULTS

We compared the proposed algorithm with two different initialization strategies, PPAR-GMM and PPAR-GMM-J, to the two algorithms proposed in [START_REF] Qiao | Mixture of probabilistic linear regressions: A unified view of GMM-based mapping techiques[END_REF], MPLR-GMM and MPLR-GMM-J. We used these four algorithms for training with the complete set of N = 9, 600 ILPD white-noise samples available from the audio-motor recordings. The algorithm performance was evaluated using the mean reconstruction error, , where σ2 d is defined by [START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF]. As shown in Fig. 4, the proposed method significantly outperforms the two others at a minor additional computational cost. This is not surprising since PPAR yields maximum likelihood estimators of the affine transformations and space partitioning obtained after several EM iterations 2 , while MPLR-GMM and MPLR-GMM-J perform just one optimization step [START_REF] Kain | Spectral voice conversion for textto-speech synthesis[END_REF]. GMM-J initialization strategy will be used from now on since it provides better results than GMM initialization. Interestingly, the only free parameter of our method, namely K, can be chosen based on a compromise between computational cost and precision, since the higher the value of K the lower the error. Next we evaluate the proposed inverse mapping method in the case of a complete, unique observation, i.e., sound source localization using a mean ILPD vectors computed from white noise. Localization error is evaluated as a function of the density of the training set. We represent the density by the angle δ, corresponding to the average pan and tilt absolute difference between a point and its nearest neighbor in the training set. The complete training set of density 2 • was uniformly decimated at random to obtain irregularly spaced and smaller training sets, while the test positions were randomly chosen from the complete set so that most of the test positions were outside the training set. For a given density δ, 10 different decimated sets were used for training, and 20 source positions were estimated for each one, i.e., 200 localization tasks. K was chosen such that there are approximately 30 training samples per affine transformation. The mean and standard deviation of the errors in azimuth and in elevation are shown in Fig. 5. The mean localization errors are always smaller than half the training set density, which illustrates the interpolation ability of our method. In addition, the error's standard deviation remains reasonable even for heavily decimated training sets, thus showing that the overall performance is not much affected by the distribution and size of the training set being used. No front-back azimuth or elevation confusions were observed, thanks to the asymmetry of the dummy head and to the spatial richness of ILPD cues. This is in contrast with most of the current binaural SSL approaches that focus only on frontal azimuth estimation.

Finally, we tested our SSL method with missing and redundant observations using randomly located sound sources emitting random utterances. The spectral power threshold was manually set quite high, so that the test ILPD spectrograms had 89.6% of missing data in average. Mean azimuth and elevation errors (Avg), standard deviations (Std), and percentage higher than 4×Avg (Out) over 500 tests are shown in the table below. The low Stds and low amount of high errors show the reliability and robustness of our method, even using real world recordings of very sparse sounds emitted from a wide range of azimuths and elevations.

CONCLUSION

In this paper, we emphasized the manifold structure of realworld audio-motor data using ILD/IPD cues, and we set the basis of a novel probabilistic framework for accurate 2D sound source localization relying on this structure. More generally, we presented a new methodology for examining whether a dataset is intrinsically locally linear based on manifold learning, and proposed an adequate probabilistic mapping method relying on this property. The advantage of our Bayesian formulation is that it could serve as a basis for a wide variety of extensions and other applications. For example, a mixture of PPAR could be used for multiple sound sources localization. In the future, we plan to thoroughly examine how changes in the recording environment affect the acoustic manifold structure. Such a study and a proper adaptation of our mapping model would allow to make our sound source localization method robust to changes in reverberations and positions in the room.

Fig. 1 :Fig. 2 :

 12 Fig.1: Two-dimensional representations of mean interaural vectors using LTSA. For visualization purposes, points corresponding to the same tilt angle (elevation) have the same color and they are linked in pan angle (azimuth) order with a black line.

Fig. 3 :

 3 Fig. 3: Latent space partitioning and locally affine mapping on a toy data set (N = 9600, K = 15, L = 2, D = 3). Colors encode regions in R D maximizing[START_REF] Mandel | An EM algorithm for localizing multiple sound sources in reverberant environments[END_REF]. Observe how well these regions (associated with affine transformations) are adjusted to the geometry of both the latent space and the observed manifold.

Fig. 4 :

 4 Fig. 4: Comparison of the proposed method (PPAR) with MPLR [14] as a function of the number of components. namely the ratio mean d [σ 2 d ]/var n [y dn ], where σ 2 d is defined by[START_REF] Stylianou | Continuous probabilistic transform for voice conversion[END_REF]. As shown in Fig.4, the proposed method significantly outperforms the two others at a minor additional computational cost. This is not surprising since PPAR yields maximum likelihood estimators of the affine transformations and space partitioning obtained after several EM iterations 2 , while MPLR-GMM and MPLR-GMM-J perform just one optimization step[START_REF] Kain | Spectral voice conversion for textto-speech synthesis[END_REF]. GMM-J initialization strategy will be used from now on since it provides better results than GMM initialization. Interestingly, the only free parameter of our method, namely K, can be chosen based on a compromise between computational cost and precision, since the higher the value of K the lower the error. Next we evaluate the proposed inverse

Fig. 5 :

 5 Fig. 5: Average and standard deviation of angular errors in white noise localization as a function of the training set's density.

  • 1.8 • 1.2% 1.3 • 1.2 • 1.2% δ = 6 • 5.6 • 5.2 • 1.2% 4.5 • 5.3 • 1.4%

	Training	Azimuth		Elevation
	density	Avg Std	Out	Avg Std	Out
	δ = 2 •	2.0

When applying PCA locally, we observed a significant eigengap between the

2nd and 3rd eigeinvalues, thus validating the choice L = 2.

Convergence was reached after 10 to 20 iterations using mean ILPD data.


Antoine Deleforge

Radu Horaud

A Latently Constrained Mixture Model for Audio Source Separation and Localization

We present a method for audio source separation and localization from binaural recordings. The method combines a new generative probabilistic model with time-frequency masking. We suggest that device-dependent relationships between point-source positions and interaural spectral cues may be learnt in order to constrain a mixture model. This allows to capture subtle separation and localization features embedded in the auditory data. We illustrate our method with data composed of two and three mixed speech signals in the presence of reverberations. Using standard evaluation metrics, we compare our method with a recent binaural-based source separation-localization algorithm.

Introduction

We address the problem of simultaneous separation and localization of sound sources mixed in an acoustical environment and recorded with two microphones. Time-frequency masking is a technique allowing the separation of an arbitrary number of sources with only two microphones by assuming that a single source is active at every time-frequency point -the W-disjoint orthogonality (W-DO). It was shown that this assumption holds, in general, for simultaneous speech signals [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. The input signal is represented in a time-frequency domain and points corresponding to the target source are weighted with 1 and otherwise with 0. The masked spectrogram is then converted back to a temporal signal. A number of methods combine time-frequency masking with localization-based clustering ( [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], [START_REF] Mouba | A source localization/separation/respatialization system based on unsupervised classification of interaural cues[END_REF], [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF]), e.g., DUET [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] which allows to separate anechoic mixtures when each source reaches the microphones with a single attenuation coefficient and delay. This mixing model is well suited for "clean" binaural recordings. In practice, more complex filtering effects exist, namely the head-related transfer function (HRTF) and the room impulse response (RIR). These filters lead to frequencydependent attenuations and delays between the two microphones, respectively called the interaural level difference (ILD) and the interaural phase difference (IPD). Some approaches attempted to account for these dependencies by learning a mapping between azimuth, frequencies and interaural cues [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF][START_REF] Roman | Speech segregation based on sound localization[END_REF][START_REF] Mouba | A source localization/separation/respatialization system based on unsupervised classification of interaural cues[END_REF]. These mappings usually consist in finding a functional relationship that best fits data obtained from an HRTF dataset. To improve robustness to RIR variations, these interaural cues can also be integrated in a mixture model, e.g., [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF].

In this paper we propose to directly learn a discrete mapping between a set of 3D point sources and IPD/ILD spectral cues. We will refer to such mappings as Source-Position-to-Interaural-Cues maps (SPIC maps). Unlike what is done in [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF][START_REF] Roman | Speech segregation based on sound localization[END_REF][START_REF] Mouba | A source localization/separation/respatialization system based on unsupervised classification of interaural cues[END_REF], the proposed mapping is built point-wise, does not rely on azimuth only and is device-dependent. We explicitly incorporate it into a novel latently constrained mixture model for point sound sources. Our model is specifically designed to capture the richness of binaural data recorded with an acoustic dummy head, and this to improve both localization and separation performances. We formally derive an EM algorithm that iteratively performs separation (E-step) followed by localization and source-parameter estimation (M-step). The algorithm is supervised by a training stage consisting in learning a mapping between potential source positions and interaural cues, i.e., SPIC maps. We believe that a number of methods could be used in practice to learn such maps. In particular we propose an audio-motor mapping approach. The results obtained with our method compare favorably with the recently proposed MESSL algorithm [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF].

Binaural Sound Representation

Spectrograms associated with each one of the two microphones are computed using short-term FFT analysis. We use a 64ms time-window with 8ms window overlap, thus yielding T = 126 time windows for a 1s signal. Since sounds were recorded at a sample rate of 16,000Hz, each time window contains 1,024 samples. Each window is then transformed via FFT to obtain complex coefficients of F = 513 positive frequency channels between 0 and 8,000Hz. We denote with s f,t the spectrogram points perceived by the left-and right-microphone respectively. The W-DO assumption implies that a single sound source k emits at a given point (f, t). The relationships between the emitted and the left and right perceived spectrogram points are:

s (L) f,t = h (L) (x k , f ) s (k) f,t and s (R) f,t = h (R) (x k , f ) s (k) f,t (1)
where x k ∈ R 3 is the 3D position of sound source k in a listener-centered coordinate frame and h (L) and h (R) denote the left and right HRTFs. The interaural transfer function (ITF) is defined by the ratio between the two HRTFs, i.e.,

I(x k , f ) = h (R) (x k , f )/h (L) (x k , f ) ∈ C. The interaural spectrogram is defined by Îf,t := s (R) f,t /s (L)
f,t , so that Îf,t ≈ I(x k , f ). Note that the last approximation only holds if there is a source k emitting at frequency-time point (f, t), and if the time delay between microphones (≈ 0.75ms) is much smaller than the Fourier transform time-window that is used (64ms). Under these conditions, at a given frequency-time point, the interaural spectrogram value Îf,t does not depend on the emitted spectrogram value s (k) f,t but only on the emitting source position x k . We finally define the ILD spectrogram α and the IPD spectrogram φ as the log-amplitude and phase of the complex interaural spectrogram Îf,t :

α f,t = 20 log | Îf,t | ∈ R, φ f,t = arg( Îf,t ) ∈ ]-π, π] (2)
As already outlined in Section 1 our method makes use of a SPIC map that is learnt during a training stage. Let X = {x n } N n=1 be a set of 3D sound-source locations in a listener-centered coordinate frame. Let a sound-source n, located at x n emit white noise and let {α n f,t } F,T f =1,t=1 and {φ n f,t } F,T f =1,t=1 be the perceived ILD and IPD spectrograms. The mean ILD µ(

x n ) = (µ n 1 . . . µ n f . . . µ n F ) ⊤ ∈ R F and mean IPD ξ(x n ) = (ξ n 1 . . . ξ n f . . . ξ n F ) ⊤ ∈ ]-π, π]
vectors associated with n are defined by taking the temporal means of α n and φ n at each frequency channel:

µ n f = 1/T T t=1 α n f,t and ξ n f = arg(1/T T t=1 e jφ n f,t ) (3)
Vector ξ is estimated in the complex domain in order to avoid problems due to phase circularity [START_REF] Nix | Sound source localization in real sound fields based on empirical statistics of interaural parameters[END_REF]. White noise is used because it contains equal power within a fixed bandwidth at any center frequency: The source n is therefore the only source emitting at each point (f, t); µ n f and ξ n f are thus approximating the logamplitude and phase of I(x k , f ). The set X of 3D source locations as well as the mappings µ and ξ will be referred to as the training data to be used in conjunction with the separation-localization algorithm described below.

Constrained Mixtures for Separation and Localization

Let's suppose now that there are K simultaneously emitting sounds sources from unknown locations {x k } K k=1 ⊂ X and with unknown spectrograms. Using the listener's microphone pair it is possible to build the ILD and IPD observed spectrograms {α f,t } F,T f =1,t=1 and {φ f,t } F,T f =1,t=1 . The goal of the sound-source separation and localization algorithm described in this section is to associate each observed point (f, t) with a single source and to estimate the 3D location of each source.

As mentioned in section 2, the observations α f,t (ILD) and φ f,t (IPD) are significant only if there is a sound source emitting at (f, t). To identify such significant observations we estimate the sound intensity level (SIL) spectrogram at the two microphones, and retain only those frequency-time points for which the SIL is above some threshold. One empirical way to choose the thresholds (one for each frequency) is to average the SILs at each f in the absence of any emitting source. These thresholds are typically very low compared to SILs of natural sounds, and allow to filter out frequency-time points corresponding to "room silence". Let M f ≤ T be the number of significant observations at f and let α f,m and φ f,m be the m-th significant ILD and IPD observations at f . Let

A = {α f,m } F,M f f =1,m=1 and Φ = {φ f,m } F,M f f =1,m=1 be the observed data.
Let z f,m ∈ {0, 1} K be the missing data, i.e., the data-to-source assignment variables, such that z f,m,k = 1 if observations α f,m and φ f,m are generated by source k, and z f,m,k = 0 otherwise. The W-DO assumption yields

K k=1 z f,m,k = 1 for all (f, m). M k = {z f,m,k } F,M f f =1,m=1 is the binary spectral mask of the k-th source. Finally, Z = {z f,m } F,M f f =1,m=1
denotes the set of all missing data. The problem of simultaneous localization and separation amounts to estimate the masking variables Z and the locations {x k } K k=1 conditioned by A and Φ, given the number of sources K. We assume that observed data are perturbed by Gaussian noise. Hence, the probability of observing α f,m conditioned by source k (z f,m,k = 1) located at x k is drawn from a normal distribution, and the probability of observing φ f,m is drawn from a circular normal distribution. The source position x k acts here as a latent constraint on ILD and IPD means:

P (α f,m |z f,m,k = 1, x k , σ f,k ) = N (α f,m |µ f (x k ), σ 2 f,k
) and ( 4)

P (φ f,m |z f,m,k = 1, x k , ρ f,k ) = N (∆(φ f,m , ξ f (x k ))|0, ρ 2 f,k ) (5)
where σ 2 f,k and ρ 2 f,k are the ILD and IPD variances associated with source k at frequency f and the ∆ function is defined by ∆(x, y) = arg(e j(x-y) ) ∈ ]-π, π]. As in [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF], (5) approximates the normal distribution on the circle ]-π, π] when ρ f,k is small relative to 2π. Preliminary experiments on IPD spectrograms of white noise showed that this assumption holds in the general case. As emphasized in [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF], the well known correlation between ILD and IPD does not contradict the assumption that Gaussian noises corrupting the observations are independent. The conditional likelihood of the observed data (α f,m , φ f,m ) is therefore given by the product of ( 4) and [START_REF] Roman | Speech segregation based on sound localization[END_REF]. We also define the priors π f,k = P (z f,m,k ) which model the proportion of the observed data generated by source k at frequency f . In summary, the model parameters are

Θ = {{x k }; {π f,k }; {σ 2 f,k }; {ρ 2 f,k }} F,K f =1,k=1
. The problem can now be expressed as the maximization of the observed-data log-likelihood conditioned by Θ. In order to keep the model as general as possible, there is no assumption on the emitted sounds as well as the way their spectra are spread across the frequency-time points. Therefore, we assume that all the observations are statistically independent, yielding the following expression for the observed-data log-likelihood:

L(A, Φ; Θ) = log P (A, Φ; Θ) = F f =1 M f m=1 log P (α f,m , φ f,m ; Θ) (6)
We address this maximum-likelihood with missing-data problem within the framework of expectation-maximization (EM). In our case, the E-step computes the posterior probabilities of assigning each spectrogram point to a sound source k (separation) while the M-step maximizes the expected complete-data log-likelihood with respect to the model parameters Θ and, most notably, with the source locations {x k } K k=1 (localization). The MAP criterion provides binary spectral masks M k associated with each source k while the final parameters {x k } K k=1 provide estimates for the source locations. The expected complete-data log-likelihood writes ( (p) denotes the p-th iteration):

Q(Θ|Θ (p-1) ) = F f =1 M f m=1 K k=1 r (p) f,m,k log π f,k P (α f,m , φ f,m |z f,m ; Θ) (7)
The E-step updates the responsibilities according to the standard formula:

r (p) f,m,k = π f,k P (α f,m , φ f,m |z f,m ; Θ (p-1) ) K i=1 π f,i P (α f,m , φ f,m |z f,m ; Θ (p-1) ) (8)
The M-step maximizes [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF] with respect to Θ. By combining (4) and ( 5) with [START_REF] Viste | On the use of spatial cues to improve binaural source separation[END_REF] the equivalent minimization criterion writes:

F f =1 M f m=1 r (p) f,m,k log σ 2 f,k ρ 2 f,k π 2 f,k + (x f,m -µ f (x k )) 2 σ 2 f,k + ∆(φ f,m , ξ f (x k )) 2 ρ 2 f,k (9)
which can be differentiated with respect to {π f,k } f , {σ f,k } f and {ρ f,k } f to obtain closed-form expressions for the optimal parameter values conditioned by x k :

π f,k = r f,k M f , with r f,k = M f m=1 r f,m,k (10)
σ 2 f,k (x k ) = 1 r f,k M f m=1 r (p) f,m,k (x f,m -µ f (x k )) 2 (11) ρ 2 f,k (x k ) = 1 r f,k M f m=1 r (p) f,m,k ∆(φ f,m , ξ f (x k )) 2 (12)
By substituting (11) and ( 12) into (9) the optimal location x k is obtained by minimizing the following expression with respect to x k :

F f =1 r f,k log 1 + (α f,k -µ f (x k )) 2 V f,k + log 1 + ∆(φ f,k , ξ f (x k )) 2 W f,k (13)
with:

α f,k = 1 r f,k M f m=1 r (p) f,m,k α f,m ; V f,k = 1 r f,k M f m=1 r (p) f,m,k (α f,m -α f,k ) 2 φ f,k = arg 1 r f,k M f m=1 r (p) f,m,k e jφ f,m ; W f,k = 1 r f,k M f m=1 r (p) f,m,k ∆(φ f,m , φ f,k ) 2
(13) is evaluated for each source location in the training dataset X (Section 2) in order to find an optimal 3D location x k . This is then substituted back in (11) and ( 12) to estimate σ f,k and ρ f,k and repeated for each unknown source k.

In general, EM converges to a local maximum of [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. The non-injectivity nature of the interaural functions µ f and ξ f and the high cardinality of Θ leads to a very large number of such maxima, especially when the training set X is large. This makes our algorithm very sensitive to initialization. One way to avoid being trapped in local maxima is to initialize the mixture's parameters at random several times. This cannot be easily applied here since there is no straightforward way to initialize the model's variances. Alternatively, one may randomly initialize the assignment variables Z and then proceed with the Mstep. However, extensive simulated experiments revealed that this solution fails to converge to the ground-truth solution in most of the cases. We therefore propose to combine these strategies by randomly perturbing both the source locations and the source assignments during the first stages of the algorithm. We developed a stochastic initialization procedure similar in spirit to SEM [START_REF] Celeux | A classification EM algorithm for clustering and two stochastic versions[END_REF].

The SEM algorithm includes a stochastic step (S) between the E-and the M-step, during which random samples R f,m,k ∈ {0, 1} are drawn from the responsibilities [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. These samples are then used instead of ( 8) during the M-step. To initialize our algorithm, we first set r (0) f,m,k = 1/K for all k and then proceed through the sequence S M* E S M, where M* is a variation of M in which the source positions are drawn randomly from X instead of solving (13). In practice, ten such initializations are used to enforce algorithm convergence, and only the one providing the best log-likelihood after two iterations is iterated twenty more times. A second technique was used to overcome local maxima issues due to the large number of parameters. During the first ten steps of the algorithm only, a unique pair of variances (σ 2 k , ρ 2 k ) is estimated for each source. This is done by calculating the means σ 2 k (x k ) and ρ 2 k (x k ) of frequency-dependent variances (11) and (12) weighted by r f,k . The optimal value x k is the one minimizing σ 2 k (x)ρ 2 k (x) evaluated over all x ∈ X . Intensive experiments showed that the proposed method converges to a global optimum in most of the cases.

Experiments, Results, and Conclusions

In order to evaluate and compare our method, a specific data set of binaural records was built1 using a Sennheiser MKE 2002 acoustic dummy-head mounted onto a robotic system with two rotational degrees of freedom, namely pan (ψ) and tilt (θ). This device, specifically designed to perform accurate and reproducible motions, allows us to collect both a very dense SPIC map for the training set (section 2) and a large test set of mixed speech point sources. The emitter (a loud speaker) is placed at approximately 2.5 meters in front of the listener. Under these conditions the HRTF mainly depends on the sound-source direction: Hence, the location is parameterized by the angles ψ and θ. All the experiments were carried out in a reverberant room and in the presence of background noise. For recording purposes, the robot is placed in 90 pan angles ψ ∈ [-90 • , 90 • ] (leftright) and 60 tilt angles θ ∈ [-60 • , 60 • ] (top-down), i.e., N = 5, 400 uniformly distributed motor states in front of the static emitter, forming the set X . Five binaural recordings are available with each motor state: Sound #0 corresponds to a 1s "room silence" used to estimate the SIL thresholds (section 3). Sound #1 corresponds to 1s white-noise used to build the training set (section 2). Sounds #2, #3 and #4 form the test set and correspond to "They never met you know" by a female (#2), "It was time to go up myself" by a male (#3), and "As we ate we talked" by a male (#4). The three sounds are about 2s long and were randomly chosen from the TIMIT database. Each record was associated to its ground-truth motor-state, thus allowing to create signals of mixed sound sources from different direction with 2 • resolution.

We generated 1000 mixtures of two and three speech signals emitted by randomly located sources. 97.7% of the individual sources were correctly mapped to 1. Comparing the mean source-to-distortion ratio (SDR) and source-to-interference ratio (SIR), in dB, for 1000 mixtures of 2 and 3 sources. Mean separation results with our approach are calculated over all sources (All) and over correctly localized sources only (Loc).

their associated position (i.e. ≤ 2 • error for both ψ and θ) in the two-source case, and 63.8% in the three-source case. The performance of separation was evaluated with the standard SDR and SIR metrics [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF]. We compared our results to those obtained with the original mixture (no mask applied), with the ground truth or Oracle mask [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF], and with the recently proposed MESSL2 algorithm [START_REF] Mandel | Model-based expectationmaximization source separation and localization[END_REF]. The Oracle mask is set to 1 at every spectrogram point in which the target signal is at least as loud as the combined other signals and 0 everywhere else. The version MESSL-G used includes a garbage component and ILD priors to better account for reverberations and is reported to outperform four methods in reverberant conditions, including [START_REF] Yılmaz | Blind separation of speech mixtures via time-frequency masking[END_REF] and [START_REF] Mouba | A source localization/separation/respatialization system based on unsupervised classification of interaural cues[END_REF]. Table 1 shows that our method yields significantly better results than MESSL-G on an average, although both algorithm require similar computational times. Notice how the localization correctness critically affects the separation performances, and decreases in the three-source case, as the number of observations per source becomes lower and the number of local maxima in (6) becomes higher. Our SDR scores strongly outperform MESSL-G in most cases, while SIR results are only slightly better when sources are more than 70 • apart in azimuth (pan angle), e.g., Fig. 1. However, they become much higher when sources are nearby in azimuth, or share the same azimuthal plane with different elevations (tilt angles). This is because MESSL relies on the estimation of a probability density in a discretized ITD space for each source, and thus does not account for more subtle spatial cues induced by the HRTF.

These results clearly demonstrate the efficiency of our method, but they somehow favor our algorithm because of the absence of RIR variations both in the training and the test data sets. The aim of experimenting with these relatively simple data has been to show that our method can conceptually separate and accurately locate both in azimuth and elevation a binaural mixture of 2 to 3 sound sources. The prerequisite is a training stage: the interaural cues associated with source positions need to be learnt in advance using white noise, and we showed that the algorithm performs well even for a very large and dense set of learnt positions. Preliminary results obtained while changing the position of the test sound source in the room suggested that our constrained mixture model coupled with frequency-dependent variances presented some robustness to RIR variations. Alternatively, one could build a training set on different premises such as seat locations in a conference room or musician locations in a concert hall, and thus directly learn the RIR during the training stage. To conclude, we proposed a novel audio source separation and localization method based on a mixture model constrained by a SPIC map. Experiments and comparisons showed that our algorithm performs better than a recently published probabilistic spectral masking technique in terms of separation and yields very good multi-source localization results. The combination of a SPIC map with a mixture model is a unique feature. In the future, we plan to study more thoroughly the behavior of our algorithm to RIR variations, and improve its robustness by extending our model to a continuous and probabilistic mapping between source positions and interaural parameters. Dynamic models incorporating moving sound sources and head movements could also be included.

  ∈ C the (f, t) point of the spectrogram emitted by sound-source k, and with s

Fig. 1 .

 1 Fig. 1. SIR as a function of azimuth (pan) and elevation (tilt) separation between two sources. Left: one source fixed at (-90 • , 0 • ) while the other takes 90 positions between (-90 • , 0 • ) and (+90 • , 0 • ). Right: one source fixed at (0 • , -60 • ) while the other takes 60 positions between (0 •, -60 • ) and (0 • , +60 • ). SIRs are averaged over 6 mixtures of 2 sources (12 targets). Topto-down: Oracle ( * ), our method (•), MESSL-G (△), and original mixture (+).

Table

		2 Sources 3 Sources
		SDR SIR SDR SIR
	Oracle Mask	11.73 19.23 9.20 16.16
	Our Approach (Loc) 5.28 8.91 2.44 3.92
	Our Approach (All) 5.19 8.84 1.72 2.74
	MESSL-G	2.83 5.74 1.48 1.47
	Original Mixture 0.00 0.45 -3.50 -2.82

Online at: http://perception.inrialpes.fr/ ~Deleforge/CAMIL_Dataset

http://blog.mr-pc.org/2011/09/14/messl-code-online/.


Michael Sapienza

Miles Hansard
email: miles.hansard@eecs.qmul.ac.uk

Radu Horaud
email: radu.horaud@inrialpes.fr

Real-time visuomotor update of an active binocular head

Keywords: Real-time vision, Active binocular vision, Visual tracking, 3D reconstruction

In order for a binocular head to perform optimal 3D tracking, it should be able to verge its cameras actively, while maintaining geometric calibration. In this work we introduce a calibration update procedure, which allows a robotic head to simultaneously fixate, track, and reconstruct a moving object in real-time. The update method is based on a mapping from motorbased to image-based estimates of the camera orientations, estimated in an offline stage. Following this, a fast online procedure is presented to update the calibration of an active binocular camera pair. The proposed approach is ideal for active vision applications because no image-processing is needed at runtime for the scope of calibrating the system or for maintaining the calibration parameters during camera vergence. We show that this homography-based technique allows an active binocular robot to fixate and track an object, whilst performing 3D reconstruction concurrently in real-time.

Introduction

The estimation of scene-structure from a binocular image pair is an important task in robot vision. Once the epipolar geometry between the two cameras is known, image-features can be matched more easily, and depth information can be recovered [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. For an active binocular robot head, which continuously fixates a moving target, a real-time method for updating the relationship between the cameras is required. Such an online method would allow independent or coupled camera vergence rotations while maintaining geometric calibration. This will enable the robot to perform 3D reconstruction while actively examining different parts of the scene, or while tracking a moving object [START_REF] Grosso | Active/dynamic stereo vision[END_REF][START_REF] Barreto | Active stereo tracking of N ≤ 3 targets using line scan cameras[END_REF][START_REF] Hansard | Cyclopean geometry of binocular vision[END_REF]. Furthermore, this online calibration update must be computationally efficient if it is to be used alongside other complex algorithms such as object tracking and 3D reconstruction.

Popeye is an active binocular robot (POP Consortium, 2008), which reproduces the sensory configuration of the human head (Fig. 1). The orientation of the robot head (pan/tilt) is controlled by motors, as is the direction of the eyes (version/vergence). The robot can therefore direct its attention towards a visual stimulus in its surroundings using combined stereo and tracking techniques [START_REF] Bellotto | Computationally efficient solutions for tracking people with a mobile robot: an experimental evaluation of Bayesian filters[END_REF]. The advantage of an active vision system is that information about the environment can be gathered more efficiently, by appropriately re-orienting the cameras. Hence visual information is combined with motor control, in a feedback loop, which enables the robot to react to a dynamic environment [START_REF] Bajcsy | Active perception[END_REF].

Previous Work

In order to preserve the calibration of the stereo system after camera vergence, one could continually recalibrate the system. This would, however, be very difficult to perform in real-time. A more efficient method would be to calibrate the system once, and then to update this calibration as the cameras move. This technique was used in [START_REF] Björkman | Real-time epipolar geometry estimation of binocular stereo heads[END_REF], where the epipolar geometry was updated by extracting information from feature correspondences in the current images. This required an iterative approach to estimate the stereo-head geometric parameters, which is computationally expensive. Furthermore, such an approach ignores kinematic information that can be used to improve the real-time performance of the system. The problem was also tackled in [START_REF] Hart | Epipolar geometry for humanoid robotic heads[END_REF] who developed an epipolar-kinematic model in which motor data is used to compute an updated representation of the epipolar geometry. The kinematics of the system are obtained by rotating the cameras and observing the relationship between the images of two views before and after the rotation. After this calibration procedure, the essential matrix [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]) can be updated, provided that the current settings of the camera motors are known. This approach allows computation of the essential matrix in real-time. Related methods can be used to coordinate an active camera with a static camera [START_REF] Horaud | Camera cooperation for achieving visual attention[END_REF].

The Popeye robot belongs to a class of active binocular robot heads with very accurate motors [START_REF] Shih | Calibration of an active binocular head[END_REF][START_REF] Aryananda | MERTZ: A quest for a robust and scalable active vision humanoid head robot[END_REF][START_REF] Beira | Design of the robot-cub (icub) head[END_REF][START_REF] Miwa | Development of a new human-like head robot we-4[END_REF]. Moreover, the design of the eye pan mechanical structures allows each camera position to be precisely adjusted to achieve close agreement between the optical and rotational centres of the cameras. This permits a direct approach to epipolar-update, as described in section 2. The accuracy and negligible backlash of the DC brushed motors ensures the repeatability of the new method, as demonstrated in section 3.

The present work is closely related to certain autocalibration procedures [START_REF] Hartley | Self-calibration of stationary cameras[END_REF][START_REF] Ruf | Visual servoing of robot manipulators Part I: Projective kinematics[END_REF][START_REF] Knight | Automated alignment of robotic pan-tilt camera units using vision[END_REF], in which constrained movements are used to estimate the camera parameters. The aim of the present work is somewhat different; we estimate a direct mapping from motor-settings to image-transformations, without explicitly estimating the calibration parameters.

The paper is organized as follows. The problem definition and main contributions of this work are stated in sections 1.2 and 1.3. Sections 2.1, 2.2 and 2.4 describe the estimation, analysis and synthesis of homographies, respectively. The new method is based on the statistical model (3), which allows for uncertainty in the homography estimates, as defined in section 2.3. Section 3 describes the results of experiments using Popeye. The conclusions of the work are stated in section 4.

Problem Definition

The problem to be solved is that of compensating for the effects of known camera-rotations on the images. This is important for both monocular and binocular tasks, as described below. The monocular geometry (which applies equally to the left and right cameras) will be described first.

The standard pinhole-camera model will be used to represent the imaging process. The scene and image points will be identified by homogeneous coordinates X (X, Y, Z, 1) and x (x, y, 1) respectively, where ' ' indicates equality up to a non-zero scale factor. If the pose of the camera, with respect to the scene coordinate-system, is represented by the 3 × 3 rotation matrix R and 3 × 1 translation-vector t, then x A(R |t)X , where A is the (invertible) uppertriangular matrix that contains the intrinsic parameters [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Now consider two views, V and V j , that differ by a rotation of the single camera described above. Specifically, suppose that the camera is aligned with the scene-coordinate system, and that image-points x ∈ V and x j ∈ V j are observed before and after a vergence rotation (monocular pan) of angle θ j . It follows that x A(I |0 )X and x j A(R j |0 )X and therefore

x j H j x where H j AR j A -1 .

(1)

The full-rank 3 × 3 matrix H j represents the homography [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]) that maps view V Fig. 2 A camera samples the bundle of rays that pass through the optical centre C . Two different views, separated by a rotation of θ j , are indicated by the planes V and V j .

Corresponding points x and x j , which are images of X , are related by a projective transformation H of the projective space P 2 .

to view V j , as shown in Fig. 2. The matrix H j will be analysed in section 2.2. Now, in the binocular case, let x and x r be corresponding points in the left and right images. The epipolar geometry of the binocular system is expressed by the constraint x r F x = 0 where F is the fundamental matrix [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. If the cameras are rotated, by angles θ j , θ rj , then the new coordinates x j , x rj will not be compatible with F . However, if homographies H j and H rj are available, then the points can be transformed back to the coordinate system of F , so that (H -1 rj x rj ) F (H -1 j x j ) = 0. Alternatively, F itself can be transformed, so that the general epipolar constraint is

x rj F j x j = 0 where F j = H - rj F H -1 j . (2)
The cases H 0 = I and H r0 = I are defined so that F 0 = F is the original fundamental matrix. The original problem, of compensating for binocular vergence, has now been reduced to that of estimating the homographies H j and H rj .

There are two possible ways to estimate these matrices. This is an essentially monocular task, and so the left/right subscripts will now be suppressed. In the image-based method, H j is estimated from 1 ≤ k ≤ N point-correspondences x k ↔ x jk . This method is accurate, but it is also slow, as an additional correspondence problem must be solved after every camera movement. Alternatively, in the motor-based method, the current motor-angle θ j is substituted into the underlying rotation matrix R j which appears in (1). This method is fast, but the results are inaccurate, owing to residual misalignment of the optical and rotational centres (see sec. 2.3).

Main Contributions

This paper shows that the image-based and motor-based approaches, as described above, can be combined. The resulting method has the accuracy of the former approach, as well as the speed of the latter. An outline of the new method is given below.

Firstly, in an offline calibration procedure, the imagebased method is used to estimate the homographies H j that are induced by a set of motor-angles θ j , where 1 ≤ j ≤ M , as described in section 2.1. Each estimated homography determines another angle φ j , which is extracted from the matrix H j , as described in section 2.2. The image-based and motor-based angles are theoretically equal, φ j = θ j . In practice, however, there is a systematic difference between these parameters, owing to the optical and mechanical effects described in section 2.3. These systematic effects can be represented by a statistical model f , and a vector of parameters η, such that φ ≈ f (θ, η).

(3)

If the parameters can be learned, then a homographyangle φ can be predicted from any current motor-angle θ, as described in section 2.3. The corresponding homography H θ can then be quickly constructed following a real-time online procedure described in section 2.4. This method is accurate, because the predicted homography H θ must be compatible with the original H j at all fitted values θ = θ j , where j indexes the M views that are used to estimate the model (3). This method is also fast, because no feature point extraction and matching is required after each camera rotation. This computational efficiency allows the system to fixate a moving target with calibrated cameras. Finally our proposed method is robust because at runtime, it does not depend on the scene texture, establishing feature point correspondences, or a generic scene assumption, unlike purely image-based methods (section 3.5). In summary, the main contribution of this paper is a novel method for real-time epipolar geometry update, which is applied to an active binocular robotic head. We demonstrate the validity of this homography-based method on the Popeye robot, which can simultaneously fixate an object, track it over time, and perform 3D reconstruction in real-time.

Methods

The update-procedure will now be presented, in order of execution. Most importantly, section 2.3 introduces a model and solution for the angle mapping (3).

Homography Estimation

The initial homographies H j are estimated by the standard direct linear transformation algorithm [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Each homography has eight degrees of freedom (3 × 3 minus scale). This means that k = 1, . . . , N correspondences x k ↔ x jk are required, where N ≥ 4. It can be shown [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]) that each correspondence defines a pair of homogeneous constraints

0 3 -w k x jk y k x jk w k x jk 0 3 -x k x jk h j = 0 2 (4)
where

x k = (x k , y k , w k ) , h j = H 11 j , H 12 j , . . . , H 33 j
, and 0 L is the column-vector of L zeros. The N instances of the 2 × 9 matrix on the left of ( 4) are stacked to form a single 2N × 9 matrix G j . A solution for H j can then be obtained from the right singular-vector of G j that corresponds to the smallest singular value.

Homography Analysis

Recall from (1) that H j AR j A -1 , where R j is a rotation matrix. This means that the homography H j is a conjugate rotation, having the same eigenvalues as R j [START_REF] Hartley | Self-calibration of stationary cameras[END_REF]. Furthermore, H j can be decomposed as

H j U D j U -1 (5)
where D j is a diagonal matrix of eigenvalues, and U j is a matrix of eigenvectors. The decomposition has the following structure [START_REF] Hartley | Self-calibration of stationary cameras[END_REF][START_REF] Ruf | Visual servoing of robot manipulators Part I: Projective kinematics[END_REF], which involves a complex-conjugate pair of eigenvalues:

D j diag λ j , λ * j , 1 where λ j = exp(iφ j ) (6) and U = (u, v , w ). ( 7
)
The angle φ j and axis w correspond to the rotation R j . The complex vectors u and v represent the circular points in the plane orthogonal to the rotation axis [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Under pure rotations, the axis of rotation is fixed, and therefore the matrix U remains unchanged. In practice, there are M eigenvector matrices U j , which are only approximately equal, owing to misalignment of the optical and rotational centres. Thus, a synthesis matrix Ū may be chosen as that associated with the largest motor-angle θ j , or estimated with a suitable averaging procedure in an offline learning stage. 1 Following this, the pre-computed eigenvector matrix Ū will be used in section 2.4 for homography synthesis in a real-time online procedure, since only λ j needs to be computed at runtime. The angle φ can be computed from the motor-rotation θ as described in the following section.

Motor-Based Parametrization

This section addresses the chief issue, which is the relationship (3) between the motor parameter θ, and the homography parameter φ. The model f in (3) must account for two systematic effects. Firstly, the existence of a 'pinhole' optical centre is only an approximation for real cameras. Secondly, the approximate optical centre may not coincide with the rotational centre of the system. This means that camera rotations will be accompanied by small translations [START_REF] Hayman | The effects of translational misalignment when self-calibrating rotating and zooming cameras[END_REF]. It is important to note that these effects need only be modelled over a limited range of angles, as determined by the mechanics of the system. Furthermore, this range has a natural origin, which corresponds to the straight-ahead camera position. All motor-counts θ are specified in relation to this origin. The preceding considerations suggest the zero-offset linear model φ ≈ η θ, which will be experimentally validated in section 3. This model states that the M available angle-pairs {θ j , φ j } are related by

φ j = η θ j + j where 1 ≤ j ≤ M. (8)
The random errors j are due to the optical and mechanical effects described above, as well as to lens distortion and feature mis-localisation. The unknown parameter η can be estimated, given one or more anglepairs. The estimation procedure requires the definition of a distance metric in which the observed discrepancies δ = φ -η θ can be minimized. Two possible metrics are |δ| R ∈ [0, ∞], the standard Euclidean metric on R 1 , and |δ| S ∈ [0, 1], an angular metric on the circle S 1 . These are respectively defined as

|δ| R = √ δ 2 and |δ| S = 1 2 (1 -cos δ). (9)
The Euclidean metric is not suitable for general angular problems because, for example, |δ| R = |δ + 2π| R .

1 The Ū matrix was obtained by computing columns ū and w as the leading eigenvectors of j u j u j and j w j w j respectively, with v = ū * .

However, using the Taylor approximation cos(δ) = 1 -

1 2 δ 2 + O(δ 4 ), it is clear that |δ| S = 1 4 |δ| 2 R + O(δ 4 ). (10)
Hence, for small values of δ, it is possible to use |δ| 2 R rather than |δ| S . This is important, because the Euclidean metric is much easier to work with. The approximation can be quantified by the relative error

E RS (δ) = abs |δ| S -1 4 |δ| 2 R |δ| S . (11)
The maximum relative error E RS (φ j -θ j ) is 0.12% for the data-set used here, and so the Euclidean metric will be used. It will be demonstrated in section 3.3 that the results of using |δ| S are, for practical purposes, identical.

The minimum of the Euclidean error can be found by setting the derivative d/dη of 1 4 N j |φ j -η θ j | 2 R to zero. This leads to the well-known estimate for the slope of a regression-line through the origin:

η = M j θ j φ j M j θ 2 j . ( 12
)
If the errors j in (8) are normally distributed, pr( j ) ∝ exp -1 2 2 j , then η in ( 12) is also the Maximum Likelihood estimate of the underlying parameter η.

It is now straightforward, given any motor angle θ, to compute a predicted homography angle

φ = η θ (13)
by analogy with the original model ( 8). This prediction, given the estimate η, does not involve any other computations.

Homography Synthesis

Suppose that a homography angle φ has been predicted from a motor-angle θ, using the fitted model (13). It is then possible, by inverting the procedure of section 2.2, to create a new homography H , as noted in [START_REF] Hartley | Multiple view geometry in computer vision[END_REF][START_REF] Knight | Automated alignment of robotic pan-tilt camera units using vision[END_REF]. First, by analogy with (6), the eigenvalue matrix is synthesized:

D = diag λ, λ * , 1 where λ = exp(iφ), (14)
and λ, λ * are complex conjugates. The coordinates u and v of the circular points (7) are independent of the rotation angle. It follows that, by analogy with (5), the synthesis matrix of eigenvectors Ū can be combined with D as follows,

H = Ū D Ū -1 (15)
Algorithm 1 Offline parameter estimation 1. Estimate homographies H j induced by a set of motorangles θ j (section 2.1).

2. Eigendecompose H j to obtain (i) eigenvalues λ j , from which the image-based rotation φ j is found, and (ii) a matrix of eigenvectors Ū (section 2.2).

3. Estimate parameters η relating image and motor-based rotation angles using (12).

Algorithm 2 Online epipolar geometry update 1. Calculate predicted homography angle φ using (13).

2. Create new homography H from synthesized eigenvector matrix Ū and eigenvalue matrix D as in ( 14).

3. Update fundamental matrix F j using the homographies H j and H rj representing the rotation of both cameras using (2).

which gives the estimated homography.

The above procedure is performed separately for the left and right cameras, so that H j and H rj are obtained from θ j and θ rj , respectively. It is emphasized that only (13,14 & 15) need to be computed at run-time, with Ū fixed in advance. The current fundamental matrix F j is then obtained from H - rj F H -1 j , as in (2). Both the offline and online parameter estimation procedures are detailed in algorithms 1 and 2 respectively.

Experiments and Results

In order to evaluate the new approach, a sequence of tests were carried out using the Popeye robot, as described in sections 3.1 and 3.2. The experiments were set in a laboratory environment in which the Popeye head was placed in a position to best view any surrounding activity. The resulting data is explored in section 3.3, and the assumptions of section 2 are validated. The performance of the new method is evaluated in section 3.4 using an image-based error measure, and compared to a purely image-based method (section 3.5).

Robot Hardware

The robot's vision system consists of two PointGrey colour cameras, which provide images of size 1024×768. The Popeye head has four rotational degrees of freedom, but only the left/right vergence motors are used for the experimental analysis. The rotations are performed by DC brushed motors, which are controlled by discrete (and repeatable) angle-steps. The mechanical system also provides a reference position (straight ahead).

Fig. 3 An illustration of the motor-camera calibration procedure. For each camera, a set of images are extracted which are representative of the viewable area. The relationship between the motor-angle θ and the homography-angle φ is modelled by the method described in section 2.3.

Calibration Procedure

A practical angular range of operation was chosen to be ±20 • around the reference position, as shown in Fig. 3. This range was split into a discrete set of nine views V j , each separated by ∆θ = 5 • . The resulting images are representative of the viewable area of the scene. Textured cards were fixed to the facing furniture, containing stable calibration features. Each camera was rotated in turn by ∆θ over the whole angular range of operation, at each step taking a snapshot of the viewable area. For each view, ten evenly distributed featurepoints were matched and verified by standard methods [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. Three calibration runs were performed for each of the two cameras, giving a total of six data-sets.

The homography matrices mapping each view V j to the fronto-parallel position V were then estimated, using the method of section 2.1. The eigen-decomposition of 2.2 was then used to extract the homography angle φ j corresponding to each motor angle θ j . The linear model f , that characterizes the relationship between the motor-angle θ and the homography-angle φ, was then estimated by the method of section 2.3.

Exploratory Analysis

The first task is to examine the relationship between the measured motor angles θ j and the estimated homography angles φ j . The mean absolute discrepancy |θ j -φ j | is 2.92 • , with a standard deviation of 1.85 • . The maximum absolute discrepancy is 6.87 • . The data from the first trial are plotted in Fig. 4. It can be seen q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Fig. 4 Visualization of the relationship between φ j (inner points) and θ j (outer points) for the left and right cameras.

The connecting lines would be purely radial if all {θ j , φ j } pairs were equal. The area of each translucent sector is proportional to the corresponding angular discrepancy |φ j -θ j |. that the discrepancy tends to increase for more extreme views, in agreement with the model φ ≈ ηθ in (8).

The data can now be used to validate the quadratic approximation |δ| R of the angular error |δ| S in (9). Note that the latter is a data-dependent sum of j = 1, . . . , M cosinusoids, 1 2 1-cos(φ j -ηθ j ) . The small range of the angular errors means that there is no significant difference between the minima of the two metrics, as shown in Fig. 5. The optimum η for the angular error (found by numerical minimization) differs from the regression estimate η in (12) by 1.9 × 10 -8 .

Having validated the error-metric, the adequacy of the simple model φ = ηθ + in (8) must also be confirmed. This can be done by plotting the predicted values φ θ = ηθ in the same form as the observed values in Fig. 4. It can be seen from Fig. 6 that, with respect to the predicted values, there is no systematic pattern in the residuals. This indicates that the linear model is adequate. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Fig. 6 Visualization of the relationship between φ j (inner points) and φ = ηθ j (outer points) for the left and right cameras. The angular errors (represented by the translucent sectors) are greatly reduced with respect to those in Fig. 4 (note that the inner points, representing φ j , are unchanged).

(a) A binocular image pair taken from the frontoparallel configuration (left camera -right camera).

(b) The left and right cameras are rotated by -10 and +12 degrees respectively, and the homographies representative of these rotations are calculated.

(c) The estimated homographies are then used to un-rotate the current images to the original frontoparallel views (note: the field of view is truncated). Lastly, the validity of the synthetic homographies can be checked visually, in the actual images. This is done by synthetically un-rotating the robot's current view, based on the known motor angles and fitted models, as illustrated in Fig. 7(a)-(c).

Quantitative Evaluation

Consider the views V and V j , separated by an angle θ j and related by homography H j . Recall that point x k in V is mapped to point x jk in V j as x jk H j x k (1). The error of each mapping H j will now be evaluated over the N point-pairs in each data-set. The points are subject to errors in both images, and so the symmetric transfer error [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]) is used:

E 2 H j = 1 2N N k=1 d 2 x k , H -1 j x jk + d 2 x jk , H j x k . (16)
Here, for example, d 2 x jk , H j x k is the squared Euclidean image-distance between the measured point x jk in V j and the point H j x k that is mapped from the reference view V. This criterion will be used to quantitatively evaluate the accuracy of the synthesized homographies that were obtained from the motor-camera model. Matched points were extracted from the three data-sets, with left and right views spanning ±20 • of pan. For each set of images, the motor-camera calibration procedure was carried out as described in section 3.2. The estimated linear relationship between the homography-angles φ j and the motor angles θ j can be seen in Fig. 8.

The model parameters were then used to synthesize the homographies needed to map the scene points from the fronto-parallel images V to the rotated views V j at angles θ j . For a comparison between our motor-image based and a purely image-based method, we additionally estimated the homographies with the standard direct linear transformation (DLT) algorithm [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. The actual and predicted scene points obtained with both methods were then compared in terms of symmetric transfer error, as defined in (16). A discussion of the results obtained with both methods is presented in the following section.

Discussion

The results obtained after computing the RMS transfer error for our motor-image method, over all six data sets, 2 was 2.09 pixels (in the 1024 × 768 images). The maximum and standard deviation of the transfer errors were 6.68 and 1.16 pixels, respectively. The purely image-based method achieved an RMS transfer error of 1.03, and a maximum and standard deviation of 3.28 and 0.58 pixels respectively. The plots in Fig. 9 show a breakdown of the errors across the two cameras, nine angles and three experimental trials.

Note that image-based results are representative of an ideal case, in which sufficient stable points and perfect correspondences are available to the algorithm. In practice, an automatic interest point detector and correspondence algorithm will have to cope with low textured scenes, and imperfect matches. Fig. 8 Homography-angles φ j plotted as a function of the motor-angles θ j . Three trials are shown (rows) for each camera (columns). The least-squares regression model is drawn through the data points. It can be seen that the homography/motor relationship is approximately linear, with a slope η < 1. The dashed line, for comparison, has unit slope.

Left Camera

The results from Fig. 9 show that in an ideal case, the purely image-based method will produce lower pixel errors. However, a collection of independently fitted homographies has many more parameters than our synthesis model, so indeed they should produce a better fit. Let N c and N r denote the number of monocular cameras and monocular camera rotations respectively. Then, the number of estimated parameters for the image-based method is 8 × N r × N c , estimated online (H has 8 degrees of freedom). For the motor-image method, the number of parameters is (1 + 8) × N c , estimated offline (1 dof from φ, and 8 from the U matrix). U has at most eight degrees of freedom, as it is obtained from the 3 × 3 matrix H , which is defined up to an overall scale. More importantly, it is clearly seen that the number of parameters estimated by our method does not depend on the number of camera rotations, N r . For example in our experiments, the number Note that H 0 = I , by definition, and so there is no error at θ = 0. The image-based method (green) represents the best case scenario where there are sufficient stable points, prefect correspondences, and necessary scene structure to avoid degenerate cases. Our motor-image approach (black) shows performance in practice and does not require point correspondences, scene texture or a particular structure.

Left Camera

of estimated parameters totalled 8 × 8 × 2 = 128 for the image based method, and (1 + 8) × 2 = 18 for the motor-camera method.

It is additionally important to point out that for our motor-image based technique, the scene-structure or texture is of no particular importance; each image-point could represent any scene-point on the corresponding ray (as indicated in Fig. 2). Hence the motor-image based results generalize immediately to any scene.

In contrast, purely image-based methods do depend on the scene texture and structure to compute feature correspondences and to avoid degenerate cases. In the case where F is re-estimated online, at least n ≥ 7 point correspondences are needed after each camera ro-tation, subject to the points being in general 3D position. This means that the points must not lie on (or near) a plane or other 'degenerate' surface, [START_REF] Hartley | Multiple view geometry in computer vision[END_REF]. The image based method to which we compare involves the re-estimation of homographies to update F . In order to compute each H , at least four (no three collinear) point correspondences per camera rotation are required. Since the motor-image method does not depend on computing H from feature point correspondences, it is not susceptible to the degeneracies described above, making it invariant to scene texture and structure. A complete theoretical error analysis of the epipolar geometry [START_REF] Csurka | Characterizing the uncertainty of the fundamental matrix[END_REF][START_REF] Brandt | On the probabilistic epipolar geometry[END_REF] is beyond the scope of this work.

In these experiments we have shown that our method accounts for small misalignments in the optical and rotational camera centres. The development of a similar approach which accounts for the purposeful misalignment of camera rotation and optical centres [START_REF] Hayman | The effects of translational misalignment when self-calibrating rotating and zooming cameras[END_REF] is left for future work.

Target tracking and 3D reconstruction

In order to validate the homography-based method in a real-world setting, the algorithm was used to extend the capabilities of the Popeye robot by allowing sparse 3D reconstruction during camera vergence movements. Whereas previously the cameras of the robot were static and the viewable area constrained to the initial calibration position of the cameras, it is now possible to allow the cameras to move and verge on a target of interest whilst maintaining its calibrated state. A simple gaze control feedback algorithm was implemented to keep each camera centred on a target of interest as it moved. In this scenario, the robot was placed in a room and its task was to keep track of a human face whilst performing sparse 3D reconstruction of this area of interest.

The possible movements allowed by the binocular robot are horizontal head pan rotation, a vertical tilt movement, and independent pan movements for the stereo cameras. The camera pan movements are considered to be independent and we consider the task of bringing the target position to the centre of the cameras as a monocular tracking task, unlike the tracking in [START_REF] Pagel | Self calibration of the fixation movement of a stereo camera head[END_REF], in which the camera movements were coupled. This leads to a simple update procedure to bring the object of interest to the centre of each camera:

δθ ∝ arctan δ x f x δψ ∝ arctan δ y f y (17)
where δθ and δψ are the damped pan and tilt eye motor rotations necessary to keep track of a moving target Fig. 10 Response of the active binocular head motors to a step input. A static face was detected and the cameras rotated to fixate on the face. The damping factor was tuned so as to give a fast rise-time with minimal overshoot.

under stability. The damping factor for each motor ξ was found by giving a step-input in the form of a fixed face target to the active head. The motor response was plotted for varying ξ in Fig. 10, and the parameters close to critical damping were chosen.

The commonly used Viola-Jones face detection algorithm [START_REF] Viola | Robust real-time face detection[END_REF] was used to return the position of a single face in each of the camera images. Feature points within this face area were found using an interest point detector [START_REF] Harris | A Combined Corner and Edge Detector[END_REF]. A more recent method inspired by human visual attention [START_REF] Frintrop | Active gaze control for attentional visual SLAM[END_REF] can easily be incorporated to determine a region of interest in the image. The detected feature points were matched after limiting the search-space, for each point, to a region around the conjugate epipolar line. The matched points were then triangulated and plotted in euclidean space as shown in Fig. 11(a)-(c). The complete system runs in real-time (approx. 6fps), with all image-processing implemented on standard PC hardware. A demonstration of the system in operation can be seen in the attached multimedia file3 .

Conclusion

The results described in section 3.4 show that the new method of homography-generation is sufficiently accurate for practical use. For example, the method can be used to predict the coordinates of image-features, after controlled rotations of the cameras. The method is also useful for binocular vision, in which the fundamental matrix is updated as described in section 1.2. Furthermore, the new method is also fast, because the mapping (c) As the target moves closer to the camera, the vergence angle between the cameras is updated to continually maintain fixation. from motor settings to image transformations is only estimated once. No additional point-correspondences are needed at run-time, and so computationally expensive feature-extraction and matching is avoided. The methods described here have enabled active tracking and image-stabilization to be performed in real-time by the Popeye robot. Future work will consider more complicated models for the mapping φ ≈ f (θ, η) which relates the motor and homography angles (3). It is expected that the method can be extended, in this way, to systems that have poor alignment between the optical and rotational centres. Finally, the idea of estimating the relationship between control-based and image-based parameters could be applied to other visual processes, such as active zoom.

Fig. 1

 1 Fig. 1 The Popeye active binocular robotic head (designed and built in the Institute for Systems & Robotics, University of Coimbra, Portugal). The binocular vergence is controlled by a pair of motors, one below each camera. The pan and tilt angles of the head are controlled by another pair of motors.

Fig. 5

 5 Fig.5Comparison of the angular (dashed) and Euclidean (solid) metrics (9), on the actual data-set. The upper plot shows periodic form of the angular error. The global minimum is well-approximated by a parabola, as in (10). The least-squares estimate (12) is the minimum of the parabola, indicated by the vertical line. For the present data, the difference between the two metrics is insignificant. The lower plot shows a magnified view.

Fig. 7

 7 Fig. 7 Visual validation of the synthesized homographies.

Fig. 9

 9 Fig.9Symmetric transfer error (16) as a function of the homography angle θ j for the motor-image method (black), and the purely image-based method (green) (images of size 1024 × 768). Three trials are shown (rows) for each camera (columns), as in Fig.8. The vertical bars indicate the spread ±σ of the errors around the corresponding mean. Note that H 0 = I , by definition, and so there is no error at θ = 0. The image-based method (green) represents the best case scenario where there are sufficient stable points, prefect correspondences, and necessary scene structure to avoid degenerate cases. Our motor-image approach (black) shows performance in practice and does not require point correspondences, scene texture or a particular structure.

  (a) A snapshot of the sparse reconstruction (left) of the plane formed by a textured poster (right) fixed to the wall of the lab. (b) The active head has detected a face and begins to track its movement whilst reconstructing matched points in this region of interest.

Fig. 11 A

 11 Fig. 11 A 3D viewer (left) and the image captured from the left camera of the binocular stereo pair (right). In all three camera images, the matched points (red dots) are plotted in euclidean space and seen with the 3D viewer, which is updated in real-time. The motion of the head, for example between 11(b) and 11(c) is well reconstructed despite the camera movement and varying vergence angle.

The statistics were computed after excluding the 'perfect' values at θ = 0.

Visit: http://www.youtube.com/watch?v=jcIK8AMoejo