Hamza Hallak Elwan

Ramin Khayatzadeh

Julien Poette

Beatrice Cabon

Relative intensity noise in optical heterodyning applied to millimeter-wave systems



Abstract-In this paper, the impact of Relative Intensity Noise (RIN) on optical heterodyne millimeter-wave (MMW) generation for Radio-over-Fiber (RoF) applications and systems is presented. The theoretical and experimental study is provided for two different techniques that generate MMW signals. The model of RIN is also presented and is in close agreement with the experiment results. The novelty and originality of this paper is the analysis of RIN in MMW frequency band, where impairments of the MMW generated carrier due to RIN are shown for the first time. The theoretical concepts and the experimental results prove that RIN appearing as side bands of MMW signal is directly generated from initial RIN at low frequency. Two significant impacts are considered: the phase noise on the beat-note carrier and the frequency response of the system, in order to extract real RIN contribution.

Index Terms-Distributed Feedback (DFB) Laser, Millimeterwave (MMW), Passively Mode Locked Laser Diode (PMLLD), Phase Noise, Relative Intensity Noise (RIN), Radio-over-Fiber (RoF).

I. INTRODUCTION

ILLIMETER-WAVE (MMW) frequency band and beyond has driven a lot of research for future wireless communication applications using Radio-over-Fiber (RoF) to distribute MMW signals [START_REF] Guo | 60-GHz Millimeter-Wave Radio: Principle, Technology And New Results[END_REF]. MMW generation can use optical modes mixing on a photo-diode (PD) such as beating of two independent running lasers [START_REF] Stoehr | 60 GHz Radio-over-Fiber Technologies For Broadband Wireless Services[END_REF]. In this approach, the generated MMW signal exhibits a large phase noise because of the process of MMW generation and both laser linewidth and wavelength fluctuations [START_REF] Johansson | Generation And Transmission Of Millimeter-Wave Data Modulated Optical Signals Using An Optical Injection Phase-Lock Loop[END_REF]. Monolithic dual Distributed Feedback (DFB) lasers are developed for obtaining a narrow linewidth beat-note with reduced phase noise [START_REF] Van Dijk | Monolithic Dual Wavelength DFB Lasers For Narrow Linewidth Heterodyne Beat Note Generation[END_REF]. In another approach for optical heterodyning, Quantum Dashed Mode Locked Lasers (QD-MLL) is used to demonstrate MMW beatnote linewidth around 100 KHz. Heterodyning advantage is that mode beating avoids any external electrical oscillator [START_REF] Brendel | Analog Link Performance of Mode-Locked Laser Diodoes in the 60 GHz[END_REF].

In MMW frequency range, the phase noise and jitter were studied for stability problems of the generated signal and some solutions proposed to overcome this degradation whereas the Relative Intensity Noise (RIN) which is also an indicator of the laser intensity stability, has never been explored yet. Usually, RIN of laser is measured from DC to some GHz [START_REF] Poette | Highly-Sensitive Measurement Technique of Relative Intensity Noise And Laser Characterization[END_REF].

The authors are with the L'Institut de Microélectronique Electromagnétisme et Photonique -Laboratoire d'Hyperfréquences et de Caractérisation, IMEP-LaHC, Minatec, 3 Parvis Louis Néel -CS 50257 -38016 Grenoble Cedex 1, France (email : hallakeh@minatec.grenoble-inp.fr).

In this paper, RIN is studied for the first time close to MMW carrier as optical heterodyne signal. Both theoretical and experimental study confirm here that the RIN phenomenon in MMW frequency band, generated from initial RIN at low frequency. Two techniques using different lasers are shown, with a very good matching between the model of RIN and the experimental results. It is worth mentioning that the results of the presented study are applicable to any kind of optical frequency generation.

Depending on the heterodyne process, RIN at MMW signal can be larger than the initial RIN at low frequency and causes spectral degradation in optical heterodyne signals. For RoF communication systems in MMW range, RIN impact on generated MMW signals can be observed by measuring Error Vector Magnitude (EVM) after suppression of phase noise effect by using an envelope detector at the receiver [START_REF] Khayatzadeh | Impact Of Amplitude Noise In Millimeter-Wave Radio-over-Fiber Systems[END_REF]. In this paper, due to system stabilization and applications in MMW, we have used two techniques for optical heterodyne generation. The first technique uses two independent DFB lasers and is the simplest and cheapest way. The fundamental features of this technique are the ability to generate signals on different frequencies and the capability to demonstrate several cases of RIN close to beat-note by varying the temperature and the relative optical power of the DFB lasers. Phase noise is large and previous results [START_REF] Khayatzadeh | Accurate Millimeter-Wave Laser Phase Noise Measurement Technique[END_REF] are then used here to decorrelate RIN and phase noise impacts on MMW signal impairments. The second technique is based on Passively Mode Locked Laser Diode (PMLLD), where the Free Spectral Range (FSR) of an optical resonator at 60.64 GHz is used to generate the MMW signal. Using PMLLD, the frequency stability of MMW carrier is improved because phase noise is much smaller than previous technique. Results are compared.

Section II discusses the theory and experimental set-up. RIN results in [9 KHz -19 GHz] band is analyzed in section III, and in [50 -62] GHz range is given in section IV. Finally conclusion is presented in section V.

II. THEORY AND EXPERIMENTAL SET-UP

Generating MMW signal by optical heterodyning requires at least two optical signals. The electrical field E i of optical signal (i) can be expressed as:

      ) ( exp ) ( 1 ) ( t i t i j t i i A t i E       (1)
where A i , δ i ,  i , and  i are the amplitude of electrical field i, amplitude noise, angular frequency, and phase noise, respectively.

For each laser, RIN expression was described [START_REF] Agrawal | Optical Transmitter[END_REF], and for

Relative Intensity Noise in Optical Heterodyning Applied to Millimeter-Wave Systems

Hamza Hallak Elwan, Ramin Khayatzadeh, Julien Poette, and Beatrice Cabon M a standard semiconductor laser, the generic expression of RIN that can be extracted from laser rate equation is:

) ( 2 4 ) ( 2 2 ) ( ) (       i elec P PSD opt P opt i RIN    (2)
where δ opt is the optical power variation spectrum, P opt is the optical power, PSD is power spectral density of the intensity noise, P elec is DC electrical power, and δ i ()=FFT(δ i (t)). The combined signal after coupler is transmitted through conventional single mode optical fiber of 1.5 m length to PD at the receiver end. In this case, the beating between only two modes from independent lasers is considered, and they are not phase locked, the fiber length will just impact on the received optical power. The photodetected current I PD is expressed as:

      ))) ( 1 ) ( 2 ( ) 1 2 cos(( ) ( 2 ) ( 1 2 1 ))) ( 1 ) ( 2 ( ) 1 2 cos(( 2 1 2 ) ( 2 2 ) ( 1 1 2 1 ) ( t t t t RIN t RIN I I t t t I I t RIN I t RIN I I I t PD I                      (3)
where I i is the DC photodetected current due to optical signal (i, i=1 or 2), and RIN i is relative intensity noise RIN(t)=FFT -1 (RIN()). Second order terms related to noiseto-noise beatings have been neglected.

The first term of (3) represents the total average current, the second term represents the noise current of two optical signals from DC to some GHz. The third term is the beat-note current created during their simultaneous detection, where the beat frequency is equal to difference between frequencies of the two optical modes (f2-f1). The last term represents noise current of the two optical signals close to beat-note. The angular frequency ( 2 - 1 ) exhibits phase noise ( 2 - 1 ) that varies due to how the two optical modes are produced, which depends on the correlation between these modes. In the last term from (3), it can be noticed that there is interdependence between the two laser amplitudes in opposition to the second term. This non-linear process could lead to an increase of intensity noise during heterodyne process comparing to the initial intensity noise. This increase is noticeable especially when the optical powers of the two lasers are different.

After substituting (3) in (2), and as there is no correlation between modes when considering two independent DFB lasers, the initial relative intensity noise (RIN ini ) at low frequency and relative intensity noise close to beat-note (RIN mmw ) are then obtained as follows:

  2 2 1 ) ( 2 2 2 ) ( 1 2 1 ) ( I I f RIN I f RIN I f ini RIN    (4)     2 2 1 ) ( 2 ) ( 1 2 1 ) ( I I f RIN f RIN I I f mmw RIN    (5)
As can be seen, they are correlated. In the last term of (3), the product in time domain of mode intensity noises with the carrier leads to the frequency conversion of the noise spectrum directly to the beat frequency. In case the optical signals have the same amplitudes (I 1 =I 2 ), then RIN ini =RIN mmw . While a power difference exists, for example (I 2 >I 1 ), the RIN i of less powerful mode (RIN 1 ) in ( 5) is increased by a factor I 2 /I 1 as compared to RIN i in (4) thus RIN mmw is larger than RIN ini by a factor I 2 /I 1 . This phenomenon has not been published yet to the best of our knowledge and is evidenced in section III with both experimental and model results.

The paper investigates RIN at MMW, for two purposes: more accurate recognition between RIN ini and RIN mmw and better amplification of the received signal, we present results in [9 KHz-19 GHz] frequency band. After PD in Fig. 1 (path 1), the bias tee is added to suppress the DC component and to extract the photodetected RIN value. Low Noise Amplifier (LNA) whose center frequency is 9.5 GHz and gain is 40 dB is used. In path 2, MMW signal is amplified using 35 dB LNA of 50-62 GHz bandwidth, then the amplified signal is down-converted by a mixer 50-65 GHz of 6 dB losses and a local oscillator (LO). Finally, a spectrum analyzer (ESA) of 20 GHz bandwidth is monitored to extract spectral power density of measured RIN.

B. RIN based on PMLLD:

The electrical field E PMLLD of any multi-mode laser such PMLLD for (M) modes can be defined as:

           i t i t RF i t i i A t PMLLD E ) ( 0 exp ) ( 1 ) (     (6)
where i is the mode number, and A i , δ i , and  i are the amplitude, amplitude noise, and phase noise of mode i, respectively.  0 and  RF are the angular frequency of the first mode and the angular frequency of beat-note correspond to FSR. In Fig. 1, the generic schematic diagram of experimental set-up based on multi-mode laser is shown (second technique) for two frequency bands. In some cases, by increasing the bias current of optical source, PSD of thermal noise is close to PSD of intensity noise, then the sensitivity of the system must be improved by using a lock-in amplifier. A digital signal generator drives one electrode of a Mach-Zenhder modulator (MZM) and some synchronization is necessary between ESA and lock-in amplifier via signal generator. The output current of PD can be described as:

              ) ( ) ( cos : , ) ( ) ( : , ) ( ) ( cos 2 ) ( ) ( t i t j t i j j i j i t j RIN t i RIN j I i I j i j i t i t j t i j j I i I i t i RIN i I i i I t PD I                          (7)
The initial relative intensity noise (RIN ini ) and relative intensity noise at beat-note (RIN mmw ) can be expressed as:

2 2 ) ( ) (                      i i I i f i RIN i I f ini RIN (8)   2 2 : , ) ( ) ( ) (                        i i I j i j i t j RIN f i RIN j I i I f mmw RIN (9)
Terms in [START_REF] Khayatzadeh | Impact Of Amplitude Noise In Millimeter-Wave Radio-over-Fiber Systems[END_REF] for M modes have same expressions as in (3) for two modes. RIN expressions of ( 8) and ( 9) are also the same as ( 4) and ( 5) respectively, but the extremely important difference between the two cases is the correlation of modes for the case of PMLLD.

III. RIN IN [9 KHZ -19 GHZ] FREQUENCY BAND

The essential features of a two different DFB lasertechnique (see Fig. 1 For verifying (4) and ( 5), the bias current of the first laser DFB 1 is set to 15 mA (P opt = -2.1 dBm) while the bias current of the second laser DFB 2 is 90 mA (P opt = +3.5 dBm), thus RIN 1 and RIN 2 levels are very different for the first and second laser respectively. The experimental and theoretical results of initial RIN ini and beat-note RIN mmw are given in Fig. 2. On RIN in results, the thermal noise impact is removed when lasers are turned off, the shot noise has been measured and then removed, and the frequency response of the global system has been removed as well after global S parameters are measured with a Vector Network Analyzer (VNA).

In Fig. 2, the initial RIN ini from 9 KHz to 7 GHz and two side-bands RIN mmw on each side of the heterodyne signal are well observed. Since the stability of this technique is poor, the phase noise (red curve) is provided to prove that impairments are not due to phase noise, but come from laser intensity noise, the beat linewidth is 150 MHz. The values of initial RIN ini and the beat-note RIN mmw using ( 4) and ( 5) are calculated: RIN ini = -119.3 dB/Hz and RIN mmw = -112.7 dB/Hz; which well agree with the experimental results. From Fig. 2 and computations, it can also be observed that the beat-note RIN mmw is higher than the initial RIN ini by a factor I 2 /I 1 = 6.7 dB because the optical power of DFB 2 is larger , the RIN of less powerful mode (RIN 1 ) in ( 5) is then increased by a factor I 2 /I 1 compared to RIN 1 in (4) as expected.

The initial RIN ini has also been modeled (thin dashed curve) from RIN expression of laser [START_REF] Agrawal | Optical Transmitter[END_REF]. The phase noise of beat signal corresponds to the convolution of the two optical individual spectrums. The noise shape can be modeled by lorentzian distribution, while other random variations of the wavelength difference are considered as Gaussian fluctuations, resulting in a Voigt profile (thick dashed curve). The model of beat-note RIN mmw (light thick curve) in Fig. 2 In Fig. 3, we intend to show that RIN mmw and not phase noise impacts mainly the MMW impairments. Since the relative intensity noise (RIN 1 ) of DFB 1 has the main impact on the initial RIN ini , two different bias currents are set for DFB 1 : 13 mA (P opt = -4.1 dBm) and 17 mA (P opt = -0.6 dBm) corresponding to relaxation frequencies of 3.3 GHz and 4.6 GHz respectively, while DFB 2 is still biased at 90 mA (P opt = +3.5 dBm). The temperature of both lasers is set to keep the beat-note around 14 GHz, results are shown in Fig. 3. From this figure, it is inferred that RIN mmw close to beat signal is directly generated from initial RIN ini and have the same profiles. The difference between initial RIN ini and RIN mmw at beat-note for the above two cases is also observed because the lowest powered mode will be increased through the heterodyne process by I 2 /I 1 . Fig. 3 depicts that despite the phase noise of beat-note for both cases are approximately identical, RIN mmw keeps the same profiles as initial RIN ini , so this further confirms that RIN phenomenon at generated MMW signal is distinct from phase noise of heterodyne signal.

IV. RIN IN [50 -62] GHZ FREQUENCY BAND

In MMW measurements and using PMLLD to generate the MMW signal at 60.64 GHz due to FSR of an optical resonator, the schematic diagram of experimental set-up is described in path 2 of Fig. 1 (second technique). A mixer and an LO whose frequency is fixed at 52.2 GHz is used for down-conversion to adjust the ESA bandwidth. The mixing process does not influence on the quality of the involved signals because the LO phase noise is very low (-140 dBc/Hz for frequency offset more than 10 MHz) [START_REF] Agilent | D PSG Microwave Analog Signal Generator Datasheet[END_REF]. PMLLD is biased at 53 mA (P opt = +1.8 dBm), and the temperature is set to 25°C. Fig. 4 represents the beat-note in [50-62] GHz before downconversion, but PSD has been measured after mixer, thus in Fig. 4 using MATLAB the shifted spectrum is presented. Fig. 4 depicts one side-band RIN mmw of MMW signal because the insertion loss of system (PD, amplifier, mixer, and ESA) is -22 dB with frequency cutoff at 62 GHz. The initial RIN ini (black curve in a small frame) and phase noise (red curve) are presented, and the linewidth of MMW carrier is 90 MHz. The models of RIN ini (red curve in a small frame), RIN mmw (light curve), and phase noise (thick dashed curve) are a good agreement with the experimental results. According to equilibrium of optical power between M modes, it is very important to observe that the beat-note RIN mmw is not increased or decreased as compared to the initial RIN ini because of the factor I j /I i =1 as in ( 8) and ( 9), the value of RIN ini and RIN mmw is -130 dB/Hz. From this figure, the relaxation frequency of RIN ini and RIN mmw is approximately 5 GHz, so MMM RIN mmw is generated from initial RIN ini . As comparison between these two techniques using two different lasers and from the experimental results presented in Fig. 2 and Fig. 4, we can conclude that the physical phenomenon of RIN is the same: generation of RIN mmw close to beat signal for any kind of optical heterodyne generation. For two DFB lasers at MMW frequency band, we measured RIN mmw , and we obtained on the two side-bands of RIN mmw for beat carrier generated around 55 GHz due to controlling the temperature of both DFB lasers. We also measured RIN mmw when increasing the bias current of PMLLD using lock-in amplifier for improving the system sensitivity, but we can not present these results due to the limited size of the paper.

V. CONCLUSION

In this paper, we have analyzed relative intensity noise of optical heterodyne signals in MMW frequency band for RoF applications and systems, and we have verified the phenomenon of RIN close to beat-note based on two different techniques, two independent DFB lasers and PMLLD. In order to show a real RIN mmw , the phase noise result is provided, and the frequency response of system has been removed.

The experimental results confirm that RIN mmw at MMW carrier is directly generated from initial RIN ini at low frequency and is distinct from phase noise of beat signal. Due to the heterodyne process the RIN mmw level can be higher than the RIN ini , and this increasing depends on a factor I j /I i between the modes power. These results are independent from how the two optical modes are generated, so they can be directly transposed to any kind of optical generation of beat frequencies. The models of RIN ini and RIN mmw are presented and are very good matching with the experimental results.

A

  . RIN based on two DFB lasers: Two independent DFB lasers, DFB 1 and DFB 2 operating in the C-band, are used to generate beat-note carrier which is adjusted and controlled by temperature and bias current of the two lasers. The polarization controller (PC) is utilized for optimizing the beating power, and (50/50) coupler is used to mix two optical signals on PD of 70 GHz bandwidth. The generic schematic diagram of experimental set-up for two separate optical sources (first technique) is shown in Fig. 1.

Fig. 1 .

 1 Fig. 1. Experimental set-up for optical heterodyning based on two independent lasers and PMLLD in [9 KHz-19 GHz] and [50-62] GHz.

  path 1) are ability to distribute beat signals on various frequencies and demonstration different cases of RIN at beat-note by changing the temperature and the relative optical power of lasers. Due to the relaxation frequency of the lasers of 4.1 GHz and the 19 GHz bandwidth limitation of the amplifier, a beat signal around 14 GHz has been chosen to highlight the impact of RIN in optical heterodyne process as shown on Fig. 2.

Fig. 2 :

 2 Fig. 2: RIN and phase noise of a heterodyne signal around 14 GHz for two independent DFB lasers.

  refers to the last term in (3), representing the phase noise of the beating signal multiplied in time domain by the sum of the individual laser intensity noises. Then, the model results of beat signal phase noise and RIN mmw are summed and presented in the blue curve. Based on the aforementioned, experimental results are well confirmed by model results and give evidence of RIN phenomenon at MMW, as presented in (3), (4), and (5).

Fig. 3 :

 3 Fig. 3: RIN of a beat-note signal around 14 GHz for two different bias current of the first laser DFB1.

Fig. 4 :

 4 Fig. 4: RIN and phase noise of a heterodyne signal around 55 GHz for two independent DFB lasers.

ACKNOWLEDGMENT

We would like to express our gratitude to III-V labs for proving the PMLLD chip.


A D Chepelianskii

F Chevy

E Raphaël

On capillary-gravity waves generated by a slow moving object

Keywords: numbers: 47, 35, -i , 68, 03, -g

We investigate theoretically and experimentally the capillary-gravity waves created by a small object moving steadily at the water-air interface along a circular trajectory. It is well established that, for straight uniform motion, no steady waves appear at velocities below the minimum phase velocity cmin = 23 cm • s -1 . We show theoretically that no such velocity threshold exists for a steady circular motion, for which, even for small velocities, a finite wave drag is experienced by the object. This wave drag originates from the emission of a spiral-like wave pattern. Our results are in good agreement with direct experimental observations of the wave pattern created by a circularly moving needle in contact with water. Our study leads to new insights into the problem of animal locomotion at the water-air interface.

Capillary-gravity waves propagating at the free surface of a liquid are driven by a balance between the liquid inertia and its tendency, under the action of gravity and surface tension forces, to return to a state of stable equilibrium [START_REF] Landau | Lifshitz Fluid Mechanics[END_REF]. For an inviscid liquid of infinite depth, the dispersion relation relating the angular frequency ω to the wave number k is given by ω 2 = gk + γk 3 /ρ, where ρ is the liquid density, γ the liquid-air surface tension, and g the acceleration due to gravity [START_REF] Acheson | Elementary Fluid Dynamics[END_REF]. The above equation may also be written as a dependence of the wave velocity c(k) = ω(k)/k on wave number: c(k) = (g/k + γk/ρ) 1/2 . The dispersive nature of capillary-gravity waves is responsible for the complicated wave pattern generated at the free surface of a still liquid by a moving disturbance such as a partially immersed object (e.g. a boat or an insect) or an external surface pressure source [START_REF] Acheson | Elementary Fluid Dynamics[END_REF][START_REF] Lighthill | Waves in Fluids[END_REF][START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Rayleigh | [END_REF]6]. Since the disturbance expends a power to generate these waves, it will experience a drag, R w , called the wave resistance [START_REF] Lighthill | Waves in Fluids[END_REF]. In the case of boats and large ships, this drag is known to be a major source of resistance and important efforts have been devoted to the design of hulls minimizing it [7]. The case of objects small relative to the capillary length κ -1 = (γ/(ρg))

1/2 has only recently been considered [8,9,10,11].

In the case of a disturbance moving at constant velocity V , the wave resistance R w cancels out for V < c min where V stands for the magnitude of the velocity, and c min = (4gγ/ρ) 1/4 is the minimum of the wave velocity c(k) given above for capillarity gravity waves [START_REF] Lighthill | Waves in Fluids[END_REF][START_REF] Lamb | Hydrodynamics[END_REF]8]. For water with γ = 73 mN • m -1 and ρ = 10 3 kg • m -3 , one has c min = 0.23 m • s -1 (room temperature). This striking behavior of R w around c min is similar to the well-known Cerenkov radiation emitted by a charged particle [12], and has been recently studied experimentally [13,14]. In this letter, we demonstrate that just like accelerated charged particles radiate electromagnetic waves even while moving slower than the speed of light [START_REF] Jackson | Classical Electrodynamics[END_REF], an accelerated disturbance experiences a non-zero wave resistance R w even when propagating below c min . We consider the special case of a uniform circular trajectory, a situation of particular importance for the study of whirligig beetles (Gyrinidae, [START_REF] Nachtigall | The Physiology of Insecta[END_REF]) whose characteristic circular motion might facilitate the emission of surface waves that they are thought to be used for echolocation [START_REF] Tucker | [END_REF][START_REF] Denny | Air and Water[END_REF]. This work is therefore restricted to the effect of a wake stationary in the rotating frame, and do not consider time dependent contributions, like vortex shedding [START_REF] Bush | [END_REF]29].

We consider the case of an incompressible infinitely deep liquid whose free surface is unlimited. In the absence of external perturbation, the free surface is flat and each of its points can be described by a radius vector r = (x, y) in the horizontal plane. The motion of a small object along the free surface disturbs the equilibrium position of the fluid, and each point of the free surface acquires a finite vertical displacement ζ(r). Rather than solving the complex hydrodynamic problem of finding the flow around a moving object, we consider the displacement of an external pressure source P ext (r, t) [START_REF] Rayleigh | [END_REF]6]. The equations of motion can then be linearized in the limit of small wave amplitudes [START_REF] Dias | For nonlinear effects[END_REF].

In the frame of this linear-response theory, it is convenient to introduce the Fourier transforms of the pressure source Pext (k, t) and of the vertical displacement ζ(k, t) [START_REF]The Fourier transform f (k, t) is related with the function f (r, t) through f (r, t) = d 2 k (2π) 2 e ik.r f (k, t)[END_REF]. It can be shown that, in the limit of small kinematic viscosity ν, the relation between ζ(k, t) and Pext (k, t) is given by [9]

∂ 2 ζ ∂t 2 + 4 ν k 2 ∂ ζ ∂t + ω 2 (k) ζ = - k Pext (k, t) ρ (1)
In this letter we assume that the pressure source has radial symmetry and that the trajectory r 0 (t) of the object is circular, namely : r 0 (t) = R (cos(Ωt), sin(Ωt)).

Here R is the circle radius, and Ω is the angular frequency. The linear velocity of the object is then given by V = R Ω. With these assumptions, the external pressure field is P ext (r, t) = P ext (|rr 0 (t)|, yielding in Fourier space Pext (k, t) = Pext (k)e -ik.r0(t) . Since the right hand side of Eq. ( 1) is periodic with frequency Ω, it is possible to find its steady state solution by expanding the right hand side into Fourier series. The problem then becomes equivalent to the response of a damped oscillator to a sum of periodic forces with frequencies nΩ, where n is an integer. The vertical deformation at any time t can then be reconstructed by evaluating the inverse Fourier transform. For the particular case of uniform circular motion, the time dependence is rather simple. Indeed, in steady state, the deformation profile rotates with the same frequency Ω as the disturbance. Therefore, in the rotating frame, ζ depends on the position r only. The analytical expression of ζ(r) in cylindrical coordinates (x, y) = r(cos φ, sin φ) is given by

ζ(r, φ) = ∞ n=-∞ e inφ k 2 dk 2πρ Pext (k)J n (kr)J n (kR) n 2 Ω 2 -ω 2 (k) + 4inνk 2 Ω (2)
where J n is n-th order Bessel function of the first kind.

The summation index n is directly related to the n-th Fourier harmonic of the periodic function e -ik.r0(t) and, since the problem is linear, the contributions of all the harmonics add together.

The knowledge of the exact structure of the wave pattern is precious, but a quantitative measurement of the wave resistance is needed in order to understand, for example, the forces developed by small animals moving at the surface of water. In the case of the circular motion under study, the wave resistance R w can be calculated from its average power P w = -d 2 r P ext (r, t) ∂ζ(r,t) ∂t by R w = P w /V . Using the Fourier expansion of ζ, one then obtains in the limit νκ/c min → 0 (for water, νκ/c min ∼ 10 -3 ):

R w (V, R) = n>0 n ρR (k n J n (k n R) Pext (k n )) 2 d ω 2 dk kn (3)
where k n is the unique solution of the equation ω(k n ) = n Ω (the notation R w (V, R) stresses the dependence of R w on the velocity magnitude and on the trajectory radius). Equation [START_REF] Lighthill | Waves in Fluids[END_REF] shows that the wave resistance R w takes the form of a sum R w = n>0 A n , where the A n are positive numbers that measure the contribution of each Fourier mode of the external pressure source (with frequency n Ω) to the wave resistance.

A numerical calculation of the wave resistance is presented in Fig. 1 for a pressure source Pext (k) = p 0 exp(-kb), where p 0 is the total force exerted on the surface and b is the typical object size [22]. As observed, Eq. 3 differs significantly from the original prediction on the wave drag in the case of a straight uniform motion with velocity V [8,11] given by

R w,l (V ) = ∞ 0 kdk 2πρ P 2 ext (k) θ(V -c(k)) V 2 1 -(c(k)/V ) 2 , (4)
where θ(.) is the Heavyside function and c(k) = ω(k)/k is the phase velocity. Most notably, the wave drag for a circular motion is non-zero for all velocities, even for V < c min where wave-resistance vanishes exactly in the case of a linear motion and this effect is far from negligible: for R/b = 10 and at velocities as slow as V /c min ∼ 0.6, the wave drag is still one fifth of that applied to an object moving linearly at V /c min = 1. The radiation of waves by an accelerated particle should not be surprising and actually is a very general phenomenon that can be observed for instance in electromagnetism (bremsstrahlung) or in general relativity (Zeldovich-Starobinsky effect [23]). Mathematically, the fact that, for a circular motion, the wave resistance is finite even below c min can be understood as follows. In the case of uniform motion, all the wavenumbers such as c(k) < V contribute to the wave drag, whereas for circular motion this is the case for only a discrete set of wavenumbers k n . While the condition c(k) < V can be satisfied only when V > c min , the equations for the wavenumber k n , ω(k n ) = nV /R, have positive solutions for any velocity V . These wavenumbers k n create finite contributions A n > 0 to the wave drag. Therefore for a circular trajectory a finite wave drag exists at any velocity V > 0; for the same reasons R w is also continuous at V = c min . Moreover, the wave resistance develops a small oscillating component as a function of the velocity V . It originates from the oscillatory behavior of Bessel func- tions and will be analyzed more thoroughly in a future publication. Finally, we note that despite these striking differences Eqn. ( 3) and (4) should coincide in the limit of a large trajectory radius R. We confirmed this behavior by checking both analytically [24] and numerically that in the limit R → ∞, R w (V, R) → R w,l (V ). However even if the circular wave drag R w (V, R) is close to R w,l (V ) starting from R/b ∼ 10, important differences remain even up to very large values of R/b such as R/b ∼ 100. Figure 2 represents the wave crest pattern (computed numerically form Eq.( 2)) at the origin of this finite wave drag. It exhibits characteristic concentric Archimedean spirals (also known as arithmetic spirals) of the form r = aφ+ r 0 . This can be understood from our theoretical results as follows. In a first estimation, one can assume that the integrals in equation Eq. ( 2) are dominated by the contribution of the poles at k = k n . Thus ζ(r) can be written as ζ(r) ∼ 1 √ r n B n e i(nφ-knr) , where we have used the asymptotic development of J n (k n r) at large distances r and B n are complex coefficients that do not depend on the position r = r(cos φ, sin φ). By separating the contribution of the different modes in the relation F (t) = -d 2 rP ext (r, t)∇ζ(r, t), one finds that B n is proportional to A n (where, as defined earlier, the positive coefficients A n measure the contribution of each Fourier mode to the wave drag: R w = n>0 A n ). One can show that in the regime of small object sizes κb ≪ 1, the proportionality constant between B n and A n depends only weakly on the Fourier mode number n; thus, one has -knr) . We have checked numerically that in the regime V < c min , the distribution of the coefficients A n is usually peaked around n ∼ κR. For example, for κR = 10 and κb = 0.1, A n is peaked around n = 10 for velocities V in the interval (c min /2, c min ). The wave-crests are given by the lines of constant phase nφ -k n r = const of the dominant mode n = κR, leading to the following expression for a:

ζ(r) ∝ 1 √ r n A n e i(nφ
a ≈ κR k(ω = κV ) (5)
where k(ω) is the inverse function of ω(k). An interesting special case of the formula Eq. ( 5) corresponds to V = c min , for which one obtains a ≈ R. The spiral predicted by Eq. ( 5) is in very good agreement with the exact numerical results (Eq. ( 2)), as can be seen in Fig. 2.

We have also compared our theoretical approach with experimental results obtained using a one millimeter wide stainless steel needle immersed in a 38 cm wide water bucket. The needle was rotated on circular trajectories of various radii and angular velocities. Since direct measurement of wave drag, and in particular comparison with theory, is non-trivial even for a linear motion [13,14], we restricted ourselves to the study of the wake itself. A typical wave pattern obtained by this method is shown on Fig. 3 for R ≈ 2.7cm and Ω ≈ 2π × 1.2 Hz (corresponding to V /c min ≈ 0.9) and unambiguously demonstrates the existence of a wake at velocities smaller than c min . The observed wave pattern is in remarkable agreement with the theoretical prediction r = aφ + b with a given by Eq. ( 5) and r 0 a free parameter corresponding to an overall rotation of the spiral [25]. For V /c min lower than 0.8, no wake was observed by naked eye. At lower rotation velocities, we probed the surface deformation by measuring the deflexion of a laser beam reflected by the air-water interface at a distance r = 11 cm from the rotation axis.

Using this scheme, we have established the existence of waves down to V /c min ≈ 0.6, and verified quantitatively that the wave packet spectrum is peaked around < ω >∼ κR Ω (see Fig. 3). Experimentally, the frequency < ω > corresponds to the period of the fast temporal oscillations of the laser deflection angle (see Fig. 3 inset). In order to compare our experimental results with our model, we note that the deflection of the laser at a point r is proportional to the derivatives 1 r ∂ζ(r,t) ∂φ and ∂ζ(r,t) ∂r .

For simplicity, we will mainly consider the angular derivative, but we have checked numerically that our result do not depend on this choice. Using Eq. ( 2) the angular derivative can be decomposed into Fourier series: ∂ζ(r,t) ∂φ = n C n e in(φ-Ωt)-iknr . The coefficients C n are proportional to the contribution of the frequency nΩ to the wave packet spectrum and we can thus calculate the mean wave packet frequency using the expression:

< ω >= Ω n>0 n|C n |/ n>0 |C n |.
As shown in Fig. 3, our model is consistent with good accurcy with the experimental data without any adjustable parameters.

Below V /c min ≈ 0.6, the signal to noise ratio of the experiment becomes to small to observe the laser deflection. Note that this value is in qualitative agreement with Fig. 1 where the wave resistance (hence the wave amplitude) has also significatively decreased with respect to its maximum value for V /c min 0.5: we indeed note that for To summarize, we have shown theoretically that a disturbance moving along a circular trajectory experienced a wave drag even at angular velocities corresponding to V < c min , where c min is the minimum phase velocity of capillary-gravity waves. Our prediction is supported by experimental observation of a long distance wake for V /c min as low as 0.6. For V /c min > 0.8, we observed by naked eye Archimedean spiral shaped crests, in good agreement with theory. These results are directly related to the accelerated nature of the circular motion, and thus do not contradict the commonly accepted threshold V = c min that is only valid for a rectilinear uniform motion, an assumption often overlooked in the literature. It would be very interesting to know if whirligig beetles can take advantage of such spirals for echolocation purposes. Although restricted to stationary wakes and thus excluded effects such as vortex shedding, the results presented in this letter should be important for a better understanding of the propulsion of water-walking insects [START_REF] Alexander | Principle of Animal Locomotion[END_REF][START_REF] Bush | [END_REF]28,29] where accelerated motions frequently occurs (e.g when hunting a prey or escaping a predator [30]). Even in the case where the insect motion is rectilinear and uniform, one has to keep in mind that the rapid leg strokes are accelerated and might produce a wave drag even below c min .

FIG. 1 :

 1 FIG. 1: (Color online) Plot of the wave resistance Rw in units of p 2 0 κ/γ, as a function the reduced velocity V /cmin = R Ω/cmin for different ratios between the trajectory radius R, and the object size b, as predicted by Eq. (3). The red curve (presenting many oscillations) corresponds to R/b = 100, while the black one (with fewer oscillations) corresponds to R/b = 10. The green curve displaying a typical discontinuity at V = cmin is the wave drag for a straight uniform motion with velocity V [8]. The object size, b, was set to b = 0.1 κ -1 .

FIG. 2 :

 2 FIG. 2: (color online) Wave radiation for V ≈ 21cm /s ≈ 0.9cmin with a radius R ≈ 2.7cm ≈ 9κ -1 Left: Color diagram of the surface deformation ζ(r) computed numerically from Eq. (2). This image represents a square region of size 400κ -1 around the center of rotation, red color corresponds to maximal ζ(r) values, while green corresponds to minimal values of ζ(r). The cross indicates the center of the trajectory and the moving object is located in the region of highest deformation. Right: Photography of the wave crests generated on a water surface by a needle rotating at a velocity . On both pictures, the black curve represents the the Archimedean spiral of radius given by Eq. (5.

FIG. 3 :

 3 FIG.3:(Color online) Inset: typical time dependence of the laser deflection angle (arbitrary units) during a rotation period Trot = 2π/Ω, the fast oscillation frequency is given by < ω >= 2π/T . Main figure: Dependence of the ratio < ω > /Ω on κR for different needle velocities. The dashed curves represent experimental results, while the continuous curve display the numerical results of our model. Red, green and blue curves (diamonds, squares and circles respectively) correspond to V /cmin = 0.69, 0.76 and 0.84. The black curve correspond to the analytical estimate < ω > /Ω = κR.

We are grateful to José Bico, Jérôme Casas, M. W. Denny and J. Keller for fruitful discussions. F.C. acknowledges support from Région Ile de France (IFRAF) and A.C. acknowledges support from Ecole Normale Supérieure Paris.


P Biasutti
email: pierre.biasutti@labri.fr

A Bugeau

J-F Aujol

M Brédif

Visibility estimation in point clouds with variable density

Keywords: 3D point cloud, visibility, visualization, LiDAR, dataset, benchmark

Estimating visibility in point clouds has many applications such as visualization, surface reconstruction and scene analysis through fusion of LiDAR point clouds and images. However, most current works rely on methods that require strong assumptions on the point cloud density, which are not valid for LiDAR point clouds acquired from mobile mapping systems, leading to low quality of point visibility estimations. This work presents a novel approach for the estimation of the visibility of a point cloud from a viewpoint. The method is designed to be fully automatic and it makes no assumption on the point cloud density. The visibility of each point is estimated by considering its screen-space neighborhood from the given viewpoint. Our results show that our approach succeeds better in estimating the visibility on real-world data acquired using LiDAR scanners. We evaluate our approach by comparing its results to a new manually annotated dataset, which we make available online.

Introduction

Over the past decade, the use of 3D point clouds as an alternative to meshes has been constantly growing. A 3D point cloud simply consists in 3D positions, sometimes associated with supplementary information such as its color, reflectance or normal. It can be considered as a sampling of continuous surface, providing a simpler representation than the full topology.

The estimation of the visibility of a point cloud consists in assigning a label to each point of the scene: visible if the point lies on an object that is directly visible from a given viewpoint, non-visible otherwise (Fig. 1). This task is a typical step for various applications in computer graphics such as in surface reconstruction [START_REF] Zach | A globally optimal algorithm for robust TV-L 1 range image integration[END_REF][START_REF] Shalom | Cone carving for surface reconstruction[END_REF][START_REF] Berger | A survey of surface reconstruction from point clouds[END_REF] in which estimating and removing points that are not visible from a given point of view improves the interpolation and the approximation of the surface to recover. In point cloud rendering and visualization [START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF][START_REF] Bouchiba | High quality and efficient direct rendering of massive real-world point clouds[END_REF], the estimation of the visibility enables better rendering performances as well as an improvement of the scene understanding. For both application scenarios, the existing methods [START_REF] Zach | A globally optimal algorithm for robust TV-L 1 range image integration[END_REF][START_REF] Shalom | Cone carving for surface reconstruction[END_REF][START_REF] Berger | A survey of surface reconstruction from point clouds[END_REF][START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF]a. b. c. d. [START_REF] Bouchiba | High quality and efficient direct rendering of massive real-world point clouds[END_REF] strongly rely on strict sampling assumptions [START_REF] Berger | A survey of surface reconstruction from point clouds[END_REF]) (e.g. on point clouds with constant density in terms of number of points per cubic meters).

Recently, the development of acquisition systems designed for acquiring urban scenes has been increasing. In particular, Mobile Mapping Systems (MMS) equipped with LiDAR (Light Detection And Ranging) [START_REF] Paparoditis | Stereopolis II: A multi-purpose and multi-sensor 3D mobile mapping system for street visualisation and 3D metrology[END_REF][START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF][START_REF] Maddern | 1 year, 1000 km: The Oxford RobotCar dataset[END_REF] have been widely used to survey cities, road, highways, etc. Those campaigns have resulted in the production of large, unorganized point clouds that provide precise 3D representations of the urban environment. Due to the acquisition method, these point clouds present high variation in their sampling and density. Most modern MMS also embed additional materials, mostly optical images, which provide complementary information on the scene. The multimodal aspect of these datasets may be leveraged to improve detection, classification and prediction techniques in urban environments [START_REF] Benenson | Ten years of pedestrian detection, what have we learned[END_REF][START_REF] Eigen | Depth map prediction from a single image using a multi-scale deep network[END_REF]. Therefore, the fusion and the registration of LiDAR and optical data became critical as the use of multi-modal data definitely increases performances of classification/prediction algorithms. Most of the recent related works strongly rely on good visibility estimates [START_REF] Mastin | Automatic registration of LIDAR and optical images of urban scenes[END_REF][START_REF] Guislain | Fine scale image registration in largescale urban LIDAR point sets[END_REF].

The majority of actual LiDAR/optical registration techniques that use visibility rely on estimation techniques that were built for point clouds with strict sampling assumptions that are not met by the LiDAR data on which they operate. On the other hand, point cloud rendering and surface reconstruction methods presented above are not designed to perform on point clouds with variable density. However, the quality of the visibility estimation is a crucial preprocessing step for multi-modal fusion applications as it drastically lowers the ambiguities from one modality to another. This work aims at studying a new method for estimating the visibility in point clouds without constant density acquired via MMS to improve the data fusion.

The paper contribution is threefold: first, we propose a novel approach for the visibility estimation in a point cloud that is robust to high density variations. This method is designed to be fully automatic and to perform well on any types of 3D point clouds. The second contribution of this article is a new visibility dataset of over 1 million annotated points for testing the performances of visibility estimation methods. This dataset, as well as the code for our method, are made publicly available online. The third contribution of this dataset is a full numerical and visual comparison between our method and other state-of-the-art methods.

The paper is organized as follows: first, a brief overview of the related work is presented. Then, the methodology of the method is explained. Finally, evaluation and results are shown and a conclusion is drawn.

Related works

There has been many contributions to the state-ofthe-art techniques for the estimation of the visibility of a point cloud given a certain viewpoint. In this section, we briefly overview the methods that are most relevant to our work.

Surface reconstruction based One intuitive way to compute the visibility of a point cloud is to reconstruct the surface. Indeed, the projection of the surface as a depth map may be used to estimate which points are not visible. Some methods do not require prior knowledge of the visibility and can therefore be used for visibility estimation. Surface smoothness approaches [START_REF] Lipman | Parameterization-free projection for geometry reconstruction[END_REF][START_REF] Xiong | Robust surface reconstruction via dictionary learning[END_REF] approximate the surface by locally defining operators that weigh surrounding points in order to estimate the local surface. This constrains the reconstructed surface to fit the point cloud as close as possible while ensuring a certain level of smoothness and preserving sharp features. To deal with large amounts of missing data, Volume smoothness techniques [START_REF] Tagliasacchi | VASE: Volume-Aware Surface Evolution for Surface Reconstruction from Incomplete Point Clouds[END_REF][START_REF] Huang | L1-medial skeleton of point cloud[END_REF] exploit the prior of smooth variation of the volume of the reconstructed surface. Unfortunately, these methods are based on strong prior of uniform sampling of the point cloud, which is not suitable for MMS LiDAR point clouds. Primitive based methods [START_REF] Schnabel | Completion and reconstruction with primitive shapes[END_REF][START_REF] Lafarge | Surface reconstruction through point set structuring[END_REF] aim at fitting geometric shapes (i.e. planes, spheres, cylinders, boxes, etc.) in order to reconstruct the scene. However, the complex shapes that can be met in real world scene often jeopardize the results of such methods. Finally, Global regularity approaches (Li et al., 2011a;[START_REF] Li | 2D-3D fusion for layer decomposition of urban facades[END_REF][START_REF] Monszpart | RAPter: rebuilding man-made scenes with regular arrangements of planes[END_REF] take advantage of the repeatability of certain parts of the scene. These methods have shown great strength for the reconstruction of individual regular shapes such as facades or roads but underperform on realistic complete scenes. Although each technique provides satisfying results on specific scenarios, surface reconstruction is a difficult problem, which often requires additional information, such as normals, sufficiently dense input and uniform sampling.

Convex hull based Some methods estimate the visibility based on the local geometry of the 3D point cloud. Based on the raw point cloud (i.e. only 3D positions), [START_REF] Katz | Direct visibility of point sets[END_REF] proposes an approach for estimating which part of the point cloud is not self-occluded given a certain viewpoint. This method admits to perform better on closed shapes. A spherical inversion is performed on the point cloud. The convex hull of the inverted point cloud augmented by the viewpoint position is computed. Then, points that are lying on the convex hull are considered visible, and the rest of the point cloud as non-visible. The acceptance of concave features is tuned by the sphere radius, which is a global parameter so that this method strongly relies on a uniform sampling of the point cloud. Later, this method was improved to handle small changes in the sampling corresponding to noisy acquisitions [START_REF] Mehra | Visibility of noisy point cloud data[END_REF], but still relies on constant density in the point cloud. Moreover, the computational cost of the convex hull [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF] can rapidly increase depending on the wanted concavity. Finally, [START_REF] Katz | Direct visibility of point sets[END_REF][START_REF] Mehra | Visibility of noisy point cloud data[END_REF] are both designed to perform on point clouds that represent closed shapes, acquired from all directions, which is not realistic in urban scenarios where MMS are not able to scan all surfaces.

Likelihood based Different methods aim at estimating the likelihood of a point to be visible given a point of view, by considering its neighborhood. The most common methods rely on the estimation of visibility cones in screen-space [START_REF] Shalom | Cone carving for surface reconstruction[END_REF], and more recently [START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF]. For each point, a visibility cone is estimated, where the apex of the cone is the given viewpoint. The aspect of the cone is directly related to the visibility. Thus, a point that belongs to a wide cone is more likely to be visible than a point that belongs to a narrow cone. However, the threshold on how open a cone should be in order to consider the point visible strongly depends on the point cloud, and can be hard to set.

In this paper, we propose an automatic screenspace method for estimating the visibility of points in a point cloud given a viewpoint. This method makes no assumptions on the sampling or the density of the point cloud and can therefore performed on any point cloud.

Visibility estimation method

The first contribution of this paper is a method for estimating visibility in a 3D point cloud that is robust to high sampling variations.

As illustrated in Fig. 2, points from two objects located at different distances from the given viewpoint overlap in the image plane once projected. In this context, a point is visible only if it lies on the closest object and occluded otherwise. From this observation, we propose an algorithm that considers the neighborhood of a point in screen-space in order to estimate whether this point lies on the closest object or not. The algorithm consists in 4-steps detailed hereafter. be computed using any K-NN alogrithm with a Euclidean distance. The use of the K-NN algorithm defined in [START_REF] Friedman | An algorithm for finding best matches in logarithmic expected time[END_REF] ensures logarithmic computation time while being parallelizable.

Visibility estimation For each point, we want to determine if it lies on the object in its neighborhood that is the closest to the viewpoint. If so, we can consider it as visible. To that end, we compare its position to the closest and the farthest point of its neighborhood. We define the visibility of each point as follows:

α p = e -(dp-d min p ) 2 (d max p -d min p ) 2 (1)
where

d min p = min q∈N (p) d q , d max p = max q∈N (p) d q
The visibility estimation of each point p ∈ P Φ is now given by α p ∈ [0, 1], where α p = 0 means that p is occluded and α p = 1 means that p is surely visible.

Binarization The visibility of a point cloud being a binary notion, we propose the following binarization of α p :

αp = 1 if α p ≥ ᾱ 0 otherwise.
(2)

with ᾱ = 1 Card(P φ ) ∑ p∈P φ
α p the mean of the estimated visibilities. Note that various ᾱ values have been tested such as ᾱ = 0.5 or the median value of the estimated visibilities as discussed Section 5. However, in our experiments on LiDAR data, the mean value remains the best threshold. When point clouds have constant densities, ᾱ = 0.99 appears to be more adequate.

Visibility estimation dataset for LiDAR point clouds

The evaluation of visibility estimation techniques has mostly been done either by visual analysis or by comparison to degradated synthetic models. In [START_REF] Shalom | Cone carving for surface reconstruction[END_REF][START_REF] Katz | Direct visibility of point sets[END_REF], visual results are displayed to show the qualitative performances of each algorithm. In [START_REF] Mehra | Visibility of noisy point cloud data[END_REF], small degradations on synthetic model are applied in order to build groundtruths. Although these methods of evaluation provide convincing results, they do not provide complete and objective quantitative measures on real data. Real data, such as LiDAR, differ from synthetic data in two aspects. The first difference is that the point cloud density is highly variable on real data depending on the distance to the sensor, while constant on synthetic data. The second one is that real urban data only acquire partial representations of each object of the scene as the sensor does not see objects from every possible viewpoints. On the other hand, the synthetic data presented in the related works [START_REF] Shalom | Cone carving for surface reconstruction[END_REF][START_REF] Katz | Direct visibility of point sets[END_REF][START_REF] Mehra | Visibility of noisy point cloud data[END_REF] are always complete 3D objects, which can be seen from any viewpoint. To our knowledge, we give here the first annoted dataset on real urban LiDAR data, which makes the second contribution of this paper.

Overview of the dataset

We propose a manually annotated dataset containing over a 1 million points with the label 1 or 0 depending on if the points are visible or not. This dataset has been obtained by manually labeling 3 point clouds acquired by the RobotCar system [START_REF] Maddern | 1 year, 1000 km: The Oxford RobotCar dataset[END_REF] at different locations, in urban environment. Two of these point clouds are acquired several meters from one another in order to test the stability of visibility estimation methods. The third point cloud corresponds to another location and covers a much wider area which enables testing the limit of methods in case of large distances (> 100m).

Annotations were done manually by comparing the projections of the point clouds to the optical images acquired at the same viewpoints. Fig. 4 presents an overview of the produced dataset. The first row illustrates each scene as acquired from the optical sensor at each viewpoint. The second row shows the projections of each point cloud in the image domain (with the calibration matrices provided by the Robotcar dataset [START_REF] Maddern | 1 year, 1000 km: The Oxford RobotCar dataset[END_REF]), where occluded points are highlighted in red. Finally, the third row shows a 3D visualization of each point cloud, with same color code than above. It illustrates the amount of points to be processed as well as the size of the scenes. The statistics of the dataset are summed up in Table 1. The dataset proposes different levels of visible / occluded points, as well as different size of the scene.

This dataset is publicly available online1 . The archive contains 3 text files in the .xyz format that correspond to each point cloud. In each file, a line corresponds to [x, y, z, u, v, label] where x, y, z are the 3D coordinates of the point, u, v are the 2D coordinates of the point when projected into φ and label is the visibility label (0 for occluded points, 1 for visible points). To ensure good understanding of each of the 3 scenes, we also provide the optical RGB image of size 1280 × 960px associated to each viewpoint.

Experiments & Results

In order to evaluate the performances of our visibility estimation method, we first perform a full numerical and visual comparison between our method and other state-of-the-art methods on (1) the proposed visibility dataset, (2) a point cloud with constant density. Next, we show application of our method to data fusion by performing point cloud colorization from RGB images. All the algorithms are run on Matlab 2018a with a 3.5Ghz CPU.

Evaluation on the Visibility Estimation Dataset

Using our new annotated dataset, we propose an evaluation with two state-of-the art methods, and with our proposed model, against a groundtruth. For each method, we set all the parameters to its optimal value (e.g. the parameters that gives best results against the groundtruth). In our case, we set N = 75 and we detail results for different ᾱ values. We measure the efficiency of each method by computing the following metric:

S(P ) = 1 Card(P ) ∑ P∈P α p × GT p (3)
where GT p corresponds to the annotation of the point P (0 or 1, occluded or visible respectively). This metric aims at capturing the percentage of correctly labeled points provided by each method. The results of this evaluation are displayed in Table 2. Table 2 demonstrates that our algorithm outperforms each compared methods for the 3 scenes. The best scores are obtained by setting the threshold ᾱ equal to the mean of visibility estimations for our method. This observation is explanable. Indeed, as mentionned above, LiDAR aquisitions only capture pieces of the scene. Thus, objects are represented by one of their face only which makes them well separated from one another. In this sense, when an object overlaps an other in screen-space, the mean of the visibility estimations usually represents the visibility of a point that would be in between those two objects. If a point has a visibility estimation above this threshold, it is likely to fall on the closest object, otherwise, it is occluded. Therefore, the mean value can be used when working on LiDAR point clouds because of the way objects are separated from one another, making the method fully automatic in this context.

We also demonstrate that our method operates faster than any other tested method with the ability of treating the whole dataset in less than a second. Moreover, the code is run on a single CPU. Among the 4 steps of the algorithm, the computation of the K-NN is the most time consuming (about 86% of the total running time). Therefore, one can expect much faster running times by operating on GPU with parallel implementation of the K-NN algorithm.

The problem of visibility estimation is a classification problem with two classes: visible and occluded points. Therefore, we enrich our evaluation by com- puting typical classification metrics for each method and display them in Table 3. For each metric, the best scores are obtained using our method. In particular, our method with ᾱ = 0.5 maximizes the true-positives and minimizes the false-negatives. On the opposite, our method with ᾱ as the median of the estimations maximizes the false-positives and true-positives. Once again, using our method with ᾱ as the mean of the estimations provides a good tradeoff between true-positives/true-negatives and false-positives/falsenegatives. On the other hand, [START_REF] Katz | Direct visibility of point sets[END_REF] and [START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF] tend to over-estimate the visibility of each point, resulting in many occluded points being label as visible. This is expressed by the very high percentage of false-positive. We computed accuracy and the F1-score of each method against the ground truth. For both, our method with ᾱ as the mean of the estimation achieves, once again, the best results.

For the task of data-fusion, it is often preferable to discard the maximum of occluded points [START_REF] Bevilacqua | Joint inpainting of depth and reflectance with visibility estimation[END_REF]. Therefore, the number of false-positives has to be kept as low as possible. In this sense, our method provides very satisfactory results, especially using the ᾱ as the mean of estimations when working on LiDAR data.

We conclude this evaluation on LiDAR data by a visual analysis of the results of the different methods. Fig. 5 shows the results of the visibility estimations visualized in 3D. For each result, the dark cone in the bottom left corner represents the viewpoint. Fig. 5(a) shows the point cloud colorized with the depth toward the viewpoint (cold colors for close points, hot colors for far points). Fig. 5(b) shows the annotated groundtruth for this scene, where red points are points that are visible from the viewpoint and dark points are supposed to be occluded. Fig. 5(c) and 5(d) are the results of HPR [START_REF] Katz | Direct visibility of point sets[END_REF] and our method (with ᾱ set as the mean of estimations) respectively. We can see that HPR [START_REF] Katz | Direct visibility of point sets[END_REF] estimates too many points as visible points, especially on the closest points. On the opposite, our method succeeds in discarding occluded points, and provides a result that is very close to the groundruth.

We also illustrate these results as seen from the associated viewpoint in Fig. 6. For better understanding purpose, Fig. 6(a) shows an image acquired from the same viewpoint. In Fig. 6(b), we only display visible points of the groundtruth. Fig. 6(c) and 6(d) shows the results of HPR [START_REF] Katz | Direct visibility of point sets[END_REF] and our method respectively. We can see once again that HPR [START_REF] Katz | Direct visibility of point sets[END_REF] labels too many occluded points as visible and fails to distinguish foreground from background objects. This is mostly due to the fact that this scene presents very high variations of density. In particular, the center of the road concentrates a very high density of points as the sensor is close from the road. Therefore, the convex-hull has to be relaxed enough to fit this region of the point cloud, which leads to visual abberations on regions with lower density. Our method succeeds better results this scene, which demonstrates its robustness against high density variations.

Evaluation on constant density point cloud

In previous section, we demonstrated that our method performs better than other methods for point clouds with high density variations. In this section, we aim at showing that our method remains competitive on constant density point clouds. The Stanford Bunny model is a point cloud (from the Stanford University CG Laboratory) that was created by merging 10 depth aquisitions of a real object and equalizing the density of the fused point cloud. The final point cloud is composed of 31655 points. As each depth acquisition only acquires points that are visible from a single viewpoint, we created a groundtruth by comparing the final point cloud to the points that were aquired at a certain viewpoint. Criterion (3) and classification metrics have been computed for our method and stateof-the-art methods. Results are displayed in Table 4.

Here, the point cloud is of constant density and represents a very smooth object as illustrated Fig. 7. This is a scenario that is perfectly adequate for the HPR algorithm, which shortly outperforms the two other methods. Our method underperforms only by about 1 percent but still remains very efficient on these types of data. Compared to the two other methods, our method fails on tangential points that are located at the boundaries of the projection of the object as it is presented on the last row of Fig. 7. This is mostly due to the fact that on tangential points, the neighborhood covers only a small area, thus the difference between foreground and background is hard to set. These artifacts are limited when using the mean of estimation as visibility threshold, but it increases false-positives. Table 4 also illustrates the classification metrics. We can see that all tested methods reach very good levels of accuracy and F1-score. Our method succeeds better when ᾱ = 0.99 than when using the mean value. Indeed, for complete object, there is no separation between foreground object and background object as was the case for LiDAR point clouds. Thus, only points with high likelihood should be kept to improve results, which justifies ᾱ = 0.99.

Finally, Table 4 assesses that all methods limit the appearance of false-positives while ensuring to gather as many visible points as possible. To that end, HPR and our method succeed the best true-positive/falsepositive ratio, which is ideal for data-fusion purposes, as discussed in next section.

Example of application to data fusion

To conclude our experiments, we show the interest of our visibility estimation for the task of data fusion.

Using the KITTI dataset [START_REF] Geiger | Vision meets robotics: The KITTI dataset[END_REF], we aim at colorizing a 3D LiDAR point cloud acquired in a street using only RGB images. Each point is projected in the image domain of the closest image (e.g. the image that was acquired at the closest position from the point). The point then takes the color of the pixel it projects into only if it is considered visible. Fig. 8 presents the result of the colorization on a point cloud composed of 3289533 points, and is colorized using 40 RGB images. Fig. 8(a) shows the colorization result where all points are considered visible. We can see that artifacts appear as the colors do not match the objects. This is particularly noticeable behind cars where the ground points take the color of the car. Fig. 8(b) displays the colorization result where the visibility is estimated using HPR [START_REF] Katz | Direct visibility of point sets[END_REF]. There, some artifacts appear behind cars as the convex-hull go through the glasses of the car. Moreover, this method discards many visible points on the ground and behind cars compared Ground truth HPR [START_REF] Katz | Direct visibility of point sets[END_REF] Cone [START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF] Ours ( ᾱ = 0.99) Ours (mean) [START_REF] Katz | Direct visibility of point sets[END_REF] for the visibility estimation, (c) colorization with our visibility estimation method. The result provided by our method presents no artifacts on occluded areas, especially behind cars compared to the two other results.

to our method (about 17% less points are colorized). Fig. 8(c) presents the colorization result using our visibility estimation method. We can see that the artifacts behind cars have completely disappeared, while keeping most of the visible points. Finally, due to the number of points, the visibility estimation using HPR [START_REF] Katz | Direct visibility of point sets[END_REF] for each viewpoint takes an average of 10.9 seconds whereas our method processes the point cloud at each viewpoint in about 1.2 seconds.

Conclusion & Perspectives

In this paper, we have proposed a novel method for the visibility estimation in a point cloud. Compared to other methods from literature, this method is very robust to high variations of density. By considering the closest neighbors of each point in screenspace, we defined a criterion in order to automatically determine the visibility of each point. We have also proposed a new annotated dataset for testing the efficiency of point cloud visibility algorithms on real LiDAR urban data. This dataset is composed of over a million manually annotated points. Finally we have compared our method to state-of-the-art methods. We have validated that our method significantly outperforms existing methods on real urban data. Although our method was specifically designed for the estimation of visibility on point cloud with various density (such as LiDAR point clouds), we have also demonstrated that it still remains competitive on point clouds with constant density.

In the future, we would like to focus on the evaluation of our method against photogrammetric point clouds, for tasks such as view generation. We also would like to improve LiDAR / Optic fusion thanks to the visibility estimation.

Figure 1 :

 1 Figure 1: Illustration of the visibility problem. (a) an optical image corresponding to a view point, (b) the 3 main structures of the scene, (c) the projection of the acquired point cloud seen from the same view point with same colors as in (b), (d) the visibility errors brought by the projection (red points should not be visible).

Figure 2 :Figure 3 :

 23 Figure 2: Illustration of notations in 3D and in the screenspace. Blue points corresponds to the points lying on the closest object while red points lies on the farthest object. Although the blue points and red points are well separated in 3D, they overlap in screen-space.

Figure 4 :

 4 Figure 4: Overview of the proposed dataset. First row: optical image corresponding to each viewpoint. Second row: point cloud once projected in the image domain, where red pixels correspond to occluded points. Third row: 3D visualization of each point cloud, with the same color code than in the second row.

Figure 5 :

 5 Figure 5: Results of visibility estimation on the first scene of our visibility estimation dataset.. (a) the point cloud where the heat of the color is proportional to the depth, (b) is the annotated point cloud (red: visible, grey: non-visible), (c) HPR result and (d) our result. The result brought by HPR estimates too many visible points, whereas our method provides a result that is very close to the groundtruth.

Figure 6 :

 6 Figure 6: Results of the visibility estimations on the first scene of the dataset in screen-space. (a) the optical image associated with the point of view, (b) visible points with respect to our annotation, (c) HPR result and (d) ours. Red points in (c) and (d) correspond to misestimated points.

Figure 7 :Figure 8 :

 78 Figure 7: Visual comparison of the visibility estimation from different methods on a point cloud with constant and high density.Each column corresponds to one method. Rows are respectively: the results in 3D, the results in 2D (seen from the viewpoint), and a zoom of the 2D result focused on the ear region. The 3 methods succeed very well in estimating the visibility. Our methods lacks of precision for points that are tangent to the viewpoint as can be seen on the last row.

Table 1 :

 1 Content of the visiblity estimation dataset

		Points	Visibility	Farthest point	Size
	Scene #1	337384	55.5%	75.8m	20.6Mb
	Scene #2	247682	57.0%	54.3m	15.1Mb
	Scene #3	463531	65.9%	179.2m	28.3Mb
	Total	1048597	59.46%	-	64Mb

Table 2 :

 2 Comparison of the scores of two state-of-the art and our visibility estimation methods on our visibility estimation dataset

		HPR (Katz et al., 2007)	Cone (Pintus et al., 2011)	Ours	Ours	Ours
	Threshold	optimal	optimal	ᾱ = 0.5	α p median	α p mean
	POV #1	74.09%	68.76%	90.15%	86.35%	90.96%
	POV #2	69.09%	61.68%	86.95%	86.78%	88.39%
	POV #3	81.55%	75.58%	82.21%	76.35%	83.75%
	Average	74.91%	68.67%	86.43%	83.16%	87.70%
	Total time	7.82s	1.53s	0.91s	1.03s	0.91s

Table 3 :

 3 Comparison of the different methods for point cloud visibility classificationHPR[START_REF] Katz | Direct visibility of point sets[END_REF] Cone[START_REF] Pintus | Real-time rendering of massive unstructured raw point clouds using screen-space operators[END_REF]

	Ours	Ours	Ours

Table 4 :

 4 Comparison of the scores of the different methods on constant density point cloud

		HPR	Cone	Ours	Ours
	Threshold	optimal	optimal	ᾱ = 0.99	α p mean
	Score (Eq. (3))	96.57%	93.75%	95.23%	93.02%
	True-positive	95.17%	88.63%	94.44%	98.07%
	False-positive	1.15%	0.88%	2.14%	6.07%
	False-negative	2.28%	5.37%	2.63%	0.91%
	True-negative	97.82%	98.33%	95.95%	88.50%
	Accuracy	98.25%	96.76%	97.56%	96.40%
	F1-score	98.23%	96.59%	97.54%	96.57%

http://www.labri.fr/perso/pbiasutt/Visibility/

Acknowledgement

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR GOTMI) (ANR-16-CE33-0010-01). This project has also received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 777826.