Xavier Alameda-Pineda

Radu Horaud

GEOMETRICALLY-CONSTRAINED ROBUST TIME DELAY ESTIMATION USING NON-COPLANAR MICROPHONE ARRAYS

In this paper we present a geometrically-constrained time delay estimation method for sound source localization (gTDE). An algebraic analysis reveals that the method can deal with an arbitrary number of non-coplanar microphones. We derive a constrained non-linear optimization problem that can be solved using local convex programming. Unlike existing techniques, which consider pairwise TDE's, the proposed method optimally estimates a set of time delays that are consistent with the source's location. Extensive simulated experiments validate the method in the presence of noise and of reverberations.

INTRODUCTION

For the last decades, source localization from time delay estimates (TDEs) has proven to be an extremely useful methodology with a variety of applications in such diverse fields as aeronautics, telecommunications and robotics. Also referred to as multilateration, this problem is highly related to the one of estimating time delays. We are particularly interested in the development of a general-purpose TDE-based method for sound-source localization in indoor environments, e.g, human-robot interaction, ad-hoc teleconferencing using microphone arrays, etc. This type of consumer-oriented applications are extremely challenging for several reasons: (i) there may be several sound sources and their number varies over time, (ii) regular rooms are echoic, thus leading to reverberations, and (iii) the microphones are often embedded in devices (robot heads, smart phones, etc.) generating highlevel noise.

The TDE problem has been very well investigated and a recent review can be found in [START_REF] Chen | Time Delay Estimation in Room Acoustic Environments: An Overview[END_REF]. The vast majority of existing approaches deals with a microphone pair but it is not straightforward to extend most of these methods to more than two microphones. Methods addressing multichannel TDE can be roughly divided into two categories: methods estimating the acoustic impulse responses and methods exploiting the redundancy among several microphones. [START_REF] Doclo | Robust Adaptive Time De-lay Estimation for Speaker Localization in Noisy and Reverberant Acoustic Environments[END_REF] is illustrative of This work was supported by the EU project HUMAVIPS FP7- ICT-2009-247525. the first category where a method based on generalized eigenvalue decomposition is proposed. The second category is represented by [START_REF] Chen | Robust time delay estimation exploiting redundancy among multiple microphones[END_REF] where a multichannel criterion based on crosscorrelation is proposed to estimate time delays using a linear microphone array. In both cases, experiments are performed on speech data in a simulated indoor environment.

As already mentioned, an alternative to TDE is multilateration, which makes assumptions about the time delay estimates. This provides a framework for casting the problem into maximum-likelihood estimation or into mean-squared error minimization (see [START_REF] Pertilä | Acoustic Source Localization in a Room Environment and at Moderate Distances[END_REF] for a review). Two recent methods deserve to be mentioned. In [START_REF] Brutti | Comparison between different sound source localization techniques[END_REF] the authors use the acoustic maps together with the GCC-PHAT technique to localize sound sources from TDE's. The model in [START_REF] Ribeiro | Using reverberation to improve range and elevation discrimination for small array sound source localization[END_REF] includes the reverberations in order to enhance the localization performance while using a uniform circular array of microphones.

In this paper we propose a method that combines multichannel time delay estimation with source localization. More precisely, we cast the simultaneous estimation of time-delays and of source localization into a constrained optimization problem. We provide a detailed algebraic analysis of the proposed formulation, thus allowing us to estimate TDE values consistent with a source location. We show that the method can be used in conjunction with an arbitrary number of noncoplanar microphones. We describe a practical algorithm for solving the optimization problem at hand and we provide extensive experimental results.

The remainder of the paper is organized as follows. The signal and geometric generative models are described in Section 2. The proposed method is derived in Section 3 while implementation details, extensive experiments, and results are provided in Section 4.

SIGNAL AND GEOMETRIC MODELS

In this section we describe the signal acquisition model and the geometric model allowing to relate time delays with the relative position between source and microphones. We introduce the following notations: the position of the sound source S ∈ R N , the number of microphones M , as well as their positions, {M m } m=M m=1 ∈ R N . Let x(t) be the signal emitted by the source. The signal received at the m-th microphone writes:

x m (t) = x(t -t m ) + n m (t), (1)
where n m is the noise associated with the m-th microphone and t m is the time-of-arrival from the source to that microphone. The microphones' noise signals are assumed to be zero-mean independent Gaussian random processes. Throughout this paper, constant sound propagation speed is assumed, denoted by ν. Hence we write t m = S -M m /ν. Using this model, the expression for the time delay between the m-th and the n-th microphones, denoted by t m,n , writes:

t m,n = t n -t m = S -M n -S -M m ν . (2)
Notice that, for a fixed value of t m,n , the sound source generating t m,n lies in one sheet of a two-sheet hyperboloid with foci M n and M m . The sign of t m,n determines in which of the two sheets lies the sound source. Another remarkable property of the set of time-delays is that they are not independent, i.e., the relation t m,n = t m,k + t k,n holds for all k, m, n.

The signal model ( 1) together with the geometric generative model ( 2) allow us to cast the TDE problem into a constrained non-linear optimization problem, as explained in the next section.

PROPOSED METHOD

The proposed solution for multichannel TDE is described in the detail below. First, we show how the criterion used in [START_REF] Chen | Robust time delay estimation exploiting redundancy among multiple microphones[END_REF] for multichannel TDE in the case of linear microphone arrays can be used when the geometry of the array is not known beforehand. Second, we present how the knowledge of the microphones' positions can be used to constrain the algorithm. Finally, we summarize the proposed solution by outlining the optimization problem to solve.

Unconstrained TDE

The criterion used in [START_REF] Chen | Robust time delay estimation exploiting redundancy among multiple microphones[END_REF] was built from the theory of linear predictors. We outline, in the following, the basic steps to obtain the criterion to optimize for unconstrained multichannel TDE.

Given the M received signals {x m (t)} m=M m=1 , we would like to estimate the time delays between them. As explained before, only M -1 of the delays are independent. Without loss of generality we choose the delays t 1,2 , . . . , t 1,m , . . . , t 1,M . We select x 1 (t) as the reference signal and set the following prediction error:

e c,t (t) = x 1 (t) - M m=2 c 1,m x m (t + t 1,m ), (3)
where c = (c 1,2 , . . . , c 1,m , . . . , c 1,M ) T is the vector of the prediction coefficients and t = (t 1,2 , . . . , t 1,m , . . . , t 1,M ) T is the vector of the prediction time delays. Notice also that the signals x m (t + t 1,m ) and x n (t + t 1,n ) are on phase. The criterion to minimize is the expected energy of the prediction error in [START_REF] Chen | Robust time delay estimation exploiting redundancy among multiple microphones[END_REF], which is equivalent to (see [START_REF] Chen | Robust time delay estimation exploiting redundancy among multiple microphones[END_REF]):

t * = arg min t J(t), (4)
where J(t) = det (R(t)) with R(t) ∈ R M ×M being the real matrix of normalized cross-correlation functions evaluated at t. That is R(t) = [ρ i,j (t)] ij with:

ρ i,j (t) = E {x i (t + t 1,i )x j (t + t 1,j )} E {x 2 i (t)} E x 2 j (t)
,

where E{x} denotes the expectation of x.

This is how the time delay estimation problem is cast into a non-linear optimization problem. The problem is multivariate due to the fact that the signal model does not encode the geometry of the array. Solving the optimization provides for a set of consistent TDEs. Moreover, in the next section we show how to constrain the multichannel TDE criterion presented above, by means of the microphones' positions.

Geometrically-Constrained TDE

These values of t that correspond to a sound source position can be fully characterized by means of the geometric generative model presented in [START_REF] Doclo | Robust Adaptive Time De-lay Estimation for Speaker Localization in Noisy and Reverberant Acoustic Environments[END_REF]. In [START_REF] Chan | A simple and efficient estimator for hyperbolic location[END_REF] a solution is proposed in the 2D case (N = 2), namely the source and the microphones are coplanar. In this paper we propose a generalization to N = 3.

The M -1 constraints defined by (2) can be rewritten in the form a linear system for the sound source position S:

MS = K -ν 2 t 2 -2νtd 1
where M ∈ R (M -1)×N is a real matrix with its m-th row, 1 ≤ m ≤ M -1, given by (M m+1 -M 1 ) T , K ∈ R M -1 is a vector with its m-th element given by M m+1 2 -M 1 2 , t 2 denotes the component-wise square power of t and d 1 is the distance from the sound source to the first microphone. If M ≥ N and if M is a full rank matrix we can write:

S = Ad 1 + B, (5)
where

A = -2νM † t and B = -M † (K -ν 2 t 2 ), with M † = M T M -1 M T being the left inverse of M. Notice
that M is full rank if and only if the M sensors do not lie in the same hyperplane of R N . By replacing S from (5) in d 2 1 = S -M 1 , S -M 1 we obtain a second-degree equation in the unknown d 1 :

d 2 1 ( A 2 -1)+2 A, B -M 1 d 1 + B -M 1 2 = 0. ( 6
)
If this equation has two (or one) real solutions, the consistent set of TDEs, i.e., vector t, corresponds to a sound source position. It can be shown that only one of the two real solutions of (6) corresponds to t (the proof of this statement is left out due to space limitations). Therefore, the geometric constraint determining the feasibility of t is the non-negativeness of the discriminant of (6):

0 ≤ ∆ = A, B -M 1 2 -B -M 1 2 ( A 2 -1). ( 7
)
Notice that ∆ depends on the known microphone locations {M m } m=M m=1 , on the constant speed ν and on the variables t. Hence, given these parameters, we write ∆ = ∆(t).

To summarize, the minimization of ( 4) becomes the following non-linear constrained optimization problem:

t * = arg min t J(t) ∆ (t) ≥ 0 (8)
Notice that the optimal solutions of the constrained and unconstrained problems are exactly the same.

IMPLEMENTATION, EXPERIMENTS & RESULTS

The minimization of ( 8) is carried out using a publicly available MATLAB implementation [START_REF] Carbonetto | MATLAB primal-dual interior-point solver for convex programs with constraints[END_REF] of the log-barrier interior point method [START_REF] Boyd | Convex Optimization[END_REF]. This method is designed for continuous convex optimization problems. On one hand, it is likely to fail in finding the global optimum of non-convex problems such as [START_REF] Carbonetto | MATLAB primal-dual interior-point solver for convex programs with constraints[END_REF]. To overcome this issue, our algorithm starts from several initial points, i.e., the set S I = {t I i } P i=1 . For each one of these initializations a local minimum is found, then the minimum over these local minima is selected. In our simulations, P = 4096 and S I consists of points placed in a regular rectangular grid. Since the t m,n delays have upper and lower bounds: M m -M n , -M m -M n respectively, the grid limits are defined by the geometry of the problem. On the other hand, the function to optimize must be continuous and the signals are discrete. We chose to compute the normalized cross-correlation function of the linear interpolation of the discrete signals.

In order to accurately evaluate and validate the proposed gTDE method 1 , we developed a formal evaluation protocol using simulated data. A 3 × 4 × 2.5 meter room, with uniform absorption coefficients, was simulated using the state-of-the-art Image-Source Model (ISM) [START_REF] Lehmann | Prediction of energy decay in room impulse responses simulated with an image-source model[END_REF] available from [START_REF] Lehmann | Matlab code for imagesource model in room acoustics[END_REF]. The main parameter of this model is T 60 , which corresponds to the time needed for an energy decay of 60 dB. We simulated four microphones placed at (in meters): M 1 = (2.35, 1.25, 1.179) T , M 2 = (2.15, 1.25, 1.179) T , M 3 = (2.25, 1.35, 1.32) T and M 4 = (2.25, 1.15, 1.32) T , i.e., forming a regular tetrahedron. The sound source was 1 The software is available at http://gtde.berlios.de.

placed at 27 different positions, namely all the possible 3-tuples S = (x, y, z) T with x ∈ {0.825, 1.5, 2.175}, y ∈ {1.1, 2, 2.9} and z ∈ {0.6875, 1.25, 1.8125} (in meters). The source emitted speech fragments randomly chosen from [START_REF] Garofolo | Timit acoustic-phonetic continuous speech corpus[END_REF]. One hundred millisecond cuts of these sounds are the input of the evaluated method. We assume that only one source emitting within each cut. Finally the sensor noise, whose power depends on the chosen SNR, is added to these cuts.

Fig. 1 shows the results obtained with four different methods: pair-wise independent estimation of the time delays (bypairs), estimation using multichannel information without minimization (init), estimation based on unconstrained time delay minimization of (4) (tde), and the proposed geometrically-constrained minimization (gtde). Fig. 1(a) plots the percentage of good (non-anomalous) estimates, i.e., with an absolute error smaller than 100 µs and Fig. 1(b) plots the standard deviation of the good estimates as a function of the signal-to-noise ratio. These results correspond to an anechoic setup, T 60 = 0. As expected, the proposed method significantly improves the percentage of good TDEs while it Additional simulations were carried out in order to precisely evaluate the performance of the gTDE method in the presence of noise and of reverberations. Figures 2 and3 show the results on time delay estimation and sound source localization for different levels of noise and reverberations. In all the plots, the x-axis represents the SNR value (dB). The color corresponds to the method used: green for bypairs and blue for gtde; and the line style corresponds to the level of reverberation: solid-circle for T 60 = 0 s and dashed-square for T 60 = 0.1 s.

Regarding the TDE, in Fig. 2(a) the percentage of nonanomalous TDE is plot and in Fig. 2(b) the standard deviation of such time delays estimation is shown. Notice how the performance of both method improves with the SNR. Also, the proposed method clearly outperforms the baseline. This is not a surprise since it is using all the available information to consistently estimate the TDEs at a time. A remarkable fact is that the proposed method under reverberant conditions performs similarly than the baseline method in the anechoic case. Also, with higher percentage of non-anomalous estimates, the gtde method has lower error standard deviation.

Concerning the localization, Fig. 3(a) plots the percentage of localization inliers and Fig. 3(b) the standard deviation of the angular error. A sound source is considered to be an inlier if the angular absolute error is less than 30 • . As in the case of the TDE, the methods' performance improve with the SNR. It is worth noticing also that the proposed method outperforms the baseline method, even when that one is under echoic conditions and this one is under anechoic conditions.

Generally speaking, the methods perform as expected. The higher the SNR value the better the methods estimate the time delays, the higher the percentage of inliers and the lower the localization error. We can also observe a clear trend with respect to the reverberation level: the methods' performance decreases with T 60 . However the SNR and the T 60 havev different effects on the function to minimize. On one side, the sensor noise decorrelates the microphones' signals leading to much more (and randomly spread) local minima and increasing the value of the true minimum. If this effect is extreme, the hope for a good estimate decreases fast. On the other side, the reverberations produce only a few strong local minima. This perturbation is systematic given the source position in the room. Hence, there is hope to learn the effect of such reverberations in order to improve the quality of the estimates. These types of perturbations (noise and reverberations) of the function to minimize have clearly different effects on the results. Notice that the reverberation level has almost no effect on the quality of the estimates when the SNR is low. However, when this random effect disappears, i.e. higher values of SNR, a systematic and significative difference appears both in the time delay and in the localization estimates. Finally, the authors would like to remark on the method's robustness. For moderate levels of noise and reverberations (T 60 ≤ 100 ms, SNR ≥ 0 dB) the method is able to localize the sound source with mean angular squared error of 6 • in more than 60 % of the cases.

CONCLUSIONS & FUTURE WORK

In this paper, a new method on time delay estimation for sound source localization working on non-coplanar microphone arrays is presented. The estimation is cast into a multivariate optimization problem. In addition, the geometric model for the time delays add a non-linear constraint. The optimal value of the problem is a consistent set of time delays, useful to localize the sound source. Experiments on simulated data show the quality of the method and validate the approach.

The are several ways to extend this work. As outlined before, it would be very useful to learn the effect the reverbera- tions have on the objective function as in [START_REF] Ribeiro | Using reverberation to improve range and elevation discrimination for small array sound source localization[END_REF]. Also, it is worth to consider the multiple source case, following approaches like [START_REF] Chen | Acoustic source localization and beamforming: theory and practice[END_REF]. Besides that, a frequency decomposition stage may be useful to avoid the analysis in non-informative frequency bands ( [START_REF] Valin | Robust 3D localization and tracking of sound sources using beamforming and particle filtering[END_REF]). Also, experiments on more reverberant data, and on real data have to be done in order to explore the real extent of these initial results. Last but not least, it would be interesting to explore cases in which the microphones' positions have some error, and see how to adapt the method to estimate the time delays and correct the microphones' positions.

  (a) Good (non-anomalous) TDEs (b) TDE standard deviation

Fig. 1 .

 1 Fig. 1. (Best seen in color) TDE method comparison based on random speech fragments. The plots show (a) the percentage of good TDEs and (c) the standard deviation as a function of SNR in an anechoic setup (T 60 = 0).

  (a) Non-anomalous TDE (b) TDE Standard Deviation

Fig. 2 .

 2 Fig. 2. Evaluation of the TDE performance of gTDE method. The x-axis corresponds to the SNR value (dB), the color to the methods (blue for gtde, green for bypairs), and the line style to the reverberation level (solid-circle for T 60 = 0 and dashed-square for T 60 = 0.1 s). (a) shows the percentage of non-anomalous TDE and (b) the standard deviation of this estimates.

Fig. 3 .

 3 Fig. 3. Evaluation of the localization performance of gTDE method. The x-axis corresponds to the SNR value (dB), the color to the methods (blue for gtde, green for bypairs), and the line style to the reverberation level (solid-circle for T 60 = 0 and dashed-square for T 60 = 0.1 s). (a) shows the percentage of localization inliers and (b) the mean squared error of localization error.

20th European Signal Processing Conference (EUSIPCO 2012) Bucharest, Romania, August 27 -31, 2012 © EURASIP, 2012 -ISSN 2076-1465


Jordi Sanchez-Riera

Xavier Alameda-Pineda

Johannes Wienke

Antoine Deleforge

Soraya Arias

Jan Čech

Sebastian Wrede
email: swrede@cor-lab.uni-bielefeld.de

Radu Horaud

Online Multimodal Speaker Detection for Humanoid Robots

In this paper we address the problem of audiovisual speaker detection. We introduce an online system working on the humanoid robot NAO. The scene is perceived with two cameras and two microphones. A multimodal Gaussian mixture model (mGMM) fuses the information extracted from the auditory and visual sensors and detects the most probable audio-visual object, e.g., a person emitting a sound, in the 3D space. The system is implemented on top of a platformindependent middleware and it is able to process the information online (17Hz). A detailed description of the system and its implementation are provided, with special emphasis on the online processing issues and the proposed solutions. Experimental validation, performed with five different scenarios, show that that the proposed method opens the door to robust humanrobot interaction scenarios.

I. INTRODUCTION

Humanoid robots acting in populated spaces require a large variety of communication skills. Perceptive, proprioceptive as well as motor abilities are mandatory to make the information flow natural between people and robots participating in interaction tasks. On the perceptive side of the communication process, the tasks are mainly detection, localization and recognition. Depending on the available sensory modalities, the robot should be able to perform tasks such as: sound/speech detection/recognition, action/gesture recognition, identity/voice recognition, face detection/recognition, etc. Moreover, if several modalities are combined, multimodal tasks such as audio-visual event recognition or audio-visual speech processing are known to be more robust than unimodal processing and hence multisensory perception can drastically improve the performance of a large variety of human-robot interaction activities.

The problem of data fusion and multi sensory integration has been recognized for a long time as being a key ingredient of an intelligent system, e.g., [START_REF] Luo | Multisensor integration and fusion in intelligent systems[END_REF]. More recently, multimodal integration has been used in action recognition applications [START_REF] Lili | A framework for human action detection via extraction of multimodal features[END_REF], [START_REF] Wu | Realistic human action recognition with audio context[END_REF]. In these papers, the authors exploit the fact that for some actions such as "talk phone" the auditory information is relevant for describing the action. Another multimodal approach was followed in [START_REF] Lacheze | Audio/video fusion for objects recognition[END_REF], in which the auditory information was used to recognize objects that can Work supported by the EC project HUMAVIPS FP7-ICT-2009-247525 Fig. 1: A typical scenario in which a companion humanoid robot (NAO) performs audio-visual fusion in an attempt to detect the auditory status of each one of the speakers in the room. The system described in this paper processes the raw data gathered with the robot's camera and microphone pairs. The system output is a speaking probability of each one of the actors together with the 3D location of the actors' faces.

be partially occluded or difficult to detect. Notice that, the visual information is also very helpful when the auditory data is strongly corrupted by noise or by multiple sound sources. Another example can be found in [START_REF] Itohara | Particle-filter based audio-visual beat-tracking for music robot ensemble with human guitarist[END_REF] where the authors combine information coming from the two modalities to perform beat tracking of a person playing the guitar. Auditory and visual information is combined together to better address the problems of beat-tracking, tempo changes, and varying note lengths. Also this different-modality combination is used for improvement of simultaneous speach signals in [START_REF] Nakadai | Improvement of recognition of simultaneous speech signals using av integration and scattering theory for humanoid robots[END_REF]. Using a pair of cameras faces are first detected and then in 3D. Using two microphones sound-source separation by ADPF (Active Direction Pass Filter) is applied. Finally these data is integrated at two levels, at the signal level and at the word level. Among all possible applications using audio-visual data, we are interested in detecting multiple speakers in informal scenarios. A typical example of such a scenario is shown in Fig. 1, in which two people are sitting and chatting in front of the robot. The robot's primary task (prior to speech recognition, language understanding, and dialog) consists in retrieving the auditory status of several speakers along time. This allows the robot to concentrate its attention onto one of the speakers, i.e., turn its head towards in the speaker's direction to optimize the emitter-to-receiver pathway, and attempt to extract the relevant auditory and visual data from the incoming signals. We note that this problem cannot be solved within the traditional human-computer interaction paradigm which is based on tethered interaction (the user must wear a close-range microphone) and which primarily works in the single-person-to-robot communication case. This considerably limits the range of potential interactions between robots and people engaged in a cooperative task or simply in a multi-party dialog. In this paper we investigate untethered interaction thus allowing a robot with its onboard sensors to perceive the status of several people at once and to communicate with them in the most natural way.

The original contribution of this paper is a complete realtime audio-visual speaker detection and localization system that is based binocular and binaural robot perception as well as on a generative probabilistic model able to fuse data gathered with camera and microphone pairs. The remainder of the paper is organized as follows. Section II describes previous work in human-robot interaction and audio-visual processing that is closely related to our approach. Next, the audio-visual fusion model is depicted in section III; this model is specifically designed to handle visual and auditory observations gathered at different sampling rates and which live in different observation spaces. Being based on a mixture model, the proposed fusion framework is able to estimate the number of audio-visual objects (speakers) in the scene and to deal with intermittent auditory data. Sections III-A and III-B briefly describe the visual processing and auditory processing modules while the audio-visual data alignment method, i.e., calibration is described in section IV. The overall system architecture is described in section V together with implementation details. Results and experiments are described in section VI. Finally, section VII draws some conclusions and discusses topics for future research.

II. RELATED WORK AND CONTRIBUTIONS

Audio-visual processing has been studied by many researchers. In [START_REF] Beal | Audio-visual sensor fusion with probabilistic graphical models[END_REF] the authors describe a speaker detection probabilistic graphical model fusing the information coming from one camera and two microphones. An EM algorithm estimates the model's parameters, i.e., the audio-visual appearance and the position of the speaker. In [START_REF] Fisher | Speaker association with signal-level audiovisual fusion[END_REF] the author proposes to use maximally informative projections to retrieve the main speaker. One camera per potential speaker and one microphone are used to gather the raw data, which is projected in order to subsequently select the speaker, based on information-theoretic criteria. This approach is well suited for applications such as video-conferencing.

A second group of methods deal with interaction in smartroom environments. These methods assume the existence of several sensors distributed in the scene. For instance, [START_REF] Shivappa | Auvio-visual fusion and tracking with multilevel iterative decoding: Framework and experimental evaluation[END_REF] uses data acquired in a room equipped with a multi-camera system and an array of microphones. A tracking system is developed to complement information for a room with a smart interaction environment. The authors of [START_REF] Zhang | Boosting-based multimodal speaker detection for distributed meeting videos[END_REF] present an application aiming at making meetings more dynamic for people who are remotely connected. Based on one camera and one microphone array, this methodology is able to detect the speaking persons in real-time.

Because of on-line and on-board processing constraints associated with humanoid platforms, the computational load and complexity are constraints that need to be taken into account. Furthermore, the robot does not have a distributed sensor network, but merely a few sensors, which are all located in its headan agent-centered sensor architecture. Hence, one should achieve a trade-off between performance and complexity.

In this paper we present both a novel method and an original system approach to tackle the problem of on-line audiovisual detection of multiple speakers using the companion humanoid robot NAO 1 . The proposed method uses data coming from a stereo pair of cameras and two microphones. Implemented on a hardware-and sensor-independent middleware, the software runs on-line with good performance. The 3D positions of the speakers' heads are obtained from the stereo image pair, and interaural time difference (ITD) values are extracted from the binaural signals. These features are then fused in a probabilistic manner in order to compute, over time, the probability of each person's speaking activity.

The approach exhibits a number of novelties with respect to previous work addressing audio-visual fusion for speaker detection: (i) visual features are obtained from a stereoscopic setup and thus represented in 3D, (ii) auditory features are obtained from only two microphones, while most of previous work uses an array of microphones, (iii) the software is reusable with other robot sensor architectures, due to the flexibility of the underlying middleware layer, and (iv) good on-line performance in a complex environment, e.g., echoic rooms, simultaneous auditory sources, background noise, uncontrolled lighting, cluttered scenes, etc.

III. AN AUDIO-VISUAL FUSION MODEL

The overall goal is to retrieve the audio-visual (AV) state of the speakers in front of the robot. That is, the number of speakers as well as their positions and their speaking state. In order to reach this goal, we adopted the framework proposed in [START_REF] Alameda-Pineda | Finding audio-visual events in informal social gatherings[END_REF]. Based on a multimodal Gaussian mixture model (mGMM), this method is able to detect and localize audiovisual events from auditory and visual observations. We chose this framework because it is able to account for several issues: (i) the observation-to-speaker assignment problem, (ii) observation noise and outliers, (iii) the possibility to weight the relevance of the two modalities, (iv) a generative

{π n } N +1 n=1 Z {µ n , σ 2 n } N n=1 X
Fig. 2: Graphical model generating the audio-visual observations. The hidden variable Z follows a multinomial distribution with parameters π 1 , . . . , π N +1 . The audio-visual observations X follow the law described by the probability density function in Equation 2.

formulation linking the audio and visual observation spaces, and (v) the possibility to deal with a varying number of speakers through a principled model selection method.

In a first stage, the low-level auditory and visual features are extracted. While the former correspond to the interaural time differences (ITDs), the latter correspond to interest points in image regions related to motion which are further reconstructed in the 3D space using a stereo algorithm. These 3D points will be referred to as the visual features.

The following direct sound propagation model:

ITD(S) = S -M 1 -S -M 2 ν , (1)
is assumed. In this equation S corresponds to the sound source positions in the 3D space, e.g., a speaker, M 1 and M 2 are the 3D coordinates of the microphones in some robot-centered frame, and ν denotes the sound speed. Equation (1) maps 3D points onto the 1D space of ITD observations. The key aspect of our generative audio-visual model [START_REF] Khalidov | Conjugate mixture models for clustering multimodal data[END_REF], [START_REF] Alameda-Pineda | Finding audio-visual events in informal social gatherings[END_REF] is that (1) can be used to map 3D points (visual features) onto the ITD space associated with two microphones, on the premise that the cameras are aligned with the microphones [START_REF] Khalidov | Calibration of a binocularbinaural sensor using a moving audio-visual target[END_REF]. Hence the fusion between binaural observations and binocular observations is achieved in 1-D.

The underlying multimodal GMM (mGMM) is a onedimensional mixture of Gaussians. Each mixture component is associated with an audio-visual object centered at µ n and with variance σ 2 n . This mixture has the following probability density function:

prob(x ; Θ) = N n=1 π n N (x ; µ n , σ 2 n ) + π N +1 U(x), ( 2
)
where N is the number of components, i.e., audio-visual objects, π n is the weight of the n th component, N (x ; µ n , σ 2 n ) is the value of the Gaussian distribution at x, U is the value of the uniform distribution accounting for outliers, and Θ = {π n , µ n , σ 2 n }. In this equation, x stands for a realization of the random variable X, shown in the corresponding graphical model on Fig. 2, that could be either an auditory observation, i.e., an ITD value or an observed 3D point, i.e., a visual feature, mapped with [START_REF]message filters/approximatetime[END_REF]. Notice that both Θ and the hidden variable Z (modeling the observation-to-object assignments) need to be estimated. This is done using an Expectation-Maximization (EM) algorithm, derived from the probabilistic graphical model. Notice that with this formulation the number of AV objects N can be estimated from the observed data by maximizing a Bayesian information criterion (BIC) score [START_REF] Alameda-Pineda | Finding audio-visual events in informal social gatherings[END_REF]. However, this implies to run the EM algorithm several times with different values of N , which is prohibitive in the case of an on-line implementation. From a practical point of view the problem of estimating N can be overcome by replacing the 3D visual points with 3D faces as described below.

A. Visual Processing

The initial implementations of the nGMM EM algorithm was using 3D points [START_REF] Khalidov | Conjugate mixture models for clustering multimodal data[END_REF], [START_REF] Alameda-Pineda | Finding audio-visual events in informal social gatherings[END_REF] as just described. Alternatively, one can replace 3D points with 3D faces, more precisely with 3D face centers which are fair approximations of 3D mouth positions, i.e., the 3D acoustic emitters. In practice we start by detecting faces in images using [START_REF] Šochman | Waldboost -learning for time constrained sequential detection[END_REF]. Face centers are then detected in the left image of the stereo camera pair. Each left-image face center is then correlated with the right image along an epipolar line in order to obtain a stereo correspondence between the face center and a associated point along the epipolar line in the right image. This allows to reconstruct a 3D point, S n , that can be viewed as 3D face center. See [START_REF] Hansard | Cyclopean geometry of binocular vision[END_REF] for more details. The use of faces drastically simplifies the complexity of the approach because a single semantically-meaningful face center replaces a cloud of points associated with a, possibly moving, 3D object. Initial means can be easily obtained from (1), i.e., µ n = ITD(S n ) while N , the number of AV objects can be easily estimated using the face detector [START_REF] Šochman | Waldboost -learning for time constrained sequential detection[END_REF].

B. Auditory Processing

As already mentioned we use ITDs, i.e., the time delay between the signals received at the left and right microphones. Notice that, due the symmetric nature of the ITD function, there is a front/back ambiguity, which is however slightly attenuated by the transfer function of the robot head. There are several methods to estimate ITDs (see [START_REF] Chan | Time-of-arrival based localization under NLOS conditions[END_REF] for a review); We chose the cross-correlation method, since it optimizes a trade-off between performance and complexity. ITD values are obtained in real-time by computing the crosscorrelation function between the left and right perceived signals during an integration time window of length W , expressed in number of time samples, or frames. The time delay τ corresponding to the maximum of the cross-correlation function in the current integration window is computed as follows:

τ = 1 F s argmax d∈[-d M ,d M ] W t=1 l(t)r(t + d) (3)
where l and r are the left and right audio signals, F s is the sampling frequency and d M denotes the maximum possible delay between microphones, i.e.,

d M = M 1 -M 2 /ν.
The time window W is a trade-off between reliability and significance. On one hand, a high W value implies more reliable ITD values, since the effect on the local maxima of the cross-correlation function is reduced. On the other hand, a small W value speeds up the computation. The parameter f denotes the shift of the sliding window used to compute the ITD. In order to extract one ITD value, two conditions need to be satisfied. First, there should be enough samples available within the integration window W . Second, the mean energy of the signals in the integration window should be higher than a given threshold E A . In this way, we avoid to compute ITD values when the audio stream contains nothing but noise. Notice the method does not assume that the perceived sound signals are associated with some semantic i.e., speech, pulse-resonance sounds, etc.

IV. SYSTEM CALIBRATION

The audio-visual fusion model outlined above, and 1 in particular, implies that the visual and auditory observations are computed in a common reference frame. This allows visual data to be aligned with auditory data. In practice it means that the cameras' extrinsic calibration parameters (position and orientation) and the microphones' positions are expressed in a common reference frame. Extrinsic camera calibration is performed using the state-of-the-art algorithm of [START_REF] Bouguet | Camera calibration toolbox for matlab[END_REF].

Audio-visual calibration can be achieved using [START_REF]message filters/approximatetime[END_REF]. A sound-source is placed in a known position S while M 1 and M 2 are unknown and hence must be estimated. The method of [START_REF] Khalidov | Calibration of a binocularbinaural sensor using a moving audio-visual target[END_REF] (i) uses an audio-visual target (a loud-speaker emitting white noise coupled with a small red-light bulb) to precisely position the sound source in the camera-pair reference frame, and (ii) estimates the unknown parameters M 1 and M 2 by considering several target positions and by solving a non-linear system of equations of the form of (1). This calibration procedure does not take into account the fact that the microphones are plugged into the robot head, as already mentioned above. To account for head effects we introduce two corrective parameters, α and β, to form of an affine transformation.

ITD AD (S) = α S -M 1 -S -M 2 ν + β, (4)
These parameters are estimated using the same audio-visual target mentioned above. The audio-visual target is freely moved in front of the robot thus following a zigzag-like trajectory. The use of white noise greatly facilitate the task of cross-correlation, i.e., there is single shapr peak, and hence, makes the ITD computation extremely reliable. The reverberant components are suppressed by the direct component of the long lasting white noise signal. However, it is possible to set up the experimental conditions such as to reduce the effects of reverberation, e.g., the room size is much larger than the target-to-robot distance. If the microphone positions are estimated in advance, the estimation of α and β can be carried out via a linear least-square estimator derived from 4. Fig. 3 shows the extracted ITDs (red-circle), mapped 3D face centers before the adjustment (blue), i.e., using (1) and after the adjustment (green), i.e., using (4). Clearly, the affine correction model allows a better alignment between the visual and auditory data.

V. SYSTEM ARCHITECTURE

We implemented the method described above as a number of components connected within the framework of robotics services bus (RSB) [START_REF] Wienke | A middleware for collaborative research in experimental robotics[END_REF]. RSB is a platform-independent event-driven middleware specifically designed for the needs of distributed robotic applications. It is based on a logically unified bus which can span over several transport mechanisms like network or in-process communication. The bus is hierarchically structured using scopes on which events can be published with a common root scope. Through the unified bus, full introspection of the event flow between all components is easily possible. Consequently, several tools exist which can record the event flow and replay it later, so that application development can largely be done without a running robot. RSB events are automatically equipped with several timestamps, which provide for introspection and synchronization abilities. Because of these reasons RSB was chosen instead of NAO's native framework NAOqi and we could implement and test our algorithms remotely without performance and deployment restrictions imposed by the robot platform. Moreover, the resulting implementation can be reused for other robots.

One tool available in the RSB ecosystem is an event synchronizer, which synchronizes events based on the attached timestamps with the aim to free application developers from such a generic task. However, several possibilities of how to synchronize events exist and need to be chosen based on the intended application scenario. For this reason, the synchronizer implements several strategies, each of them synchronizing events from several scopes into a resulting compound event containing a set of events from the original scopes. We used two strategies for the implementation. The ApproximateTime strategy is based on the algorithm available in [START_REF]message filters/approximatetime[END_REF] and outputs sets of events containing exactly one event from each scope. The algorithm tries to minimize the time between the earliest and the latest event in each set and hence well-suited to synchronize events which originate from the same source (in the world) but suffered from perception or processing delays in a way that they have non-equal timestamps. The second algorithm, TimeFrame, declares one scope as the primary event source and for each event received here, all events received on other scopes are attached that lie in a specific time frame around the time stamp of the source event.

ApproximateTime is used in our case to synchronize the results from the left and right camera as frames in general form matching entities but due to independent grabbing of both cameras have slightly different time stamps. Results from the stereo matching process are synchronized with ITD values using the TimeFrame strategy because the integration time for generating ITD values is much smaller than for a vision frame and hence multiple ITD values belong to a single vision result.

A. Modular Structure

The algorithm is divided into four components which are described in the pipeline shown in Fig. 5. Each component is represented with a different color and white is reserved for modules provided by the RSB middleware as synchronization functionalities. Red relates to the auditory processing (see section III-B). Green to the visual features (see section III-A). Purple is for the audio-visual fusion method (described in section III) and finally blue is the visualization tool (described at the end of this section).

In detail, the visual part has five different modules. Left video and Right video stream the images received from left and right cameras. The Left face detection module extracts the faces from the left image. These are then synchronized with the right image in Face-image Synchronization, using the ApproximateTime strategy. The 3D Head Position module computes the 3D head (or face) centers.

The auditory component consists of three modules. Interleaved audio samples coming from the two microphones are streamed by the Audio module. These are de-interleaved by Sound formatting and stored into two circular buffers; for the left and right microphone's signals. Finally, the module called ITD extraction is in charge of compute the ITD values. Both visual and auditory features flow until the Audiovisual synchronization module; the TimeFrame strategy is used here to find the ITD values coming from the audio pipeline associated to the 3D head positions coming from the visual processing. These synchronized events feed the Audio-visual fusion module, which is in charge of estimating the emitting sound probabilities p n .

Finally, we developed the module Visualization, in order to get a better insight of the proposed algorithm. A snapshot of this visualization tool can be seen in Fig. 4. The visualization plot consists of three parts. The top-left part displays a bounding box around each detected face overlaid onto the original left image. In addition to the bounding boxes, a solid circle is plot on a face whenever it emits a sound. The intensity encodes the emitting sound probability, the higher it is, the darker the circle. The top-right part, framed in green, is a bird-view of the 3D reconstructed scene, in which the detected 3D faces are shown with circles. The bottom-left part, with a red frame, represents the ITD space in which both the mapped face centers (ellipses) and the histogram of ITD values are plot.

B. Implementation Details

Several considerations need to be done regarding the details of the algorithm implementation in order to guarantee the repeatability of the experiments. When computing the ITD values, a few parameters need to be set. As explained previously (see section III-B), there is a trade off when setting the integration window W and the frame shift f . A good compromise between low computational load, high rate, and reliability of ITD values was found for W = 150 ms and f = 20 ms. Finally, we set the activity threshold to E A = 0.001. Notice that this parameter could be controlled by a higher level module which would learn the characteristics of the scene and infer the level of background noise. When computing the probabilities p n , the variances of the mGMM are set to σ 2 n = 10 -9 , we found this value big enough to take into account the noise in the ITD values and small enough to discriminate speakers that are close to each other. The entire pipeline was running on a laptop with an i7 processor at 2.5 GHz. Last but not least, we used a new audio-visual head for NAO (see Fig 6). The novelty of this head is the synchronized camera pair delivering a VGA stereoscopic video flow.

VI. EXPERIMENTAL VALIDATION

To validate our algorithm we performed a set of experiments with five different scenarios. The scenarios were recorded in a room around 5 × 5 meters, and designed to test the algorithm in different conditions in order to identify the limitations of the proposed approach. Each scenario is composed by several sequences in which people count from one up to sixteen. Except for the first scenario, composed by one sequence due to its simplicity, the rest of scenarios were recorded several times. Moreover, a video is recorded to show the different scenarios and have a visual validation of the results.

In scenario S1, only one person is in the room sitting in front of the robot and counting. In the rest of the scenarios (S2-S5) three persons are in the room. People are not always in the field of view (FoV) of the cameras and sometimes they move. In scenario S2 three persons are sitting and counting alternatively one after the other. The configuration Fig. 6: Within this work we used a new audiovisual head that is composed of a synchronized camera pair and two microphones. This "orange" head replaces the former "blue" head and is fully interfaced by the RSB middleware described in this paper. of scenario S3 is similar to the one of S2, but one person is standing instead of sitting. These two scenarios are useful to determine the precision of the ITDs and experimentally see if the difference of height (elevation) affects the quality of the extracted ITDs. The scenario S4 is different from S2 and S3 because one of the actors is outside the FoV. This scenario is used to test if people speaking outside the FoV affect the performance of the algorithm. In the last scenario (S5) the three people are in the FoV, but they count and speak independently of the other actors. Furthermore, one of them is moving while speaking. With S5, we aim to test the robustness of the method to dynamic scenes.

In Fig. 7 we show several snapshots of our visualization tool. These frames are selected from the different scenarios aiming to show both the successes and the failures of the proposed system. Fig. 7a shows an example of perfect alignment between the ITDs and the mapped face, leading to a high speaking probability. A similar situation is presented in Fig. 7b, in which among the three people, only one speaks. A failure of the ITD extractor is shown in Fig. 7c, where the actor in the left is speaking, but no ITDs are extracted. In Fig. 7d we can see how the face detector does not work correctly: one face is missing because the actor is too far away and the other's face is partially occluded. Fig. 7e shows a snapshot of an AV-fusion failure, in which the extracted ITDs are not significant enough to set a high speaking probability. The Fig. 7f, Fig. 7g and Fig. 7h show the effect of reverberations. While in Fig. 7h we see that the reverberations lead to the wrong conclusion that the actor on the right is speaking, we also see that the statistical framework is able to handle reverberations (Fig. 7f and Fig. 7g), hence demonstrating the robustness of the proposed approach. for the five scenarios. The columns represent, in order: the amount of correct detections (CD), the amount of false positives (FP), the amount of false negatives (FN) and the total number of counts (Total).

The scenarios are manually annotated such that we get the ground truth. In order to systematically evaluate the proposed system we adopted an overlap-based strategy. The ground truth of actor n is split in speaking intervals I k n indexed by k and silent intervals J l n indexed by l. For clarity purposes let us denote by p n (t) the detected speaking state of actor n at time t. For each of the speaking intervals I k n we compute c k n = t∈I k n p n (t)/|I k n |. If c k n ≥ 0.5 we count one correct detection, otherwise we count one false negative. We also compute ck n = t∈J l n p n (t)/|J l n |. In case cl n ≥ 0.5 we count one false positive. In summary, if the speaker is detected during more than half of the speaking time, we count on correct detection (CD), otherwise a false negative (FN). And if it is detected more than half of the speaking time, we count a false positive (FP).

Table I shows the results obtained with this evaluation strategy on the presented scenarios. First of all we notice the small amount of false negatives: the system misses very few speakers. A part from the first scenario (easy conditions), we observe some false positives. These false positives are due to reverberations. Indeed, we notice how the percentage of FP is severe in S5. This is due to the fact that high reverberant sounds (like hand claps) are also present in the audio stream of this scenario. We believe that an ITD extraction method more robust to reverberations will lead to more reliable ITD values, which in turn will lead to a better speaker detector. It is also worth to notice that actors in different elevations and non-visible actors do not affect the performance of the proposed system, since the results obtained in scenarios S2 to S3 are comparable.

VII. CONCLUSIONS AND FUTURE WORK

We presented a system targeting speaker detection working on the humanoid robot NAO in regular indoor environments. Implemented on top of a platform-independent middleware, the system processes the audio-visual data flow from two microphones and two cameras at a rate of 17 Hz. We proposed a statistical model which captures outliers from the perception processes. The method runs in normal echoic rooms with just two microphones mounted inside the head of a companion robot with noisy fans. We demonstrated good performance on different indoor scenarios involving several actors, moving actors and non-visible actors. This works contributes to a better understanding of the audiovisual scene using a robocentric set of sensors mounted in an autonomous platform, such as NAO, under the constraints of an on-line application.

Fig. 3 :

 3 Fig. 3: The effect of the corrective parameters α and β onto audio-visual calibration. ITD values estimated as peaks of the auto-correlation function are shown with red circles. 3D face centers are mapped onto the ITD space without the corrective parameters (shown in blue) and once these parameters have been estimated (shown in green).

Fig. 4 :

 4 Fig. 4: Snapshot of the visualization tool. Top-left (blueframed): The original left image overlaid with one bounding box per detected face. In addition, an intensity-coded circle appears when the face emits a sound. The darker the circle, the higher the speaking probability. Top-right (greenframed): A bird-view of the 3D scene, in which each circle corresponds to a detected head. Bottom-left (red-framed): The ITD space. The projected faces are displayed with an ellipse while extracted ITD values are shown as bars in a histogram.

Fig. 5 :

 5 Fig. 5: Modular structure of the processing pipeline. There are five types of modules: streaming & synchronization (white), visual processing (green), auditory processing (red), audio-visual fusion (purple) and the visualization module (blue).

Fig. 7 :

 7 Fig. 7: Snapshots of the visualization tool. Frames are selected among the five scenarios such as to show both the method's strengths and weaknesses. (a) Good results on S1. (b) Good results on S2, three people. (c) The ITD extractor does not work correctly, thus missing the speaker. (d) Misses of the face detection module. (e) The audio-visual fusion fails to set a high probability to the current speaker. (f,g) The audio-visual fusion model is able to handle reverberations. (h) The reverberations are too close to the mapped head, leading to a wrong decision.

TABLE I :

 I Quantitative evaluation of the proposed approach

http://www.aldebaran-robotics.com

It is worth noticing that the module limiting the performance of the system is the ITD extraction, due to the room reverberations. We will work on making this module more robust to this kind of interferences. Moreover, audio-visual tracking capabilities are also a desirable property for any robot, since they provide for temporal coherence of the scene. In a more developed stage, it would be desirable that NAO is able to choose regions of interest, so that it could perform active learning, and enhance its audio-visual skills.


Jordi Sanchez-Riera
email: jordi.sanchez-riera@inria.fr

Jan Čech
email: jan.cech@inria.fr

Radu Horaud
email: radu.horaud@inria.fr

Robust Spatiotemporal Stereo for Dynamic Scenes

Stereo matching is a challenging problem, especially in the presence of noise or of weakly textured objects. Using temporal information in a binocular video sequence to increase the discriminability for matching has been introduced in the recent past, but all the proposed methods assume either constant disparity over time, or small object motions, which is not always true. We introduce a novel stereo algorithm that exploits temporal information by robustly aggregating a similarity statistic over time, in order to improve the matching accuracy for weak data, while preserving regions undergoing large motions without introducing artifacts.

Introduction

The use of a stereo video presents advantages over single-frame stereo. The extra information may help the matching process to improve accuracy and to enforce temporal consistency of disparity maps in the case of weakly textured regions or when the matching is ambiguous. We see three main approaches in the literature.

The first group computes the scene flow. These methods simultaneously estimate depth and motion. The is a formulation of the coupled estimation problems of disparity (between a stereo pair) and optical flow (between consecutive frames) which mutually constrain each other. This task is traditionally solved by variational [START_REF] Basha | Multi-view scene flow estimation: A view centered variational approach[END_REF], MRF [START_REF] Liu | Disparity estimation in stereo sequences using scene flow[END_REF], or seed-growing [START_REF] Cech | Scene flow estimation by growing correspondence seeds[END_REF] methods.

In the second group, there are methods which rely on independent motion estimates to improve the stereo matching. In [START_REF] Bleyer | Temporally consistent disparity maps from uncalibrated stereo videos[END_REF], the disparity maps are filtered by a median filter along pixel trajectories obtained by an external optical flow module. Independent optical flow is also used in [START_REF] Zhu | Spatial-temporal fusion for high accuracy depth maps using dynamic MRFs[END_REF], where the authors propose an MRF framework with an extra term which penalizes discrepancies in photoconsistency of the (optical flow related) neighbourhood in the current-, previous-, and next frames. Similar MRF formulation [START_REF] Larsen | Temporally consistent reconstruction from multiple video streams using enhanced belief propagation[END_REF] additionally disconnects the edges to prevent over-smoothing in case of large motion and failure of the optical flow estimates.

The third group is composed of methods of spatiotemporal stereo that do not estimate the motion explicitly, but exploit a local spatiotemporal neighbourhood of pixels to increase discriminability of the similarity statistics. Paper [START_REF] Davis | Spacetime stereo: A unifying framework for depth from triangulation[END_REF] projects an artificial pattern varying over time onto the static scene and temporally aggregate the statistic. The similarity statistic (based on bilateral filtering) is temporally aggregated also in [START_REF] Richardt | Real-time spatiotemporal stereo matching using the dual-cross-bilateral grid[END_REF], such that adjacent frames are weighted by a Gaussian kernel to cope with a small motion. In [START_REF] Zhang | Spacetime stereo: Shape recovery for dynamic scenes[END_REF] the authors study how spatiotemoral windows are deformed due to surface slant and motion and propose an optimization framework to find the distortion parameters and invariant similarity statistic. Alternatively, the same insensitivity is achieved in [START_REF] Sizintsev | Spatiotemporal stereo via spatiotemporal quadric element (stequel) matching[END_REF][START_REF] Sizintsev | Spatiotemporal oriented energies for spacetime stereo[END_REF] by representing the image using Gabor filter responses and the similarity statistic is computed in a closed form. However, all these methods assume that the disparity between frames changes only slowly. In reality this assumption is not valid near object boundaries or in the presence of rapidly moving objects, which causes serious artifacts.

The main contribution of this paper is a method of spatiotemporal stereo which benefits from aggregating the similarity statistic over a time window. Unlike previous work on spatiotemporal stereo algorithms, the proposed method is robust to abrupt temporal changes in disparity due to large motions. The main idea of our algorithm is to automatically detect image regions corresponding to such changes, such that the aggregation of the similarity statistic over the time window is disconnected for these regions. The algorithm is implemented within the efficient seed growing procedure [START_REF] Cech | Efficient sequential correspondence selection by cosegmentation[END_REF].

Method description

Let us have two synchronized and epipolar-rectified video streams I l (x, y, t) and I r (x, y, t). Variables x and y index respective horizontal and vertical image coordinates in pixels, t indexes the time, i.e. a frame of the streams. The streams are related by a disparity function d(x, y, t) which assigns the correspondences between pixels in the left and right image I l (x, y, t) ≈ Ir(x + d(x, y, t), y, t).

(

A matching algorithm must measure a certain similarity statistic between potentially corresponding pixels to establish the correspondences. The simplest statistic is a difference of pixel intensities, which is however ambiguous. More discriminable statistics use a small neighbourhood (a window) around the potential correspondence. Then, these algorithms locally approximate the disparity function. For instance, paper [START_REF] Zhang | Spacetime stereo: Shape recovery for dynamic scenes[END_REF] uses a linear approximation. In a small spatiotemporal neighbourhood N around location (x 0 , y 0 , t 0 ), e.g. a 3D window of 5 × 5 pixels over 3 frames, the disparity function

is d(x, y, t) ≈ d(d i , d 0 , d 1 , d 2 , d t ) = d i + d 0 + d 1 (x - x 0 ) + d 2 (y -y 0 ) + d t (t -t 0 ).
Then they use an optimized statistic to measure a photometric consistency of the potential correspondence for candidate disparities

d i TSSD(x0, y0, t0, di) = min d 0 ,d 1 ,d 2 ,d t X (x,y,t)∈N " I l (x, y, t)-Ir(x + d(di, d0, d1, d2, dt), y, t) " 2 (2)
to compensate the distortion which occurs due to subpixel displacement d 0 , surface slant d 1 , d 2 , and temporal disparity change d t .

However, there are several sources of errors in this approach: (i) Tendency to get stuck in a local extrema; (ii) Not a significant gain in discriminability 1 over the case where d 0 = d 1 = d 2 = d t = 0, since the statistic is improved by the optimization for both correct and incorrect matches; (iii) When the assumption on the linearity of the disparity function within the local spatiotemporal neighbourhood is violated (e.g. abrupt change in disparity), the method fails dramatically.

Therefore we adopted a different approach in the aggregation. As an elementary similarity statistic, we use Moravec's normalized cross correlation [START_REF] Moravec | Towards automatic visual obstacle avoidance[END_REF],

NCC(x0, y0, t0, di) = 2 cov `Wl (x0, y0, t0), Wr(x0 + di, y0, t0) var `Wl (x0, y0, t0) ´+ var `Wr(x0 + di, y0, t0) ´+ ǫ , (3)
where W l (x 0 , y 0 , t 0 ) is a spatial window (a subimage) of N × N pixels centered at position (x 0 , y 0 ) of the frame t 0 of the left stream. Similarly W r , and ǫ is a machine epsilon to prevent instability of the statistic in case of low intensity variance. The statistic has consequently low response in textureless regions [START_REF] Cech | Efficient sequential correspondence selection by cosegmentation[END_REF].

Then the NCC statistic is aggregated over a symmetric time window of 2T + 1 frames, such that

TNCC(x0, y0, t0, di) = 1 2T + 1 t 0 +T X t=t 0 -T NCC(x0, y0, t, di). (4)
Apparently, the TNCC statistic is decayed when the disparity changes significantly within the temporal window. Notice that the motion in general is not harmful, but the motion changing the disparity is. A typical distribution of NCC(t) statistic over the time t of the window for a correct match (x 0 , y 0 , t = 0, d i ) is the following: If the disparity is constant over time, all perframe correlations for t = {-T, . . . , T } are high. If the disparity changes slowly, the correlation is slightly lower more faraway from the central frame. However, when the disparity changes rapidly, the correlation off the central frame drops quickly, since the other correlations measure a photometric consistency of locations which are not corresponding any more.

On the other hand, a potential mismatch (i.e. wrong correspondence) has the distribution of per-frame correlations over the time window such that the correlations are low, but due to random fluctuations or texture selfsimilarity there may be high responses for any frame of the temporal window. The temporal aggregation in (4) averages out these excesses and decreases their correlations and hereby increases the discriminability.

However, it is important to detect phenomena corresponding to large changes in disparity and in these cases to use the central correlation only without any aggregation which would cause artifacts. Therefore, we propose a robust temporal normalized cross correlation

RTNCC(x0, y0, t0, di) = ( NCC(x0, y0, t0, di) if `NCC(t0)-NCC(t0 ±1) ´≥ α, TNCC(x0, y0, t0, di) otherwise. ( 5
)
This means that RTNCC uses the correlation (3) of the central frame, if it is higher than correlations of adjacent frames by threshold α. Otherwise, RTNCC uses the average correlation TNCC over the entire temporal window (4). For simplicity of notation, we omitted all other indexes x 0 , y 0 , d i in the condition in [START_REF] Davis | Spacetime stereo: A unifying framework for depth from triangulation[END_REF].

Matching algorithm. To establish the matching between stereo images, the proposed RTNCC statistic was integrated in a seed growing procedure [START_REF] Cech | Efficient sequential correspondence selection by cosegmentation[END_REF]. This algorithm uses seed correspondences obtained by match- Besides low computational complexity, the advantage of the algorithm in our context is the input of seeds. Namely, we observed the condition in (5) of RTNCC is reliable in textured regions only. Nevertheless, the seed correspondences are points with the Harris property and for them the decision works well. Therefore, we propose to take this decision for the initial seeds only. Each seed then propagates a flag indicating whether the aggregation in RTNCC is used or not and this flag is inherited by its 'offspring' seeds in the growing process.

Experiments

We performed a set of experiments to demonstrate that the proposed algorithm can cope with weak or ambiguous data and images corrupted by noise, without introducing artifacts of smoothing boundaries of rapidly moving objects. We compare the proposed method (RT-NCC) with two baseline instances of the growing algorithm: (i) the algorithm which uses the spatial neighbourhood only for matching (NCC), and (ii) the algorithm which trivially uses the spatio-temporal neighbourhood, such that all per-frame correlations are averaged (TNCC). The other comparisons are with two state-of-the art spatio-temporal methods: (i) temporal SSD optimization [START_REF] Zhang | Spacetime stereo: Shape recovery for dynamic scenes[END_REF] integrated in the growing algorithm (TSSD), and (ii) the stequel matching algorithm [START_REF] Sizintsev | Spatiotemporal stereo via spatiotemporal quadric element (stequel) matching[END_REF] (Sizintsev09).

For all experiments, we used 5 × 5 pixel windows as the spatial neighbourhood of all statistics, parameter α in RTNCC [START_REF] Davis | Spacetime stereo: A unifying framework for depth from triangulation[END_REF] was empirically set to 0.8. For the short synthetic sequence, we set temporal window half-size to T = 2, while for the real one it was set to T = 7. Ground-truth experiment. To quantitatively compare the different algorithms, we tested on a stereo sequence with ground-truth disparity maps used in [START_REF] Cech | Scene flow estimation by growing correspondence seeds[END_REF], Fig. 2(a). It consists of three objects: a background plane, a sphere, and a thin bar. The plane and sphere move slowly, while the thin bar moves rapidly (about 30 pixels per frame) from right to left crossing the entire scene. It is textured randomly with a white noise.

We measured ratio of correctly matched pixels in non-occluded regions, i.e. number of all pixels without mismatches (error ≥ 1 pixel) and unmatched pixels divided by the total number of pixels. The input images were perturbed with zero mean additive Gaussian noise with successively increasing standard deviation σ. The noise has equal variance as the signal for σ = 1.

In Fig. 1, we can see the algorithm using NCC performs very well without noise. The texture is optimal and hence the correlation is very high and unambiguous. However, it degrades with noise, Fig. 2(c), producing mismatches as small spatial only image windows do not correlate well. The TNCC degradates slowly with increasing level of noise, however the ratio of correct matches is lower since it tends to completely miss the rapidly moving bar, Fig. 2(d). The temporal aggregation helps to filter out the noise in slowly moving regions, but the aggregation is harmful for the bar where the disparity changes abruptly over time, since for this region TNCC of the false background wins over TNCC of the true bar. Similarly, the other two methods Siz-intsev09 and TSSD perform well filtering the noise but both of them have serious problems with the rapidly moving bar where the disparity changes abruptly over time, Fig. 2(e), 2(f). The proposed RTNCC performs the best. It is as good as the NCC for low noise, and it is always superior to other methods with increasing levels of noise. It has fewer mismatches in slowly moving regions, but at the same time it preserves the rapidly moving bar without artifacts, Fig. 2(g). It correctly aggregates over time using TNCC in regions where it helps, and for other regions it uses the spatial statistic NCC. The map in Fig. 2(h) shows which case in (5) was used in results of RTNCC. Pixels matched using the temporal aggregation are indicated by gray colour, pixels matched by spatial statistic by black colour. We can see, it correctly used the spatial statistic NCC for the region of the bar, while for other pixels, it correctly used the temporal aggregation TNCC.

Real outdoor scene. To show the validity of the proposed algorithm on a real outdoor scene we tested under the DAGM Exposure Changes dataset2 (DAGM). The stereo camera is in a car driving in a highway in difficult lighting conditions, sudden exposure changes. Cars going in the opposite lane moves very fast, see Fig. 3(a). We show results for the frame, where the car is passing under the bridge. The texture of the road almost disappears. It causes the spatial statistic NCC to fail, producing many mismatches, as in Fig. 3(c). The spatiotemporal version TNCC works better, Fig. 3(d). Much information is retained due to the temporal aggregation. Notice that the disparity of the road remains constant over time and this is also the case of the car going in the same direction, since the distance to it is more or less constant. However, the problem is, that the car going in the opposite direction (in red circle), whose relative velocity is very high, is missed by the TNCC. This is the same effect as the case of the rapidly moving bar in Fig. 2(d). TSSD has similar difficulties there, Fig. 3(f). Surprisingly, algorithm Sizintsev09 has severe problems with all rapidly moving pixels in the scene,

Figure 1 .

 1 Figure 1. Quantitative evaluation.

  Sizintsev09 (g) RTNCC (h) RTNCC agg.

Figure 2 .

 2 Figure 2. Synthetic dataset of [4]. Frame 6, noise level σ = 0.5. Disparity maps.

Figure 3 .

 3 Figure 3. DAGM Exposure Challenge dataset. Disparity maps.

The discriminability of the similarity statistic is proportional to a probability that the statistic has better response for the true correspondence than for the incorrect ones.

www.dagm2011.org/adverse-vision-conditions-challenge.html including those where the disparity remains constant. It produces large artifacts in regions near the camera. The proposed RTNCC works well, Fig.

3(e). It is significantly superior to NCC and all objects, including the car in the opposite direction, are preserved.

Conclusions

We presented a spatiotemporal correlation statistic that increases the discriminability by aggregating over time and hereby produces higher quality matching results. We showed the proposed method is robust to a rapid motion in the scene, which is a situation where the state-of-the-art algorithms are prone to artifacts.

Acknowledgements. This research was supported by EC project FP7-ICT-247525-HUMAVIPS.