Nicolas Bonneel

George Drettakis

Nicolas Tsingos

Isabelle Viaud-Delmon

Doug James

Fast Modal Sounds with Scalable Frequency-Domain Synthesis

Keywords: I.6.8 [Simulation and Modeling]: Types of Simulation-Animation, I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically based modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Virtual Reality Sound synthesis, modal synthesis, real-time audio rendering, physically based animation Amplitude frequency time Amplitude

Audio rendering of impact sounds, such as those caused by falling objects or explosion debris, adds realism to interactive 3D audiovisual applications, and can be convincingly achieved using modal sound synthesis. Unfortunately, mode-based computations can become prohibitively expensive when many objects, each with many modes, are impacted simultaneously. We introduce a fast sound synthesis approach, based on short-time Fourier Tranforms, that exploits the inherent sparsity of modal sounds in the frequency domain. For our test scenes, this "fast mode summation" can give speedups of 5-8 times compared to a time-domain solution, with slight degradation in quality. We discuss different reconstruction windows, affecting the quality of impact sound "attacks". Our Fourier-domain processing method allows us to introduce a scalable, real-time, audio processing pipeline for both recorded and modal sounds, with auditory masking and sound source clustering. To avoid abrupt computation peaks, such as during the simultaneous impacts of an explosion, we use crossmodal perception results on audiovisual synchrony to effect temporal scheduling. We also conducted a pilot perceptual user evaluation of our method. Our implementation results show that we can treat complex audiovisual scenes in real time with high quality.

Introduction

The rich content of today's interactive simulators and video games includes physical simulation, typically provided by efficient physics engines, and 3D sound rendering, which greatly increases our sense of presence in the virtual world [START_REF] Larsson | Better presence and performance in virtual environments by improved binaural sound rendering[END_REF]. Physical simulations are a major source of audio events: e.g., debris from explosions or impacts from collisions (Fig. 1). In recent work several methods have been proposed to physically simulate these audio events notably using modal synthesis [O' Brien et al. 2002;[START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF][START_REF] James | Precomputed acoustic transfer: Output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF]]. Such simulations result in a much richer virtual experience compared to simple recorded sounds due to the added variety and improved audio-visual coher- Frames of some of the test scenes our method can render in real time. Bottom: (Left) Time-domain modal synthesis requires summing all modes at every sample. (Right) Our frequency-domain modal synthesis exploits the inherent sparsity of the modes' discrete Fourier transforms to obtain lower costs per frame.

ence. The user's audiovisual experience in interactive 3D applications is greatly enhanced when large numbers of such audio events are simulated. Previous modal sound approaches perform time-domain synthesis [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]. Recent interactive methods progressively reduce computational load by reducing the number of modes in the simulation [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF][START_REF] Van Den Doel | Interactive simulation of complex audiovisual scenes[END_REF]]. Their computational overload however is still too high to handle environments with large numbers of impacts, especially given the limited budget allocated to sound processing, as is typically the case in game engines.

Interactive audiovisual applications also contain many recorded sounds. Recent advances in interactive 3D sound rendering use frequency-domain approaches, effecting perceptually validated progressive processing at the level of Fourier Transform coefficients [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]]. For faster interactive rendering, perceptually based auditory masking and sound-source clustering can be used [Tsin-gos et al. 2004;[START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF]]. These algorithms enable the use of high-quality effects such as Head-Related Transfer Function (HRTF) spatialization, but are limited to pre-recorded sounds. While the provision of a common perceptually based pipeline for both recorded and synthesized sounds would be beneficial, it is not directly obvious how modal synthesis can be efficiently adapted to benefit from such solutions. In the particular case of contact sounds, a frequency-domain representation of the signal must be computed on-the-fly, since events causing the sounds are triggered and controlled in real-time using the output of the physics engine.

Our solution to the above problems is based on a fundamental intuition: modal sounds have an inherently sparse representation in the frequency domain. We can thus perform frequency-domain modal synthesis by fast summation of a small number of Fourier coefficients (see Fig. 1). To do this, we introduce an efficient approximation to the short-time Fourier Transform (STFT) for modes. Compared to time-domain modal synthesis [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]], we observe 5-8 times speedup in our test scenes, with slight degradation in quality. Quality is further degraded for sounds with faster decays and high frequencies.

In addition to the inherent speed-up, we can integrate modal and recorded sounds into a common pipeline, with fine-grain scalable processing as well as auditory masking and sound clustering. To compute our STFT we use a constant exponential approximation; we also reconstruct with an Hann window to simplify integration in the full pipeline. However, these approximations reduce the quality of the onset of the impact sound or "attack". We propose a method to preserve the attacks (see Sect. 4.3), which can be directly used for modal sounds; using it in the combined pipeline is slightly more involved. We use the Hann window since it allows better reconstruction with a small number of FFT coefficients compared to a rectangular window. A rectangular window is better at preserving the attacks, but results in ringing. Our attack-preserving approach starts with a rectangular subwindow, followed by appropriate Hann subwindows for correct overlap with subsequent frames.

In contrast to typical usage of pre-recorded ambient sounds, physics-driven impact sounds often create peaks of computation load, for example the numerous impacts of debris just after an explosion. We exploit results from human perception to perform temporal scheduling, thus smoothing out these computational peaks. We also performed a perceptually based user evaluation both for quality and temporal scheduling. In summary, our paper has the following contributions:

• A fast frequency-domain modal synthesis algorithm, leveraging the sparsity of Fourier transforms of modal sounds.

• A full, perceptually based interactive audio rendering pipeline with scalable processing, auditory masking and sound source clustering, for both recorded and modal sounds.

• A temporal scheduling approach based on research in perception, which smooths out computational peaks due to the sudden occurrence of a large number of impact sounds.

• A pilot perceptual user study to evaluate our algorithms and pipeline.

We have implemented our complete pipeline; we present interactive rendering results in Sect. 7 and in the accompanying video for scenes such as those shown in Fig. 1 and Fig. 6.

Previous Work

In what follows we will use the term impact sound to designate a sound generated as a consequence of an event reported by the physics engine (impact, contact etc.); we assume that this sound is synthesized on-the-fly.

There is extensive literature on sound synthesis and on spatialization for virtual environments. We discuss here only a small selection of methods directly related to our work: contact sound synthesis in the context of virtual environments, audio rendering of complex soundscapes, and finally crossmodal perception, which we use to smooth computational peaks.

Contact sound synthesis

Modal synthesis excites the precomputed vibration modes of the objects by the contact force to synthesize the corresponding audio signal. In our context, this force is most often provided by a real-time physics engine. The acoustic response of an object to an impulse is given by:

s(t) = ∑ k a k e -α k t sin(ω k t), (1)
where s(t) is the time-domain representation of the signal, ω k is the angular frequency and α k is the decay rate of mode k; a k is the amplitude of the mode, which is calculated on the fly (see below), and may contain a radiation factor (e.g., see [START_REF] James | Precomputed acoustic transfer: Output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF]).

Eq. 1 can be efficiently implemented using a recursive formulation [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]] which makes modal synthesis attractive to represent contact sounds, both in terms of speed and memory. We define m k (t) as follows for notational convenience:

m k (t) = e -α k t sin(ω k t) (2)
The modal frequencies ω k 's and decay rates α k 's can be precomputed by simulating the mechanical behavior of each specific object independently. Such simulations can be done through finite elements methods [O' Brien et al. 2002], spring-mass system [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF], or even with analytic solutions for simple cases [van den [START_REF] Van Den Doel | Interactive simulation of complex audiovisual scenes[END_REF]. The only quantities which must be computed at run-time are the gains a k since they depend on the contact position on the objects, the applied force, and the listening position. Other methods based on measurements of the modes of certain objects are also possible [START_REF] Pai | Scanning physical interaction behavior of 3d objects[END_REF], resulting in the precomputation of the gains.

Audio rendering for complex soundscapes

There has been some work on modal sound synthesis for complex scenes. In [van den [START_REF] Van Den Doel | Interactive simulation of complex audiovisual scenes[END_REF]] a method is presented handling hundreds of impact sounds. Although their frequency masking approach was validated by a user study [START_REF] Van Den Doel | Measurements of perceptual quality of contact sound models[END_REF], the mode culling algorithm considers each mode independently, removing those below audible threshold. [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF] proposed a method based on mode pruning and sound sorting by mode amplitude; no perceptual validation of the approximation was presented however. For both, the granularity of progressive modal synthesis is the mode; in the examples they show, a few thousand modes are synthesized in real time.

For pre-recorded sampled sounds, [START_REF] Tsingos | Perceptual audio rendering of complex virtual environments[END_REF] have proposed an approach based on precomputed perceptual data which are used to cull, mask and prioritize sounds in realtime. This approach was later extended to a fully scalable processing pipeline that exploits the sparseness of the input audio signal in the Fourier domain to provide scalable or progressive rendering of complex mixtures of sounds [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF][START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF]]. They handle audio spatialization of several thousands of sound sources via clustering.

One drawback related to precomputed metadata is that sounds synthesized in real time, such as modal sounds, are not supported.

Cross-modal perceptual phenomena Audio-visual tolerance in asynchrony can be exploited for improved scheduling in audio rendering. The question of whether visual and auditory events are perceived as simultaneous has been extensively studied in neuroscience. Different physical and neural delays in the transmission of signals can result in "contamination" of temporal congruency. Therefore, the brain needs to compensate for temporal lags to recalibrate audiovisual simultaneity [START_REF] Fujisaki | Recalibration of audiovisual simultaneity[END_REF]. For this reason, it is difficult to establish a time window during which perception of synchrony is guaranteed, since it depends both on the nature of the event (moving or not) and its position in space (distance and direction) [START_REF] Alais | Synchronizing to real events: subjective audiovisual alignment scales with perceived auditory depth and speed of sound[END_REF]. Some studies report that delaying a sound may actually improve perception of synchrony with respect to visuals [START_REF] Begault | Auditory and non-auditory factors that potentially influence virtual acoustic imagery[END_REF]]. One study [START_REF] Guski | Audiovisual phenomenal causality[END_REF] (among others [START_REF] Sekuler | Sound alters visual motion perception[END_REF][START_REF] Sugita | Audiovisual perception: Implicit estimation of sound-arrival time[END_REF]), reports that a temporal window of 200 msec represents the tolerance of our perception for a sound event to be considered the consequence of the visual event. We will therefore adopt this value as a threshold for our temporal scheduling algorithm.

Our Approach

Overview The basic intuition behind our work is the fact that modal sounds have a sparse frequency domain representation. We will show some numerical evidence of this sparsity with examples, and then present our fast frequency-domain modal synthesis algorithm.

To achieve this we introduce our efficient STFT approximation for modes, based on singular distributions. We then discuss our full perceptual pipeline including scalable processing, auditory masking and sound source clustering. We introduce a fast energy estimator for modes, used both for masking and appropriate budget allocation. We next present our temporal scheduling approach which smooths out computation peaks due to abrupt increases in the number of impacts. After discussing our implementation and results, we describe our pilot perceptual user study, allowing us to evaluate the overall quality of our approximations and the perception of asynchrony. Analysis of our experimental results gives an indication of the perceptual validity of our approach.

Fourier-domain mode mixing Traditional time-domain modal synthesis computes Eq. 1 for each sample in time. For frequencydomain synthesis we use the discrete Fourier transform (FFT) of the signal (we show how to obtain this in Sect. 4.1). If we use a 1024-sample FFT we will obtain 512 complex coefficients or bins representing the signal of a given audio frame (since our signals are real-valued we will only consider positive frequencies). For each such frame, we add the coefficients of each mode in the frequency domain, and then perform an inverse FFT (see Fig. 1) once per frame, after all sounds have been added together. The inverse FFT represents a negligible overhead, with a cost of 0.023 msec using an unoptimized implementation [START_REF] Press | Numerical recipes in C: The art of scientific computing[END_REF]]. If the number of coefficients contributed by each mode is much less than 512, frequency-domain mixing will be more efficient than an equivalent time-domain approach. However, this will also result in lossy reconstruction, requiring overlapping frames to be blended to avoid possible artifacts in the time-domain. Such artifacts will be caused by discontinuities at frame boundaries resulting in very noticeable clicks. To avoid these artifacts, a window function is used, typically a Hann window. Numerous other options are available in standard signal processing literature [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]]. Our method shares some similarities with the work of [START_REF] Rodet | Spectral envelopes and inverse FFT synthesis[END_REF] which uses inverse FFTs for additive synthesis.

In what follows we assume that our audio frames overlap by a 50% factor, bringing-in and reconstructing only 512 new time-domain samples at each processing frame using a 1024-sample FFT. We implemented the Hann window as a product of two square roots of a Hann window, one in frequency to synthesize modes using few Fourier coefficients and the other in time to blend overlapping frames [Zölzer 2002]. At 44.1kHz, we thus process and reconstruct our signal using 512/44100 = 11msec-long frames.

Efficient Fourier-Domain Modal Synthesis

We provide some numerical evidence of our intuition, that most of the energy of a modal sound is restricted to a few FFT bins around the mode's frequency. We constructed a small test scene, containing 12 objects with different masses and material properties. The scene is shown in the accompanying video and in Fig. 6. We computed the energy with the signal reconstructed using all 512 bins, then progressively reconstruct with a small number of bins distributed symmetrically around the mode's frequency, and measured the error. We compute percent error averaged over all modes in the scene, for 1 bin (containing the mode's frequency), then 3 bins (i.e., together with the 2 neighboring bins on each side), then both these together with the 2 next bins, etc. Using a single bin, we have 52.7% error in the reconstructed energy; with 3 bins the error drops to 4.7% and with 5 bins the error is at 1.1%. We thus assume that bins are sorted by decreasing energy in this manner, which is useful for our scalable processing stage (see Sect. 5).

This property means that we should be able to reconstruct modal sounds by mixing a very small number of frequency bins, without significant numerical error; however, we need a way to compute the STFT of modes efficiently.

One possible alternative would be to precompute and store the FFTs of each mode and then weight them by their amplitude at runtime. However, this approach would suffer from an unacceptably high memory overhead and would thus be impractical. The STFT of a mode sampled at 44.1kHz requires the storage of 86 frames of 512 complex values, representing 352 Kbytes per mode per second. A typical scene of two thousand modes would thus require 688 Mb.

In what follows we use a formulation based on singular distributions or generalized functions [START_REF] Hormander | The Analysis of Linear Partial Differential Operators I[END_REF]], allowing us to develop an efficient approximation of the STFTs of modes.

A Fast Short-time FFT Approximation for Modes

We want to estimate the short-time Fourier transform over a given time-frame of a mode m(t) (Eq. 2), weighted by a windowing function that we will denote H(t) (e.g., a Hann window). We thus proceed to calculate the short time transform s(λ) where λ is the frequency, and t 0 is the offset of the window:

s(λ) = F λ { m(t + t 0 )H(t) }.
(3)

The Fourier transform F λ { f (t) } that we used corresponds to the definition:

F λ { f (t) } = ∞ -∞ f (t) e -i λ t dt (4)
Note that the product in the time domain corresponds to the convolution in the frequency domain (see Eq. 19 in the Appendix). We can use a polynomial expansion (Taylor series) of the exponential function:

e α(t+t 0 ) = e αt 0 ∞ ∑ n=0 c n (αt) n , (5)
where c n = 1/n!. Next, the expression for the Fourier transform of a power term is a distribution given by:

F λ { t n } = 2πi n δ (n) (λ), (6)
where δ is the Dirac distribution, and δ (n) its n'th derivative. From Eqs. 5 and 6, we have the expression for the Fourier transform of the exponential:

F λ e α(t+t 0 ) = e αt 0 ∞ ∑ n=0 c n α n 2πi n δ (n) (λ). (7)
The Fourier transform of a sine wave is a distribution given by:

F λ { sin(ω(t + t 0 )) } = iπ e -iωt 0 δ(λ+ω) -e iωt 0 δ(λ-ω) . (8)
We also know that δ is the neutral element of the convolution (see Eq. 17 in the Appendix). Moreover, we can convolve the distributions of the exponential and the sine since they both have compact support. From Eqs. 7 and 8, we finally have:

F λ { m(t + t 0 ) } = πe αt 0 ∞ ∑ n=0 c n α n i n+1 • e -iωt 0 δ (n) (λ+ω) -e iωt 0 δ (n) (λ-ω) . (9)
Convolution of Eq. 9 with a windowing function H leads to the desired short time Fourier transform of a mode. Using the properties of distributions (Eq. 17, 18 in the Appendix), and Eq. 9, we have:

s(λ) = 1 2 e αt 0 ∞ ∑ n=0 c n α n i n+1 e -iωt 0 F λ (H) (n) (λ+ω) -e iωt 0 F λ (H) (n) (λ-ω) . (10) F λ (H) (n) (λ+ω)
is the n-th derivative of the (complex) Fourier transform of the window H, taken at the value (λ+ω).

Note that Eq. 9 is still a distribution, since we did not constrain the mode to be computed only for positive times, and the mode itself is not square-integrable for negative times. However, this distribution has compact support which makes the convolution of Eq. 10 possible [START_REF] Hormander | The Analysis of Linear Partial Differential Operators I[END_REF]].

For computational efficiency, we truncate the infinite sum of Eq. 10, and approximate it by retaining only the first term. The final expression of our approximation to the mode STFT is thus:

s(λ) ≈ 1 2 e αt 0 c 0 i e -iωt 0 F λ (H)(λ+ω) -e iωt 0 F λ (H)(λ-ω) . (11)
Instead of c 0 = 1 (Eq. 5), we take c 0 to be the value of the exponential minimizing t 0 +∆t t 0 (e -αtc 0 ) 2 dt, where ∆t is the duration of a frame, resulting in a better piecewise constant approximation:

c 0 = e -αt 0 -e -α(t 0 +∆t) α ∆t . ( 12
)
This single term formula is computationally efficient since the Fourier transform of the window can be precomputed and tabulated, which is the only memory requirement. Moreover both complex exponentials are conjugate of each other meaning that we only need to compute one sine and one cosine.

Speedup and Numerical Validation

Consider a scene requiring the synthesis of M modes. Using a standard recursive time-domain solution [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF], and assuming 512-sample frames, the cost of the frame is M × 512 × C mt . The cost C mt of evaluating a mode in the time domain is 5 or 6 multiplies and adds, using the recursive formulation of Eq. 6 in [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]]. In our approach, assuming an average of B bins per mode, the cost of a frame is M × B × C ST FT plus the cost (once per frame) of the inverse FFT.

The cost C ST FT of evaluating Eq. 11, is about 25 operations. With a value of B = 3 (it is often lower in practice), we have a potential theoretical speedup factor of 30-40 times. If we take into account the fact that we have a 50% overlap due to windowing, this theoretical speedup factor drops to 15-20 times.

We have used our previous test scene to measure the speedup of our approach in practice, compared to [van den Doel and [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF].

When using B = 3 bins per mode, we found an average speedup of about 8, and with B = 5 bins per mode about 5. This reduction compared to the theoretical speedup is probably due to compiler and optimization issues of the two different algorithms.

Finally, we examine the error of our approximation for a single sound. We tested two different windows, a Hann window with 50% overlap and a Rectangular window with 10% inter-frame blending.

In Fig. 2 we show 3 frames, with the reference [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]] in red and our approximation in blue, taking B = 5 bins per mode.

Taken over a sequence including three impacts (a single pipe in the Magnet scene, see Sect. 7), for a total duration of about 1 sec., the average overall error for the Rectangular window is 15% with 5 bins, and 21% with 3 bins (it is 8% if we use all 512 bins). This error is mainly due to small ringing artifacts and possibly to our constant exponential approximation, which can be seen at frame ends (see Fig. 2). Using the Hann window, we have 35-36% error for both 512 and 5 bins, and 36% with 3 bins. This would indicate that the error is mainly due to the choice of window. As can be seen in the graph (Fig. 2 (right)) the error with the Hann window is almost entirely in the first frame, at the onset, or "attack", of the sound for the first 512 samples (i.e., 11 msec at 44.1kHz). The overall quality of the signal is thus preserved in most frames; in contrast, the ringing due to the rectangular window can result in audible artifacts. For this reason, and to be compatible with the pipeline of [START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF]], we chose to use the Hann window.

Limitations for the "Attacks" of Impact Sounds

Contrary to time domain approaches such as [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF][START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF], the use of the Hann window as well as the constant exponential approximation (CEA) degrades the onset or "attack" for high frequency modes. This attack is typically contained in the first few frames, depending on decay rate.

To study the individual effect of each of the CEA and the Hann window in our reconstruction process, we computed a time-domain solution using the CEA for the example of a falling box (Fig. 3(left))

and a time domain solution using a Hann window to reconstruct the signal (Fig. 3(right)). We plot the time-domain reference [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]] in red and the approximation in blue. As we can see, most of the error in the first 7msec is due to the Hann window whereas the CEA error remains lower. The effect of these approximations is the suppression of the "crispness" of the attacks of the impact sounds.

Use of the Hann window and the CEA as described previously has the benefit of allowing seamless integration between recorded and impact sounds, as described next in Sect. 5. In complex soundscapes containing many recorded sounds, this approximation may be acceptable. However, in other cases the crispness of the attacks of the modal sounds can be important. To better preserve the attacks of the impacts sounds, we treat them in a separate buffer and split the first 1024-sample frame into four subframes. Each subframe has a corresponding CEA and with a specialized window. In what follows, we assume that all contact sounds start at the beginning of a frame.

We design a windowing scheme satisfying four constraints: 1) avoid "ramping up" to avoid smoothing the attack, 2) end with a 512 sample square root of Hann window to blend with the buffer for all frames other than the attack, 3) achieve perfect reconstruction, i.e., all windows sum to one, 4) require a minimal number of bins overall, i.e., use Hann windows which minimize the number of bins required for reconstruction et al. 1999].

The first subframe is synthesized using a rectangular window for the first 128 samples (constraint 1) followed by half of a Hann window ("semi-Hann" from now on) for the next 128 samples and zeros in the remaining 768 samples; this is shown in blue in Fig. 4. The next two subframes use full 256-sample Hann windows, starting at samples 128 and 256 respectively (red and green in Fig. 4). The last subframe is composed of a semi-Hann window from samples 384 to 512 and a square root of a semi-Hann window for the last 512 samples for correct overlap with the non-attack frames, thus satisfying constraint 2 (black in Fig. 4). All windows sum to 1 (constraint 3), and Hann windows are used everywhere except for the first 128 samples (constraint 4). These four buffers are summed before performing the inverse FFT, replacing the original 1024 sample frame by the new combined frame. We use 15 bins in the first subframe. The increase in computational cost is negligible, since the additional operations are only performed in the first frame of each mode: for the example shown in the video the additional cost of the mixed window for attacks is 1.2%. However, the integration of this method with recorded sounds is somewhat more involved; we discuss this in Sect. 9.

A Full Perceptually Based Scalable Pipeline for Modal and Recorded Sounds

An inherent advantage of frequency domain processing is that it allows fine-grain scalable audio processing at the level of an FFT bin. In [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]], such an approach was proposed to perform equalization and mixing, reverberation processing and spatialization on prerecorded sounds. Signals are prioritized at runtime and a number of frequency bins are allocated to each source, thus respecting a predefined budget of operations. This importance sampling strategy is driven by the energy of each source at each frame and used to determine the cut-off point in the list of STFT coefficients.

Given our fast frequency-domain processing described above, we can also use such an approach. In addition to the fast STFT synthesis, we also require an estimation of energy, both of the entire impact sound, and of each individual mode. In [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]] sounds were pre-recorded, and the FFT bins were precomputed and presorted by decreasing energy.

For our scalable processing approach, we first fix an overall mixing budget, e.g., 10,000 bins to be mixed per audio frame. At each frame we compute the energy E s of each impact sound over the frame and allocate a budget of frequency bins per sound proportional to its energy. We compute the energy E m of each mode for the entire duration of the sound once, at the time of each impact, and we use this energy weighted by the mode's squared amplitude to proportionally allocate bins to modes within a sound. After experimenting with several values, we assign 5 bins to the 3 modes with highest energy, 3 bins for the next 6 and 1 bin for all remaining modes. We summarize these steps in the following pseudo-code.

1. PerImpactProcessing(ImpactSound S) // at impact notification 2. foreach mode of S 3. compute total energy E m 4.

Sort modes of S by decreasing E m 5.

Compute total energy of S for cutoff 6.

Schedule S for processing 1. ScalableAudioProcessing() // called at each audio frame 2. foreach sound S 3.

Compute E s 4.

Allocate FFT bin budget based on E s 5.

Modes m 1 , m 2 , m 3 get 5 bins 6.

Modes m 4m 9 get 3 bins 7.

1 bin to remaining modes until end of budget 8. endfor

Efficient Energy Estimation

To allocate the computation budget for each impact sound, we need to compute the energy, E s , of a modal sound s in a given frame, i.e., from time t to time t + ∆t:

E s = t + ∆t t s 2 (x)dx. ( 13
)
From Eq. 1 and 2, we express E s as:

E s = < s, s > = M ∑ i=0 M ∑ j=0 a i a j < m i , m j > . ( 14
)
For a given frame, the scalar product < m i , m j > has an analytic expression (see Eq. 22 in the additional material). Because this scalar product is symmetric, we only have to compute half of the required operations.

In our experiments, we observed that most of the energy of an impact sound is usually concentrated in a small number N of modes (typically 3). To identify the N modes with highest energy, we compute the total energy, E m , of each mode as:

E m = ∞ 0 sin(ωx)e -αx 2 dx = 1 4 ω 2 α(α 2 + ω 2 ) (15)
After computing the E m 's for each mode, we weight them by the square of the mode's amplitude. We then sort the modes by decreasing weighted energy. To evaluate Eq. 14 we only consider the N modes with highest energy. We re-use the result of this sort for budget allocation.

We also compute the total energy for a sound, which is used to determine its duration, typically when 99% of the energy has been played back. We use Eq. 14 and an expression for the total energy (rather than over a frame), given in the additional material (Eq. 21).

Numerical Validation

We use the test scene presented in Sect. 4.1 (Fig. 6) to perform numerical tests with appropriate values for N.

We evaluated the average error of the estimated energy, compared to a full computation. Using 3 modes, for all objects in this scene, the error is less than 9%; for 5 modes it falls to 4.9%.

A Complete Combined Audio Pipeline

In addition to frequency-domain scalable processing, we can also use the perceptual masking and sound-source clustering approaches developed in [START_REF] Tsingos | Perceptual audio rendering of complex virtual environments[END_REF][START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF]]. We can thus mix pre-recorded sounds, for which the STFT and energy have been precomputed, with our frequency domain representation for modal sounds and perform global budget allocation for all sounds. As a result, masking between sounds is taken into account, and we can cluster the surviving unmasked sources, thus optimizing the time for per-sound source operations such as spatialization. In previous work, the masking power of a sound also depends on a tonality indicator describing whether the signal is closer to a tone or a noise, noisier signals being stronger maskers. We computed tonality values using a spectral flatness measure [START_REF] Tsingos | Perceptual audio rendering of complex virtual environments[END_REF] for several modal sounds and obtained an average of 0.7. We use this constant value for all modal sounds in our masking pipeline.

Temporal Scheduling

One problem with systems simulating impact sounds is that a large number of events may happen in a very short time interval (debris from an explosion, a collapsing pile of objects, etc.), typically during a single frame of the physics simulation. As a result, all sounds will be triggered simultaneously resulting in a large peak in system load and possible artifacts in the audio ("cracks") or lower audio quality overall. Our idea is to spread out the peaks over time, exploiting results on audio-visual human perception.

As mentioned in Sect. 2, there has been extensive study of audiovisual asynchrony in neuroscience which indicates that the brain is able to compensate for the different delays between an auditory and a visual event in causal inference. To exploit this property, we introduce a scheduling step at the beginning of the treatment of each audio frame. In particular, we maintain a list of sound events proposed by the physics engine (which we call TempSoundsList) and a list of sounds currently processed by the audio engine (CurrSound-sList). At the beginning of each audio frame, we traverse Temp-SoundsList and add up to 20 new sounds to CurrSoundsList if one of the following is true:

• CurrSoundsList contains less than 50 sounds

• Sound s has been in TempSoundsList for more than T msec.

The values 20 and 50 were empirically chosen after experimentation. We use our perceptual tolerance to audio-visual asynchrony by manipulating the threshold value T for each sound. In particular we set T to a desired value, for example 200 msec corresponding to the results of [START_REF] Guski | Audiovisual phenomenal causality[END_REF]. We then further modulate T for sounds which are outside the visible frustum; T is increased progressively as the direction of the sound is further from the center of the field of view. For sounds completely behind the viewer the delay T is set to a maximum of 0.5 seconds. Temporal scheduling only occurs for impact sounds, and not for recorded sounds such as a gun firing which are time-critical and can have a remote effect.

Our approach reduces the number and density of computational peaks over time. Although some peaks still occur, they tend to be smaller and/or sharper, i.e., occur over a single frame (see Fig. 5). Since our interactive system uses buffered audio output, it can sustain such sparse peaks over a single frame, while it could not sustain such a computational overload over several consecutive frames.

Implementation and Results

Our system is built on the Ogre3D1 game engine, and uses the PhysX2 real-time physics engine simulator. Throughout this paper, we use our own (re-)implementations of [van den [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF]], [O' Brien et al. 2002] and [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF]. For [van den Doel and Pai 2003], we use T = 512 samples at 44.1KHz; the size of the impact filter with a force profile of cos(2πt/T ) is T = 0.37ms or 16 samples (Eq.17 of that paper).

For objects generating impact sounds, we precompute modes using the method of O 'Brien et al. [2002]. Sound radiation amplitudes of each mode are estimated with a far-field radiation model (Eq. 15, [START_REF] James | Precomputed acoustic transfer: Output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF]).

Audio processing was performed using our in-house audio engine, with appropriate links between the graphics and audio. The audio engine is described in detail in [START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF]]. In total we run three separate threads: one for each of the physics engine, graphics and audio. All timings are reported on a dual-processor, dual-core Xeon running at 2.3Ghz.

Interactive Sessions Using the Pipeline

We have constructed four main scenes for demonstration purposes, which we refer to as "Oriental", "Magnet", "Truck" and "Boxes"; we show snapshots of each in Fig. 1 (Truck and Oriental) and 6. The goal was to construct environments which are similar to typical simulators or games settings, and which include a large number of impact sounds, as well as several prerecorded sounds. All sounds were processed in the Fourier-domain at 44.1kHz using 1024-tap FFTs and 50% overlap-add reconstruction with Hann windowing. The "Oriental" and "Box" scenes contain modal sounds only and thus use the attack preserving approach (Sect. 4.3). Hence, our audio thread runs at 83Hz; we then output reconstructed audio frames of 512 samples. The physics thread updates object motion at 140Hz, and the video thread runs at between 30-60Hz depending on the scene and the rendering quality (shadows, etc.).

The Magnet scene contains prerecorded industrial machinery and closing door sounds, while the Truck scene contains traffic and helicopter sounds. Basic scene statistics are given in

Quality and Performance

We performed two comparisons for our fast modal synthesis: the first was with the "standard" time-domain (TD) method of [van den Doel and Pai 2003] using recursive evaluation, and the second with the mode-culling time-domain approach of [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF], which is the fastest TD method to date. We used the "Oriental" scene for this comparison, containing only modal sounds.

Comparison to "standard" TD synthesis

Comparison to mode-culling TD synthesis

To compare to modeculling, we apply the mode truncation and quality scaling stages of [START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF] at each audio frame. We then perform fast frequency domain synthesis for the modes which have not been culled. For the same mode budget our frequency processing allows a speedup of 4-8 times; the difference in speedup with the "standard" TD synthesis is due to implementation issues. The quality of the two approaches is slightly different for the same mode budget, but in both cases subjectively gives satisfactory results.

Full Combined Pipeline

The above comparisons are restricted to modal sounds only. We also present results for other scenes, augmented with recorded sounds. These are rendered using our full pipeline, at low overall budgets but with satisfactory quality.

We present statistics for our approach in Table 2, using a budget of 8000 bins. First we indicate the cost (in milliseconds) of each component of our new combined pipeline: mixing, energy computation, masking and clustering, as well as the total cost. As we can see there is a non-negligible overhead of the pipeline stages; however the benefit of being able to globally allocate budget across modal and recorded sounds, and of course all the perceptually based accelerations, justifies this cost.

The number of sounds at each frame over the course of the interactive sequences shown in the video varied between 195 and 970. If no masking or culling were applied there would be between 30,000 to 100,000 modes to be played per audio frame on average in these sequences. We use 15,000 to 20,000 frequency bins in all interactive sessions. The percentage of prerecorded sounds masked was around 50% on average and that of impact sounds was around 30%.

Pilot Perceptual Evaluation

Despite previous experimental studies for perceptually based audio rendering for pre-recorded sounds [START_REF] Moeck | Progressive perceptual audio rendering of complex scenes[END_REF][START_REF] Tsingos | Perceptual audio rendering of complex virtual environments[END_REF], and the original neuroscience experiments for asynchrony [START_REF] Guski | Audiovisual phenomenal causality[END_REF], we consider it imperative to conduct our own pilot study, since our context is very different. We have two conditions in our experiment: the goal of the first condition is to evaluate the overall audio quality of our approximations while that of the second is to evaluate our temporal scheduling.

Experiment Setup and Procedure

In our experiment we used the Magnet and Truck (Sect. 7) environments, but with fewer objects, to avoid making the task too hard for the subjects. We used two 6 second pre-recorded paths for each of the two scenes. To ensure that all the stimuli represent the exact same sequence of events and to allow the presentation of a reference in real-time, we synchronize all threads and store the output audio and video frames of our application to disk. Evidently any other delay in the system has to be taken into account. We choose the parameters of our simulation to be such that we do not perceive "cracks" in the audio when running interactively with the same budget settings. Video sequences are then played back during the study.

For the reference sequences, all contact sounds are computed in the time-domain and no perceptual processing is applied when mixing with the recorded sounds in the frequency domain. We use non-individualized binaural rendering using Head Related Transfer Functions (HRTFs) chosen from the Listen database3 .

The interface is a MUSHRA-like [ITU 2001[ITU -2003] ] slider panel (see the accompanying video), in which the user can choose between a reference and five different approximations (A, B, C, D, E), each with a different budget of frequency bins. The subject uses a slider to rate the quality of each stimulus. One of the five stimuli is a hidden reference. The radio button above each slider allows the subject to (re)start the corresponding sequence. For the synchronization condition, we choose one budget which has a good rating in the quality condition (typically C 3 in Table 3), and delay the audio relative to graphics by a variable threshold T . The budgets and thresholds used are shown in Table 3. A tick box is added under each sound for synchrony judgment. The experiment interface can be seen in the accompanying video.

The subject listens to the audio with headphones and is instructed to attend to both visuals and audio. There are 8 panels, corresponding to all the conditions; stimuli are presented in random order. Detailed instructions are given on screen to the user, who is asked to rate the quality and tick if asynchrony between audio and visuals is detected. Rating of each panel is limited to 3 minutes, at which point rating is disabled.

Analysis of the Experiments

We ran the experiment with 21 subjects who were members of our research institutes, and were all naive about the goal of our experiments. We show the average quality ratings and the percent perceived asynchrony averages for the experiment in Table 4.

As we can see, for the Magnet scene, the budget of 4,000 bins was sufficient to give quality ratings of 83-85% very close to the hidden reference rated at 84-89%. An analysis of variance (ANOVA) [START_REF] Howell | Statistical Methods for Psychology[END_REF]] with repeated measures on quality ratings shows a main effect of budget on perceived quality (F(4,20)=84.8, p<0.00001). For the Truck scene the quality rating for 8,000 bins was lower. This is possibly due to the fact that the recorded sounds require a significant part of the frequency bin budget, and as a result lower the overall perceived quality.

In terms of asynchrony, the results have high variance. However, it is clear that audiovisual asynchrony was perceived less than 25% of the time, for delays under 200msec.

Overall, we consider these results to be a satisfactory indication that our approximations work well both in terms of progressive processing and for our temporal scheduling algorithm. In particular, there is a strong indication that increasing the budget does result in perceptually improved sounds, and that only a small percentage of users perceive asynchrony with temporal scheduling with delays less than 200ms.

Discussion and Conclusions

We have presented a new frequency-domain approach to modal sound rendering, which exploits sparseness of the Fourier Transform of modal sounds, leading to an 4-8 speedup compared to time-domain approaches [START_REF] Van Den Doel | Modal synthesis for vibrating objects[END_REF][START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF], with slight quality degradation. Furthermore, our approach allows us to introduce a combined perceptual audio pipeline, treating both prerecorded and on-the-fly impact sounds, and exploiting scalable processing, auditory masking, and clustering of sound sources. We used crossmodal results on perception of audiovisual synchronization to smooth out computational peaks which are frequently caused by impact sounds, and we performed a pilot perceptual study to evaluate our combined pipeline.

Use of the Hann window allows direct integration of modal and recorded sounds (see Sect. 5), but leads to lower quality attacks of impact sounds. We have developed a solution to this problem, splitting the treatment of the first frame of each attack into four subframes with appropriate windows. This solution can be easily used in scenes exclusively containing modal sounds. For the pipeline combining recorded and modal sounds, and in particular for clustering, we would need a separate buffer for attacks in each cluster thus performing twice as much post-processing (HRTF processing etc.). The rest of the pipeline would remain unchanged.

We developed an initial solution for rolling sounds in our pipeline, using a noise-like sound profile, similar in spirit to [START_REF] Van Den Doel | FoleyAutomatic: physically-based sound effects for interactive simulation and animation[END_REF]]. To simulate the continuous rolling excitation, we precompute a noise profile in the Fourier domain, and perform dynamic low-pass filtering based on velocity. Convolution with the excitation is a simple product in the Fourier domain, making our approach efficient. The accompanying video contains a first example. Nonetheless, a general solution to continuous excitation in our framework requires mixing delayed copies of past frames, incurring additional costs. We expect masking and progressive processing to limit this overhead, similar to the reverberation in [START_REF] Tsingos | Scalable perceptual mixing and filtering of audio signals using an augmented spectral representation[END_REF]].

Another limitation is the overhead of our pipeline which is not negligible. For practical usage of real-time audio rendering, such as game engines, we believe that the benefits outweigh this drawback.

In addition to the perceptually based accelerations, we believe that the ability to treat recorded and synthesized sounds in a unified manner is very important for practical applications such as games.

Currently, when reducing the budget very aggressively the energy computation can become a dominant cost. It is possible to precompute a restricted version of the energy, if we assume that forces are always applied in the normal direction. This is the case for recording-based systems (e.g., [START_REF] Van Den Doel | The sounds of physical shapes[END_REF]). However, this reduces the flexibility of the system.

Figure 1 :

 1 Figure 1: Frequency-domain fast mode summation: Top: Frames of some of the test scenes our method can render in real time. Bottom: (Left) Time-domain modal synthesis requires summing all modes at every sample. (Right) Our frequency-domain modal synthesis exploits the inherent sparsity of the modes' discrete Fourier transforms to obtain lower costs per frame.

Figure 2 :

 2 Figure 2: Comparison of Reference with Hann window and with Rectangular window reconstruction using a 1024-tap FFT.

Figure 3 :

 3 Figure 3: Comparison of Constant Exponential Approximation in the time domain (TD) and Hann window reconstruction in the TD, with the reference for the sharp sound of a falling box.

Figure 4 :

 4 Figure4: Four sub-windows to better preserve the sound attack.

Figure 5 :

 5 Figure5: Effect of temporal scheduling; computational peaks are delayed, the slope of increase in computation time is smoothed out and the duration of peaks is reduced. Data for the Magnet sequence, starting just before the group of objects hits the floor.

  In terms of quality, we tested examples of our frequency-domain synthesis with 3 and 5 bins per mode, together with a time-domain reference, shown in the accompanying video. The quality for 5 bins is very close to the reference. The observed speedup was 5-8 times, compared to [van den Doel and Pai 2003].

Figure 6 :

 6 Figure 6: From left to right, snapshots of the large Magnet scene, the Boxes scene, the test scene used for numerical validation and our prototype rolling demo. Please see and listen to the accompanying video showing these scenes.

Table 1

 1 , for the demo versions of the scenes shown in the video.

	Scene	O	T	P	M i	M/o
	Oriental 173 730K	0	665	214
	Boxes	200 200K	0	678	376
	Magnet	110 300K 16 971	164
	Truck	214 600K 15 268	221

Table 1 :

 1 Basic statistics for example scenes. O: number of objects simulated by the physics engine producing contact sounds, T : total number of triangles in the scene, P: number of pre-recorded sounds in the scene and M i : maximum number of impact sounds played per frame. M/o: average number of modes/object.

Table 2 :

 2 Cost in milliseconds of each stage of our full pipeline.

	Scene	Total Mixing Energy Masking Clustering
	Magnet	3.2	1.3	0.6	1.0	0.3
	Truck	2.7	1.5	0.9	0.3	0.1

Table 3 :

 3 Budget and delay values used for the perceptual experiments. C i and T i are the budget and delay conditions used.

			Percent Perceived Quality		% Perceived Asynchrony
	Scene	C 1	C 2	C 3	C 4	Re f	T 1	T 2	T 3	T 4
	Magnet1	51.1	64.7	78.0	83.1	84.9	0	24	48	71
	Magnet2	48.8	70.1	76.5	85.2	88.9	10	5	0	10
	Truck1	26.3	41.8	54.2	66.0	87.3	24	24	14	38
	Truck2	24.8	28.6	42.7	66.0	89.0	14	14	38	29

Table 4 :

 4 Results of our experiments: average quality and percent perceived asynchrony for the 2 scenes and 2 paths.

http://www.ogre3d.org

http://www.ageia.com

http://recherche.ircam.fr/equipes/salles/listen/

Acknowledgments This research was partly funded by the EU FET project CROSSMOD (014891-2 http://www.crossmod.org). We thank Autodesk for the donation of Maya and Alexandre Olivier-Mangon and Fernanda Andrade Cabral for modeling the examples. We thank the anonymous reviewers whose suggestions made this a much better paper. The last author acknowledges a fellowship from the Alfred P. Sloan Foundation.

Appendix

A.1 Some Elements of Distribution Theory

Applying a distribution T to a smooth test function with local support f , implies the following operation:

A commonly used distribution is the Dirac distribution (note that this is not the Kronecker delta) which has value 0 everywhere, except at 0. < δ k , f >= f (k) is commonly used in signal processing (Dirac combs). We use the following properties of distributions:

where f (n) denotes the n th derivative of f .

F ( f (t)g(t)) = 1 2π


Mihir Jain

Hervé Jégou

Patrick Gros

Asymmetric Hamming Embedding Taking the best of our bits for large scale image search

This paper proposes an asymmetric Hamming Embedding scheme for large scale image search based on local descriptors. The comparison of two descriptors relies on an vectorto-binary code comparison, which limits the quantization error associated with the query compared with the original Hamming Embedding method. The approach is used in combination with an inverted file structure that offers high efficiency, comparable to that of a regular bag-of-features retrieval systems. The comparison is performed on two popular datasets. Our method consistently improves the search quality over the symmetric version. The trade-off between memory usage and precision is evaluated, showing that the method is especially useful for short binary signatures.

INTRODUCTION

Large scale image search is still a very active domain. The task consists in finding in a large set of images the ones that best resemble the query image. Typical applications include finding searching images on web [START_REF] Wu | Bundling features for large scale partial-duplicate web image search[END_REF], location [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF] or particular object [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] recognition, or copy detection [START_REF] Law-To | Video copy detection: a comparative study[END_REF]. Earlier approaches were based on global descriptors such as color histograms or GIST [START_REF] Oliva | Modeling the shape of the scene: a holistic representation of the spatial envelope[END_REF]. These are sufficient in some contexts [START_REF] Douze | Evaluation of GIST descriptors for web-scale image search[END_REF], such as copy detection, where most of the illegal copies are very similar to the original image. However, global description suffer from well-known limitations, in particular they are not invariant to significant geometrical transformations such as cropping. Here we focus on the bag-of-words (BOW) framework [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] and its extension [START_REF] Jégou | On the burstiness of visual elements[END_REF], where local descriptors are extracted from each image [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and used to compare images.

The BOW representation of images was proved be very discriminant and efficient for image search on millions of images [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF][START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF]. Different strategies have been proposed to improve it. For instance, [START_REF] Nistér | Scalable recognition with a vocabulary tree[END_REF] improves the efficiency in two ways. Firstly, the assignment of local descriptors to the so-called visual words is much faster thanks to the use of a hierarchical vocabulary. Secondly, by considering large vocabularies (up to 1 million visual words), the size of the inverted lists used for indexing is significantly reduced. Accuracy is improved by a re-ranking stage performing spatial verification [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], Pre-print (author) version. To be published in ACM Multimedia 2011. and by query expansion [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF], which exploits the interaction between the relevant database images.

Another way to improve accuracy consists in incorporating additional information on descriptors directly in the inverted file. This idea was first explored in [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF], where a richer descriptor representation is obtained by Hamming Embedding (HE) and weak geometrical consistency [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF]. HE, in particular, was shown successful in different contexts [START_REF] Wu | Bundling features for large scale partial-duplicate web image search[END_REF], and improved in [START_REF] Jégou | On the burstiness of visual elements[END_REF][START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF]. However, this technique has a drawback: each local descriptor is represented by relatively large signatures, typically ranging from 32 [START_REF] Wu | Bundling features for large scale partial-duplicate web image search[END_REF] to 64 bits [START_REF] Jégou | On the burstiness of visual elements[END_REF].

In this paper, we propose to improve HE in order to better exploit the information conveyed by the binary signature. This is done by exploiting the observation, first made in [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF], that the query should not be approximated. We therefore adapt the voting method to better exploit the precise query location instead of the binarized query vector. This requires, in particular, two regularization steps used to adjust the dynamic of the local query distribution. This leads to an improvement over the reference symmetric HE scheme. As a complementary contribution, we evaluate how our approach trades accuracy against memory with smaller number of bits. To our knowledge, such a comparison has never been published.

The paper is organized as follows. The datasets representing the application cases and the evaluation protocol are introduced in Section 2. Section 3 briefly describes the most related works: BOW and HE. Our asymmetric method is introduced in Section 4. Finally, experiments in Section 5 compare the performance of our asymmetrical method with the original HE, and provides a comparison with the state of the art on image search. Its shows a significant improvement: we obtain a mean average precision of 70.4% on the Oxford5K Building dataset before spatial verification, i.e., +4% compared with the best concurrent method.

EVALUATION DATASETS

This section introduces the datasets used in our experiments, as well as the measures of accuracy used to evaluate the different methods. These datasets reflect two application use-cases for which our method is relevant, namely place and object recognition. They are widely used to evaluate image search systems. Oxford5K and Paris. These two datasets of famous building in Oxford and Paris contain 5,062 and 6,412 images, respectively. We use Paris as an independent learning set to estimate the parameters used by our method. The quality is measured on Oxford5K by mean average precision (mAP), as defined in [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF]: for each query image we obtain a precision/recall curve, and compute its average precision (the area under the curve). The mAP is then the mean for a set of queries. INRIA Holidays. This dataset contains 1491 images of personal holiday photos, partitioned into 500 groups, each of which it is a modified inverted file. Each inverted list is associated with a visual word. Each local descriptor is stored in a cell containing both the image identifier and a binary signature (from 4 to 64 bits in our paper). This indexing structure is also used in our AHE method: only the score similarity is modified.

represents a distinct scene, location or object. The first image of each group is the query image and the correct retrieval results are the other images of the group. Again, the search quality is measured by the mAP, see [START_REF] Philbin | Object retrieval with large vocabularies and fast spatial matching[END_REF][START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] for details. A set of images from Flickr is used for learning the vocabulary, as done in [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF].

Flickr1M. In order to evaluate the behavior of our method on a large scale, we have used a set of up to one million images. More precisely, we have used the descriptors shared online 1 by Jegou et al., which were downloaded from Flickr and described by the same local descriptor generation procedure as the one used for Holidays. This dataset is therefore merged with Holidays and the mAP is measured using the Holidays ground-truth, as in [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF].

The recall@R measure is used for this large scale experiment. It measures, at a particular rank R, the ratio of relevant images ranked in top R positions. [START_REF] Douze | Evaluation of GIST descriptors for web-scale image search[END_REF] states that it is a good measure to evaluate the filtering capability of an image search system, in particular if the large scale image search is followed by a precise spatial geometrical stage, as classically done in the literature.

RELATED WORK

Bag-of-features representation

The BOW framework [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] is based on local invariant descriptors [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF][START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] extracted from covariant regions of interest [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF]. It matches small parts of images and can cope with many transformations, such as scaling, local changes in illumination and cropping.

The feature extraction is performed in two steps: detecting regions of interest with the Hessian-Affine detector [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF], and computing SIFT descriptors for these regions [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF]. We have used the features provided by the authors for all the datasets.

The fingerprint of an image is obtained by quantizing the local descriptors using a nearest-neighbor quantizer, produced the socalled visual words (VW). The image is represented by the histogram of VW occurrences normalized with the L2 norm. A tf-idf weighting scheme is applied [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] to the k components of the resulting vector. The similarity measure between two BOW vectors is, most often, cosine similarity. The visual vocabulary of the quantizer is produced using k-means. It contains a large number k of visual words. In this paper, k = 20, 000 for the sake of consistency with [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] and [START_REF] Jégou | On the burstiness of visual elements[END_REF]. Therefore, the fingerprint histograms are sparse, making queries in the inverted file efficient.

Hamming embedding

The Hamming Embedding method of [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] is a state of the art method extension of BOW, where a better representation of the im-1 http://lear.inrialpes.fr/people/jegou/data.php ages is obtained by adding a short signature that refines the representation of each local descriptor. In this approach, a descriptor x is represented by a tuple (q(x), b(x)), where q(x) is the visual word and b(.) is a binary signature of length m computed from the descriptors to refine the information provided by q(x). Two descriptors are assumed to match if and ht is a fixed threshold. The image score is obtained as the sum [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] or weighted sum [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF] of the distances of the matches satisfying (1), then normalized as BOW.

q(x) = q(y) h (b(x), b(y)) = i=1..m |bi(x) -bi(y)| ≤ ht , (1)
As BOW, the method uses an inverted file structure, which is modified to incorporate the binary signature, as illustrated by Figure 1. The matches associated with the query local descriptors are retrieved from the structure and used as follows.

• Find the nearest centroid of the query descriptor x, producing quantized indexes q(x), i.e., the visual word (VW). The entries of the inverted list associated with q(x) are visited.

• Compute the vector Q × x = b * = [b * 1 (x), . . . , b * m (x)
] associated with descriptor x, where Q is a rotation matrix for which only the first m rows are kept so that Q × x ∈ R m .

• The binary signature is obtained by comparing each component b * i , i = 1..m with a threshold τ q(x),i . This amounts to selecting bi = 1 if b * i -τ q(x),i > 0, else bi = 0. The thresholds τc,i are the median values of b * i measured on an independent learning set for all VWs c and all bit components i.

• Only the database descriptors satisfying Equation 1 make a vote for the corresponding image, i.e., they vote only if their Hamming distance is below a pre-defined threshold ht. The vote's score is 1 in [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF]. Scoring with a function of the distance improves the results [START_REF] Jégou | On the burstiness of visual elements[END_REF]. We therefore adopt this choice.

• All images scores are finally normalized.

Additionally, we consider in this paper two techniques [START_REF] Jégou | On the burstiness of visual elements[END_REF] that improve the results. First, multiple assignment (MA) reduces the number of matches that are missed due to incorrect quantization indexes. Second, the so-called burstiness (denoted by "burst") handling method regularizes the score associated with each match, to compensate the bursty statistics of regular patterns in images.

ASYMMETRIC HAMMING EMBEDDING

This section introduces our approach. It is inspired by the works of Dong [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF] and Gorda [START_REF] Gordo | Asymmetric distances for binary embeddings[END_REF], where the use of asymmetric distances was investigated in the context of Locality Sensitive hashing. This method has to be significantly adapted in our context. Using the distance to hyperplanes may suffice for pure nearest neighbor search, where the objective is the find the Euclidean k-nearest neighbor of a given query [START_REF] Dong | Asymmetric distance estimation with sketches for similarity search in high-dimensional spaces[END_REF]. However, in our case, this is not sufficient, because computed distances are used as match quality measurements. Our goal is therefore to provide a soft weighting strategy that better exploits the confidence measures of all matches to produce the aggregated image scores.

Intra-cell distance regularization. We first adapt the local distances so that they become more comparable for different visual words. This is done, in our AHE scheme, by computing the standard deviation σc of the distance b * i -τc,i to the separating hyperplanes for each of the k visual words c. This estimation is carried out using a large set of vectors from an independent learning set. We used 50M Flickr descriptors for Holidays and all the descriptors from Paris for Oxford5K. The standard deviation is either computed component-wise (one per bit dimension) of for the whole cell. In our case we chose the simple choice of estimating a single parameter per cell used for all bits (isotropic assumption). As observed in Figure 2, the standard deviations significantly vary from one cell to another. It is then worth obtaining more comparable values when using distances as quality measurements. Distance to hyperplanes. In the symmetric version, the query projected by Q is binarized and compared with the database descriptors. We instead compute the distance between the projected query (b * (x) = Q × x) and the database binary vectors that lie in the same cell (associated with q(x)). The "distance" between the i th component of b * (x) and the binary vector b(y) is given by

d i a (b * i (x), bi(y)) = |b * i (x) -τ q(x),i | × |bi(x) -bi(y)| . ( 2
)
This quantity is zero when x is on the same side of the hyperplane associated with the i th component. The distances are added for all the m components to get an asymmetric "distance" between a query descriptor x and a database descriptor y, defined as

ha (b * (x), b(y)) = 1 σ q(x) × i=1..m d i a (b * i (x), bi(y)) . (3)
The descriptors are assumed to match if ha (b * (x), b(y)) ≤ ht, as for the symmetric version. For a given query x, the values |b * i (x) -τ q(x),i | are precomputed before being compared to the database vectors. The similarity is penalized according to the distance from the hyperplane in the embedded Hamming space, providing improved measurements, as illustrated by Figure 3. In the symmetric case, it does not matter how far b * i lies from τ q(x),i . In contrast, the distance is a continuous function in our method. Score weighting. Similar to what is done in [START_REF] Jégou | On the burstiness of visual elements[END_REF] for the symmetric case, the distance obtained by Equation 3 is used to weight Oxford5K Holidays BOW [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] 40.3 -BOW+soft MA [START_REF] Sivic | Video Google: A text retrieval approach to object matching in videos[END_REF] 49.3 -HE+MA5 [START_REF] Jégou | On the burstiness of visual elements[END_REF] 61.5 77.5 HE+burst [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] 64.5 78.0 HE+burst+MA5 [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] 67. the voting score. In [START_REF] Jégou | On the burstiness of visual elements[END_REF] weights are obtained as a gaussian function of Hamming distance. Here the weights are simply the difference ht -ha (b * (x), b(y)) between the threshold and the normalized "distance". We also apply the burstiness regularization method of [START_REF] Jégou | On the burstiness of visual elements[END_REF]. As we will show in Section 5, its impact is very important in our case because the aforementioned variance regularization does not sufficiently balance the amount of score received by the different query vectors, leading individual descriptors from the query image to have a very different impact in the final score. The burstiness regularization effectively addresses this issue.

EXPERIMENTS

Search quality: HE vs AHE. Figure 4 evaluates our AHE method introduced in Section 4 against the original HE one, for varying numbers of bits. For both HE and AHE, we report the results obtained with the best threshold. As shown by Figure 5, the performance is stable around this best value.

Using the asymmetric version significantly improves the results, especially for short signatures. As stated in Section 4, the burstiness regularization of [START_REF] Jégou | On the burstiness of visual elements[END_REF] is important in our case: without it AHE only achieves a slight improvement for short signatures.

Observe the important trade-off between the search accuracy and the signature length: using more bits clearly helps. However it is important to keep this signature short so that database images remain indexed in memory. Multiple Assignment helps in both cases, at the cost of increased query time. Unless specified, we used MA to the 5 nearest visual words, denoted by MA5.

Comparison with the state of the art. Table 1 compares our results with, to our knowledge, the best ones reported in the literature. Our approach clearly outperforms the state of the art on both Holidays and Oxford5K. Interestingly, for symmetric HE, our results (in italics) on Oxford5K are better than those reported in [START_REF] Jégou | On the burstiness of visual elements[END_REF] with a geometry check. This is because on these datasets, the features used in [START_REF] Jégou | On the burstiness of visual elements[END_REF] are not as good as those we used here. We only include in our comparison the results reported with learning done on an independent dataset itself. Some papers shows that learning on the test set itself improves the results, as to be expected. But as stated in [START_REF] Jégou | Improving bag-of-features for large scale image search[END_REF] such results do not properly reflect the expected accuracy when using the system on a large scale.

Large scale experiments. As shown in Figure 6, the filtering capability of AHE is better than HE: the recall@R measure is almost as good for AHE with 16 bits as HE with 32 bits. Equivalently, the performance is much better for a given memory usage.

The complexity of the method is increased compared to the original symmetric method. In both HE and AHE, the vector has to be projected. The main difference appears in the similarity computation, which is a simple XOR operation following by a bit count in HE, while we need to add pre-computed floating point values to get ha in Equation 3. As a result, on 1 million images the search time is roughly 1.7 times slower in the asymmetric case, on average, compared to [START_REF] Jégou | On the burstiness of visual elements[END_REF]: on average searching in one million images (using 64 bits) with AHE it takes 2.9s on one processor core, against 1.7s for HE.

Figure 1 :

 1 Figure1: Overview of the HE indexing structure: it is a modified inverted file. Each inverted list is associated with a visual word. Each local descriptor is stored in a cell containing both the image identifier and a binary signature (from 4 to 64 bits in our paper). This indexing structure is also used in our AHE method: only the score similarity is modified.

  where h(b, b ) is the Hamming distance between binary vectors b = [b1, .., bm] and b = [b 1 , .., b m ],

Figure 2 :

 2 Figure 2: Empirical probability distribution function of σc measured for all visual words. The large variation of density across cells shows the need for a per-cell variance regularization.

Figure 3 :

 3 Figure 3: Illustration of HE and AHE for binary signatures of length 2. In the symmetric case, only three distances are possible (0, 1 or 2) between query and database descriptor y. AHE gives a continuous distance (reflected by the intensity of blue).

Figure 5 :

 5 Figure 5: Impact of the threshold on accuracy (m = 32 bits). Note that the ranges for ht differ for HE (Hamming distance) and AHE (derived from normalized distance to hyperplanes).

Figure 4 :

 4 Figure 4: HE vs AHE: Trade-off between memory usage (per descriptor) and search quality.

Figure 6 :

 6 Figure 6: Search quality on a large scale (1 million images): Holidays merged with Flickr1M.

Table 1 :

 1 State of the art. For[START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF], we report in italics the results obtained with the best descriptors ([START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] reports inferior results with different descriptors and with geometrical information only).

		4	79.6
	AHE+burst	66.0	79.4
	AHE+burst+MA5	69.8	81.9
	AHE+burst+MA10	70.4	81.7

Acknowledgments:

This work was realized as part of the Quaero Project, funded by OSEO, French State agency for innovation.

CONCLUSION

This paper shows that a vector-to-binary code comparison significantly improves the state-of-the-art Hamming Embedding technique by reducing the approximation made on the query. This is done by exploiting the vector-to-hyperplane distances. The improvement is obtained at no additional cost in terms of memory. As a result, we improve the best results ever reported on two popular image search benchmarks before geometrical verification.


Héctor De Arazoza

Rachid Lounes
email: rachid.lounes@parisdescartes.fr.

Andres Sánchez

Jorge Barrios

Ying-Hen Hsieh
email: hsieh@mail.cmu.edu.tw.

MODELING DETECTION OF HIV IN CUBA

A nonlinear compartmental model is developed for the HIV detection system in Cuba with different types of detections, some random and others non-random. We analyze the dynamics of this system, compute the reproduction numbers, and use the data from the Cuban HIV/AIDS epidemic between 1986-2008 to fit the model. We obtain estimates for the detection-related parameters during two separate time periods to reflect the timeline of the implementation of various types of searches. The reproduction numbers for each time period are also computed from the sets of values of the parameters. We found that random screening is most important as a mean of surveillance. Moreover, local asymptotic stability for the Disease Free Equilibrium can be achieved if (i) random screening is sufficiently effective and (ii) infection by detected HIV-positive individuals is minimal. Our results highlight the importance of education for the known infectious for the purpose of preventing further infection. Fitting the 1986-2008 HIV data to obtain the model parameter estimates indicates that the HIV epidemic in Cuba is currently approaching an endemic equilibrium.

Introduction

The Cuban HIV/AIDS program was established in 1983 and the detection of the first HIV-positive person in Cuba took place in December 1985. The first AIDS case was diagnosed in Cuba in April 1986, this signaled the official start of the AIDS epidemic in the country. The Cuban HIV/AIDS epidemic has the lowest prevalence rate in the Caribbean region [START_REF] Inciardi | HIV/AIDS in the Caribbean Basin[END_REF]. The UNAIDS Epidemiological Fact Sheet on HIV and AIDS for Cuba reports an HIV prevalence of less than 0.1 % for adults [11]. The Cuban HIV/AIDS program includes a detection system that allows for detection of HIV-positive cases from several sources. Some of these sources were started at the beginning of the program, while others were introduced later and some have been discontinued in time. Since 1993, this detection system is composed, among others, of 6 major sources. These were screenings of all blood donors, persons that were treated for other sexually transmitted infections, persons admitted to a hospital with a suspicion of HIV infection or subject to specific procedures like dialysis, persons that volunteered to be tested, persons having received a recommendation for HIV testing from his/her general practitioner (family doctor) and through sexual partner tracing. Other minor sources include testing of all pregnant women and prison inmates [START_REF] De Arazoza | The HIV/AIDS epidemic in Cuba: Description and tentative explanation of its low HIV prevalence[END_REF]. From 1986 up to 2008 more than 34 million tests were performed, in recent years since 2002, the total number of tests performed has stabilized to 1.6-1.7 million test every year [START_REF] Minsap | CUBA:PLAN ESTRATEGICO NACIONAL ITS/VIH/SIDA 2007-2011[END_REF]. From 1986, in keeping with the "partner notification program", a person tested to be HIV-positive is invited to give names and contact details of his/her sexual partners during the past two years. These partners are then traced and a recommendation for voluntary HIV testing is made. The detection system has changed over the years. For example, during the period 1986-1999, the partner notification and contact tracing detected 30.7% of the new HIV infections while the Family Doctors Program detected only 8.9% of new infections reported. From 2000 and up to September 2008, the family doctors were responsible for the detection of 31.7% of the new cases while contact tracing detected 20.4% of the new cases. From this information we can deduce that the detection system is not static in terms of the contribution to the detection of new HIV cases. Before 1999 the most important contribution was from the contact tracing and partner notification program, while after 1999 the family doctors play a dominant role in the detection system. We will focus on modeling this change in the detection system by dividing the time period into two periods and introducing model parameters that will differentiate the detection by the family doctors from contact tracing and other (random) searches.

The model

As noted earlier, the primary objective of the Cuban Program to control the HIV/AIDS epidemic is the active search of persons infected with HIV long before they show any signs of AIDS. Our focus is not to model how new infections by HIV are generated, but how the HIV-infected persons are detected. We will consider the following variables and parameters. Model Variables: X(t) the number of HIV infected persons that do not know they are infected (at time t), Y (t) the number of HIV infected persons that know they are infected, Z(t) the number of persons with AIDS. Model Parameters:

1. λ the rate of recruitment of new HIV infected persons infected by X, 2. λ ′ the rate of recruitment of new HIV infected persons infected by Y , 3. k 1 , k 2 and k 3 rates at which the unknown HIV infected persons are detected by the system through different methods of detection. 4. β the rate at which the undetected HIV-positive persons develop AIDS, reciprocal of the mean incubation 5. β ′ the rate at which the detected HIV-positive persons develop AIDS, the reciprocal of the mean time it takes to go from Y to Z, 6. µ the mortality rate of the sexually active population, 7. µ ′ the mortality rate of the population with AIDS.

The dynamics is described by the following system:

dX dt = λX + λ ′ Y -k 1 X -k 2 XY -k 3 X 2 -βX -µX, dY dt = k 1 X + k 2 XY + k 3 X 2 -β ′ Y -µY, (1)
dZ dt = βX + β ′ Y -µ ′ Z.
We denote γ = β + µ, σ = λk 1γ, γ ′ = β ′ + µ and we consider the system only in the region D = {X ≥ 0, Y ≥ 0, Z ≥ 0}. It is clear that D is positively invariant under the flow induced by [START_REF] Anderson | Potential of community-wide chemotherapy or immunotherapy to control the spread of HIV-1[END_REF]. We give the following remarks regarding model (1):

1. There are three ways for individuals to go from "unknown HIV infected" (X) to "known HIV infected" (Y ). One is through the nonlinear term k 2 XY for contact tracing, where the individual is found through his contacts with persons that are known to live with HIV. The term k 3 X 2 models the detection through family doctors. The third way they can be detected is through the term k 1 X which models all the other "random" ways of searching for seropositives. It is important to note that 1/k 1 can be viewed as the mean time from infection to detection for the persons found through a random screening. We can look at all the ways that an "unknown" can go from X to Y as a term of the form "F (X, Y ) X", where F (X, Y ) is a recruitment function from the class X into the class Y . For our study we take F (X, Y ) as a polynomial of degree 1:

k 1 + k 2 Y + k 3 X. 2.
We assume that the known HIV-infected persons are infectious, but at a much lower rate than those that do not know they are infected due to education or change of behavior. 3. We assume that once a person develops AIDS it is no longer infectious.

We consider the following two cases: Case 1. If λ ′γ ′ = 0, the system has a disease-free equilibrium P 0 = (0, 0, 0) if σ +k 1 = λ-γ = 0 and a set of endemic equilibria

P * = (X * , Y * , Z * ) of the form (λ ′ -k 2 X * )Y * + (σ + k 3 X * )X * = 0 and Z * = βX * +β ′ Y * µ ′ if σ + k 1 = λ -γ = 0. Case 2. If λ ′ -γ ′ =
0, the system has two equilibria, P 0 = (0, 0, 0) is the diseasefree equilibrium and

P * = (X * , Y * , Z * ) is the endemic equilibrium. Set R 1 = λ ′ γ ′ , R 2 = λ γ and R 0 = k 1 k 1 + γ R 1 + γ k 1 + γ R 2 . (2)
Then

X * = γ ′ (k 1 + γ)(R 0 -1) γk 2 (R 2 -1) + γ ′ k 3 (1 -R 1 ) , Y * = X * γγ ′ R 2 -1 1 -R 1 , Z * = βX * + β ′ Y * µ ′ .
The endemic equilibrium is feasible if and only if R 0 -1, R 2 -1 and 1 -R 1 have the same sign. The parameters R 0 , R 1 and R 2 play a significant role in the analysis of the behaviour of trajectories for (1). Since Case 1 requires specific values of the parameters and hence is of little practical importance, we will suppose that λ ′γ ′ = 0. The analysis of the stability of the equilibria points and the asymptotic behaviour of the trajectories is given in Table 1.

Table 1: Asymptotic states for the model. "GAS" denotes equilibrium is globally asymptotically stable, "LAS" denotes it is locally asymptotically stable and "N E" denotes equilibrium does not exist in the domain

. R0 R1 R2 P0 P * (X, Y, Z) → < 1 < 1 < 1 GAS NE P0 < 1 < 1 > 1 GAS NE P0 < 1 > 1 < 1 LAS unstable P0 or unbounded > 1 < 1 > 1 unstable GAS P * > 1 > 1 < 1 unstable NE unbounded > 1 > 1 > 1 unstable NE unbounded

Conclusions

In this work, we have identified three threshold parameters of epidemiological importance, namely R 0 , R 1 , and R 2 , for the dynamic behavior of the system in question. A summary of their respective roles in the asymptotic state of the system is given in Table 1. Previous modeling studies of HIV with secondary reproduction numbers includes [START_REF] Hsieh | Behavior Change and Treatment of Core Group and Bridge Population: Its Effect on the Spread of HIV/AIDS[END_REF], [START_REF] Hsieh | A Class of Models for HIV Contact Tracing in Cuba: Implications for Intervention and Treatment. Deterministic and Stochastic Models for AIDS Epidemics and HIV Infection with Interventions[END_REF] and [START_REF] Hsieh | HIV Model with Secondary Contact Tracing: Impact of the Partner Notification Program in Cuba[END_REF]. The biological interpretations of these parameters is the following. R 0 gives the number of infections by an infective who is detected via random screening; R 1 is the number of infections caused by an infective after he/she has been tested positive for HIV; and R 2 is the number of infections caused by an infective that is not detected during asymptomatic period, i.e., if an infective either develops AIDS-related illness or pass away before progression to AIDS. We note that R 0 can be considered as the basic reproduction number, if we do not consider the detection of HIV-positive individuals in Cuba as part of disease surveillance instead of an intervention measure, especially in light of the realistic model assumption that those detected to be HIV-positive can still infected others, albeit at a lower level, before the onset of AIDS-defined illnesses [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF]. Moreover, R 0 determines whether the DFE (P 0 ) is locally asymptotically stable, as is the typical role of a basic reproduction number. R 1 and R 2 are secondary reproduction numbers, which help us determine the asymptotic behavior of the system, often relating to the endemic equilibrium.

In the present model, whether R 1 is larger than unity determines the existence of an unstable endemic equilibrium

P * when R 0 < 1. When R 0 > 1, R 1 < 1
in combination with R 2 > 1 ensure the existence of a globally asymptotically stable equilibrium P * . Serving as threshold parameters for the model system, these reproduction numbers have an obvious and epidemiologically meaningful interpretation. For the proposed model of detection and surveillance of HIV epidemic in Cuba, our results indicate that random screening is most important as a mean of surveillance, since the number of infections due to an infective detected through random determines whether the DFE is (locally) asymptotically stable.

In other words, if the averaged total number of infections by an infective detected through random screening exceeds one, then there will always be an epidemic.

On the other hand, local asymptotic stability for the DFE can be achieved (i.e., R 0 can be brought down to less than one) if: (i) random screening is sufficiently effective (k 1 large), and (ii) new infections by detected HIV-positive individuals is minimal (λ ′ small). Our results further highlight the importance of the education for the known infectious, in light of the second threshold parameter R 1 .

If the known infectious do not change, through education, their behavior and continue to practice risky sexual habits and their infection rate is still high (i.e., λ ′ sufficiently large) so that the average number of infections by a known infective (R 1 ) exceeds unity, then the endemic equilibrium is always unstable and there is always a possibility for the total number of infectious (i.e., X + Y ) to increase without bound. This is true even when R 0 < 1, provided that the initial population sizes are outside the domain of attraction of the DFE. This scenario of adverse impact of public health measures, which had been shown previously to be theoretically possible in [START_REF] Anderson | Potential of community-wide chemotherapy or immunotherapy to control the spread of HIV-1[END_REF], [START_REF] Hsu | Modeling intervention measures and public response during SARS outbreak[END_REF], is only possible if k 1 is sufficiently small compared to γ, knowing that R 0 is a convex combination of R 1 and R 2 (see 2). In other words, an ill-designed detection system might adversely lead to the epidemic increasing without bound if (i) random screening is not comprehensive enough (k 1 too small); (ii) lack of an education program to change behavior (λ ′ too high); and (iii) the prevalence is too high when the system is first implemented (initial infective populations outside of the domain of attraction of the DFE). This result further highlights the importance of universal testing in highprevalence regions [START_REF] Granich | Universal voluntary HIV testing with immediate antiretroviral therapy as a strategy for elimination of HIV transmission: a mathematical model[END_REF], [START_REF] Hsieh | Correspondence to "Universal voluntary HIV testing and immediate antiretroviral therapy[END_REF]. On the other hand, if through an adequate education to change the behavior of the known infectious (so that λ ′ is sufficiently low) the average number of infections by a known infective is less than one, then either the DFE P 0 or the endemic equilibrium P * is globally asymptotically stable, leading to a more manageable epidemic for the public health purposes, even if the disease is not eradicated.

Application to the Cuban HIV/AIDS data.

The modeling results are clearly relevant to our understanding of the current state of the HIV epidemic in Cuba. We will use the model [START_REF] Anderson | Potential of community-wide chemotherapy or immunotherapy to control the spread of HIV-1[END_REF] to fit the data of the known HIV positives and AIDS cases in Cuba. We have divided the period 1986-2008 into two different time periods, namely 1986-1999 and 1999-2008, to take into account the introduction of the family doctors in the detection system.

In each period we have parameters that can be estimated from the data. These are β, β ′ , µ, µ ′ and X(0). Y (0) and Z(0) are known. In previous work [START_REF] De Arazoza | A non linear model for a sexually transmitted disease with contact tracing[END_REF], [START_REF] Hsieh | HIV Model with Secondary Contact Tracing: Impact of the Partner Notification Program in Cuba[END_REF] estimates for λ and λ ′ have been obtained. The family doctors program had actually started after 1990, but only as a pilot project where the family doctors typically did not prescribe HIV testing. It is only after 1999 that detection through the family doctors started to take on an important role in the yearly detection figures, arriving at more than 30% of new detections in a year. Another significant difference between the two periods chosen is λ ′ . In the first years (1986-1999) the sanatorial system played an important role in preventing HIV transmission from persons that had been detected. Hence we suppose that there is practically no transmission from the persons living with HIV (Y ) and assume λ ′ = 0 for the first period. We fit the model to the data to obtain values for k 1 , k 2 and k 3 by minimizing an error function. As traditional optimization methods failed to work properly we used a genetic algorithm approach to find an initial point for starting the optimization method using a gradient method. To compute standard errors for the parameters, 200 fitting runs were made using different values of the known parameters, taken randomly from their confidence interval.

Using PET a software written on MATLAB, we obtain the least-square estimates for the unknown model parameters (k 1 , k 2 and k 3 ), for each of the two periods of the Cuban HIV epidemic from 1986 to 2008, by fitting the Cuban HIV data of the persons known to live with HIV to the model as described previously. Here the resulting numbers for 1999 obtained from the first stage of estimation using data from 1986-1999 were used as initial values for the second stage of estimation using the 1999-2008 data. The estimated mean values of the parameters k 1 , k 2 , and k 3 for each of the two periods with 95% confidence intervals, obtained from the 100 best fits, are given in Table 2. We can also compute the theoretical values of the number of the unknown persons living with HIV, X(t), from the estimation results. For the unknown X we have between 2300 and 2400 at the end of 2008. By comparing the estimated results for the two periods, we conclude that detection by random screening (k 1 ) improved significantly after 1999, perhaps reflecting the steeper increase in reported cases after 2000 [START_REF] De Arazoza | The HIV/AIDS epidemic in Cuba: Description and tentative explanation of its low HIV prevalence[END_REF], while detection via contact tracing (k 2 ) was at a similar level throughout the whole course of the epidemic. Detection by family doctors (k 3 ) was slightly higher than that of contact tracing after 1999 but of similar magnitude. Both the analytical result (of the dynamics) and the datafitting parameter estimates indicate that random screening was the most effective route of detection, while contact tracing and family doctors played mainly secondary roles, as had been previously proposed in [START_REF] Hsieh | Correspondence to "Universal voluntary HIV testing and immediate antiretroviral therapy[END_REF]. The estimated values of the parameters in the model also allow us to calculate the three reproduction numbers, R 0 , R 1 , and R 2 , with the 95% confidence intervals for each of the two time periods which are given in Table 3 Table 3: Estimated mean values with the 95% confidence intervals for R 0 , R 1 , and R 2 . UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals. As we can see, for both time periods the Cuba HIV epidemic is in the case of R 0 > 1, R 1 < 1 and R 2 > 1. P 0 is unstable and P * is globally asymptotically stable with the trajectories approaching P * asymptotically. Hence we can conclude that the HIV epidemic in Cuba is approaching (in the long term) an endemic steady state which we can estimate from our parameters using the expression obtained from our modeling for P * = (X * , Y * , Z * ) given in Equation (2). That is, assuming no drastic changes in the prevention, transmission, detection, or treatment of HIV in Cuba in the long term future, there could be, eventually, around 2700 persons living with HIV that do not know they are infected. Further noting that the theoretical number of unknown persons living with HIV from the model is a little under 2400, we speculate with optimism that, at the endemic steady state, the number of persons living with HIV that represent the main core for the transmission of the epidemic in Cuba will not increase drastically, in the long term future.

The assumption that the parameters do not change in the long term is not a very real one, taking into account that every year new information is gained on the virus and that there is ample research on treatment, vaccines and other aspects that affect the dynamics of the epidemic. For example existing therapy reduces the probability of transmission for HIV, in terms of our model this means that coefficient λ ′ is reduced and this will make R 1 smaller, changing the value of the asymptotic point. But therapy delays the onset of AIDS, this means that β ′ gets smaller, also changing R 1 but making it bigger. So long term predictions (as in the case of asymptotic behavior) make sense only as an indicator of how is the epidemic going to behave. In this sense in table 1, we see that if R 1 > 1 trajectories could become unbounded, and this would mean that the epidemic is out of control. It is very important to manage the value of R 1 .

In this work we chose to use parameters that are constant or step functions (with

Table 2 :

 2 Estimated mean values with the 95% confidence intervals for parameters k 1 , k 2 and k 3 for both periods. UCI and LCI denote the respective upper and lower bounds for the 95% confidence intervals. .

	1986-1999 (k3 = 0)		1999-2008
	k1	k2	k1	k2	k3
	Mean 0.1347 2.381 × 10 -5	0.2195 2.728 × 10 -5 6.801 × 10 -5
	LCI 0.1345 2.366 × 10 -5	0.2192 2.725 × 10 -5 6.783 × 10 -5
	UCI 0.1349 2396 × 10 -5	0.2198 2.732 × 10 -5 6.818 × 10 -5

2 steps for now) that also produce a change in the model itself, changing the recruitment function from X to Y , F (X, Y ), from a linear polynomial in Y for the period 86-99, to a linear polynomial in X and Y for the period 99-2008.

Acknowledgments: This work was carried out during visits to the University of Paris Descartes by YHH and H de A. H de A, YHH and RL received support from the French "Agence National pour la Recherche" project "Viroscopy". H de A also received support as visiting professor from the Laboratoire Paul Painlevé of the University Lille 1, France. H. de A., A.S. and J.B. received support from the Spanish AECID, from their projects PCI D/023835/09 and D/030223/10. For all the support we have received from all the different sources we are grateful.