Jean-Baptiste Charlier
email: j.charlier@brgm.fr

Roger Moussa

Stephanie Pinson

Didier Narbaïs

Jean-François Desprats

Stéphanie Pinson

Recent advances on the characterization of karst-river interactions during floods

Keywords: surface/groundwater interactions, geomorphological index, model, runoff

come

Flood hazard and groundwater resource management in karst catchments require a better understanding of hydrogeological processes and notably groundwater/surface water interactions. This is a great challenge due to the heterogeneity of karst aquifers and due to the diversity of karst types (small perched units, wide sedimentary basins, barren/covered karst, etc.). The aim of this communication is to present recent developments on the characterization of karst-river interactions through two main approaches: a spatial approach aiming at localizing karst areas promoting surface flows, and a temporal approach aiming at modelling lateral flows from karst units in rivers.

The first spatial approach is based on the GIS index IDPR (Index of Development and Persistency of River networks, developed by BRGM©), quantifying the hydrological connectivity to the hydrographic network. From the standard version of the IDPR over France (25 m resolution), we propose to compare IDPR calculations differentiating intermittent and perennial reaches of rivers in order to detect infiltrations zones that contributive temporary to rivers, as many karst units.

The second approach is a modelling framework based on the inverse problem for the diffusive wave model, to simulate lateral flow during floods on a reach between two stations. Knowing the upstream and downstream hydrographs, we can model the lateral one, given information on the hydrological processes involved in the intermediate catchment. Applying such approach on river reaches crossed by karst units is a new way to quantify river losses and gains, characterizing localized recharge and aquifer drainage to rivers, respectively.

We propose to illustrate these two approaches through several case studies in France, where a better characterization of karst-river interactions and flood risk management are critical issues.


Vladislav Popov
email: uladzislau.papou@centralesupelec.fr

Shah Fabrice Boust
email: fabrice.boust@onera.fr

Nawaz Burokur
email: sburokur@parisnanterre.fr

Sondra

Shah Nawaz Burokur

Controlling Diffraction Patterns with Metagratings

come

For long time, the microwave community has approached a particular problem of anomalous reflection by means of reflectarray antennas [1,2]. In such antennas, a linear phase variation is created along the surface, allowing one to reflect incident waves to a desirable angle. With the development of nanofabrication technologies and metasurfaces, the concept of reflectarrays was transposed to infra-red and optical frequency domains [3,[START_REF] Stanislav | Metasurfaces: From microwaves to visible[END_REF]. A metasurface is represented by a 2D dense distribution of subwavelength scatterers and a reflectarray is a particular case of a metasurface which can generally be used for various applications other than anomalous reflections. However, reflectarrays suffer from low efficiencies for angles of anomalous reflection approximately greater than 45 degrees [START_REF] Asadchy | Perfect control of reflection and refraction using spatially dispersive metasurfaces[END_REF]. Extensive research in the area established a strong theoretical ground in the form of equivalence principle [START_REF] Pfeiffer | Metamaterial Huygens' surfaces: Tailoring wave fronts with reflectionless sheets[END_REF] for the design of wavefront manipulation devices based on the use of metasurfaces. As such, multichannel reflection with metasurfaces was demonstrated both theoretically and experimentally in [START_REF] Asadchy | Flat engineered multichannel reflectors[END_REF]. Recently, a metasurface performing highly efficient anomalous reflection at steep angle has been demonstrated in [START_REF] Díaz-Rubio | From the generalized reflection law to the realization of perfect anomalous reflectors[END_REF] on the basis of the concept of metasurfaces possessing strong spatial dispersion [START_REF] Asadchy | Perfect control of reflection and refraction using spatially dispersive metasurfaces[END_REF][START_REF] Epstein | Synthesis of passive lossless metasurfaces using auxiliary fields for re-flectionless beam splitting and perfect reflection[END_REF]. Unfortunately, a theoretical framework to design strongly spatial dispersive metasurfaces has not been developed yet, making the design of a sample time consuming [if it is possible at all] as it requires 3D fullwave simulations. In spite of advances in the field of metasurfaces, drawbacks concerning design complexity and material losses still exist, rendering implementation of high performance devices very challenging in some frequency ranges [START_REF] Ratni | Reconfigurable metamirror for wavefronts control: applications to microwave antennas[END_REF].

In this study we elaborate on the recent concept of metagratings [START_REF] Ra | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF] for the manipulation of reflected waves. Basically, a metagrating is a set of 1D arrays of scatterers such as polarization line currents, separated by a distance of the order of the operating wavelength λ. Polarization line currents are used to cancel scattering in undesirable diffraction orders. Metagratings allow one to significantly decrease the number of constitutive scatterers in contrast to metasurfaces where scatterers are tightly packed in the plane. This reduction can be very attractive to reduce the fabrication complexity as well as the joule losses.

On the theoretical level, metasurface and metagrating are described differently. As a metasurface is composed of deeply subwavelength tightly packed elements, one can introduce averaged surface impedances. Meanwhile, a metagrating is treated as an array of polarization line currents separated by distances much larger than their sizes. Even though, there can be many polarization line currents in a supercell the separation between the currents remains on the order of operating wavelength and one would speculate by introducing average surface impedances.

It has been already shown that having just a single line current per period allows one to cancel specular reflection and perform perfect beam splitting and anomalous reflection [START_REF] Ra | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Epstein | Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis[END_REF][START_REF] Rabinovich | Analytical design of printed-circuit-board (pcb) metagratings for perfect anomalous reflection[END_REF]. In [START_REF] Alex | Perfect anomalous reflection with a bipartite huygens' metasurface[END_REF], the authors numerically and experimentally demonstrated the possibility to perform highly efficient broadband anomalous reflection with a Huygens' metasurface having just two meta-atoms per supercell necessary for cancelling specular reflection. Basically, the same functionality was demonstrated in [START_REF] Ra | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Rabinovich | Analytical design of printed-circuit-board (pcb) metagratings for perfect anomalous reflection[END_REF] where a single meta-atom per supercell and the substrate thickness were used as degrees of freedom instead of two metaatoms per super cell. In this sense, the work in Ref. [START_REF] Alex | Perfect anomalous reflection with a bipartite huygens' metasurface[END_REF] is very similar to the ones on metagratings. Chalabi et al. also demonstrated the possibility to perform nearperfect anomalous reflection using two line currents per super cell [START_REF] Chalabi | Efficient anomalous reflection through near-field interactions in metasurfaces[END_REF] that are necessary for eliminating reflection in the zeroth and minus first diffraction orders. Recently, an implementation of a graphene-based tunable metagrating operating in the THz frequency range was suggested in [START_REF] Ra | Reconfigurable metagratings[END_REF].

In the present work, we study a general case of metagratings having N polarization line currents per super cell. This generalization is a necessary step for controlling diffraction patterns when the number of plane waves scattered in the far-field is greater than three. Although, the authors of Ref. [START_REF] Alex | Perfect anomalous reflection with a bipartite huygens' metasurface[END_REF] discussed the number of metaatoms per super cell necessary for controlling arbitrary number of plane waves diffracted in the far-field, a clear theory for designing a N-meta-atoms Huygens' metasurface was not elaborated.

Gaining control over many diffraction orders can be particularly interesting for implementing tunable devices and performing multichannel reflection. Indeed, having many identical wires but being able to control polarization currents in each of them allows one to perform all possible transformations of the diffraction pattern with the same device, where the only restriction remains the device size. Moreover, metagratings can operate in a broad frequency range as usually does not require resonance response of meta-atoms. Broadband response of metagratings with a single and couple of polarization currents per supercell was demonstrated in [START_REF] Rabinovich | Analytical design of printed-circuit-board (pcb) metagratings for perfect anomalous reflection[END_REF] and [START_REF] Chalabi | Efficient anomalous reflection through near-field interactions in metasurfaces[END_REF], respectively.

As a physical system, we consider a 1D periodic array of polarization electric line currents placed over a grounded dielectric substrate of thickness h and excited by an incident harmonic TE-polarized plane wave at angle θ where exp[jωt] time dependence is assumed. The array has period L and consists of super cells each having N equally separated line currents by the distance d = L/N . The schematics of the considered system is presented in Fig. 1 (a). A line current is imagined as a tightly packed row of point dipoles orientated in the same direction, see Fig. 1 (b). Practically, one can realize the dipoles as the loaded rods considered in Fig. 1 (c) and(d).

In the presence of the grounded substrate the excitation field takes the following form

E exc x (y, z ≤ -h) = e -jβ0z + R T E 0 e jβ0(z+2h) e -jk sin[θ]y . (1)
Electric line currents in the array are represented as current densities J nq (r) = I q exp[-jk sin[θ]nL]δ(y -y nq , z + h)x 0 where δ(y, z) is the Dirac delta function, y nq = nL+ (q -1)d, n and q take integer values from -∞ to +∞ and from 1 to N , respectively. The term exp[-jk sin[θ]nL] represents the phase variation of the currents introduced by the incident wave. Radiation of the array of electric line currents is represented by a series of Hankel functions [START_REF] Leopold | Radiation and scattering of waves[END_REF][START_REF] Tretyakov | Analytical modeling in applied electromagnetics[END_REF] by the array outside the substrate can be written as

E x (y, z < -h) = - kη 2L +∞ m=-∞ ρ (I) m (1 + R T E m ) β m e -jξmy+jβm(z+h) , (2)
E y = E z = 0. Corresponding magnetic fields can be found from the Maxwell equations. The series represent superpositions of plane waves having tangential component of wave vector equal to ξ m = k sin(θ) + 2πm/L, the longitudinal component is given by

β m = k 2 -ξ 2 m outside the substrate [k = ω √ εµ and k s = ω √ ε s µ s are
respectively the wave numbers outside and inside the substrate]. Thus, R T E m is Fresnel's reflection coefficient from the grounded substrate of a plane wave having tangential component of the wave vector equal to ξ m . Each current contributes to the amplitudes of the plane waves via the introduced quantity ρ

(I) m ρ (I) m = N q=1 I q exp[jξ m (q -1)d]. (3)
One can recognize in Eq. ( 3) a discrete Fourier transformation.

In general case when a plane wave illuminates a metagrating one can find r + l + 1 scattered plane waves in the far-field, where r and l are largest integers satisfying the conditions β r > 0 and β -l > 0. However, we can arbitrary control all of the r + l + 1 plane waves if the number N of line currents in a super cell is equal to r + l + 1. Indeed, amplitude A T E m of the m th plane wave depends on ρ (I) m which is determined by the currents I q [see Eq. ( 2)]

A T E m = - kη 2L ρ (I) m (1 + R T E m )e jβmh β m + δ m0 R T E 0 e 2jβ0h , (4)
where δ m0 is the Kronecker delta accounting for the incident wave reflected from the substrate. By setting all the amplitudes A T E m (m ∈ [-l, r]), one can find necessary I q from the corresponding ρ (I) m which are related via Eq. ( 3). Thus, by designing currents I q one can perform all possible transformations of the diffraction pattern, e.g. beam splitting, anomalous reflection, multichannel reflection, etc. When implementing line currents as thin perfectly conducting wires one can obtain necessary currents I q by loading wires with suitable impedance densities Z q . The last are found from

Z q I q = E (exc) x (y 0q , -h) -Z in I q - N p=1 Z (m) qp I p (5)
where the right-hand side simply represents total elec-tric field at the location of the qth wire in the zeroth supercell [y 0q = (q -1)d, z = -h]. Here, we also introduce notations for the input impedance density of wires Z in = kηH

(2) 0 [kr 0 ]/4, with H

(2) 0 [kr 0 ] being the Hankel function of the second kind and r 0 the radius of the wires, and for the mutual impedance densities Z (m) qp which account for interaction between the wires and between the wires and the grounded substrate. When a wire is realized as a perfectly conducting strip of width w, the radius is r 0 = w/4 [START_REF] Tretyakov | Analytical modeling in applied electromagnetics[END_REF]. See Supplementary Material at [19] for the derivation of Eq. ( 5) and expressions of the mutual impedance densities.

Generally, currents found from (4) correspond to active and lossy loads Z q calculated from (5). From a practical point of view, we are interested only in passive and lossless metagratings where [Z q ] = 0, i.e. which cannot radiate energy by themselves and do not require engineered joule losses. A metagrating should redistribute the energy of the incident wave between r + l + 1 diffracted in the far-field plane waves. Then, the power conservation condition when assuming a unity amplitude of the incident wave reads as

r m=-l α m = 1, α m = A T E m 2 β m β 0 , (6)
where α m is the part of the incident energy going in the m th diffraction order.

In contrast to the case of metagratings having a single line current in a super cell [START_REF] Ra | Metagratings: Beyond the limits of graded metasurfaces for wave front control[END_REF][START_REF] Epstein | Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis[END_REF][START_REF] Rabinovich | Analytical design of printed-circuit-board (pcb) metagratings for perfect anomalous reflection[END_REF], when it comes to greater number of line currents per supercell there are no exact analytical formulas for reactive load impedance densities necessary for obtaining some diffraction pattern. To approach this problem we develop a very simple real valued genetic algorithm [START_REF] Wright | Genetic algorithms for real parameter optimization[END_REF] which allows one to find reactive Z q with given impedance reactivity accuracy p for a desired diffraction pattern obtained with given transformation accuracy α. The impedance reactivity accuracy is defined in accordance with the following inequality N q=1 |Re[Z q ]/Z q | 2 < p. A diffraction pattern is set by assigning to all α 0 m certain values. Phases

φ m = arg[A T E
m ] are assumed to be not important and assigned randomly. Transformation accuracy α means that one is satisfied with a transformation when the part of the incident energy going in the m th diffraction order is within the range α m = α 0 m ± α. Still, at each step the genetic algorithm deals with α m > 0 constrained by the energy conservation condition [START_REF] Pfeiffer | Metamaterial Huygens' surfaces: Tailoring wave fronts with reflectionless sheets[END_REF].

When designing a metagrating, one should also take care of choosing parameters of the substrate. First of all, when substrate's thickness is varied the value of the excitation field Eq. (1) on a metagrating passes through zeros as illustrated in Fig. 2 (a). Clearly, a metagrating cannot be excited when the excitation field is zero on its plane. And secondly, the reflection coefficient R T E m as a function of h has poles when m is such that k < ξ m but k s > ξ m [Fig. 2 (b)]. The poles correspond to excitation of waveguide modes inside the substrate. Thus, assuming ε s and µ s of the substrate are set, one should choose the thickness: (i) corresponding to vicinity of the maximum of the excitation field on a metagrating and (ii)

|R T E m (h)| = ∞.
One can realize the line currents as dense 1D arrays of loaded dipoles [separated by distance B λ and having lumped load equal to Z] as in Figs. 1 (b) and(c). Then, the load impedance density is simply Z/B. A capacitive load can be realized as a printed circuit capacitance as illustrated in Fig. 1 (d) for which Z = -jηκ/(Aε ef f ) [when other parameters B and w are fixed], κ is the proportionality factor, ε ef f is approximated as (1+ε s )/2 and µ s is assumed equal to 1. The proportionality factor κ was introduced in [START_REF] Epstein | Unveiling the properties of metagratings via a detailed analytical model for synthesis and analysis[END_REF] for the case of a single line current per supercell. However, it turns out that in the general case of many line currents per supercell one can successfully use the same proportionality factor for all currents, i.e. independently on q. Thus, when load impedance densities Z q are found from the genetic algorithm, one can easily calculate arms lengths of necessary printed capacitors as A q = -κ/( [Z q ] λ η ε ef f ), κ = λκ/B. In order to validate the developed theoretical basis, we perform 3D full-wave simulations with COMSOL Multiphysics. We demonstrate three examples of metagratings designed to operate at 10 GHz [λ ≈ 30 mm] and perform different transformations of the diffraction pattern as shown in Fig. 3. A polarization line current is implemented as a 1D array of capacitively loaded perfectly conducting strips [as it schematically shown in the top row of Fig. 3]. In all the examples normally incident plane wave is assumed, i.e. θ = 0.

When performing a large-angle anomalous reflection with a metagrating, overall there are three diffraction orders and therefore only three polarization line currents per supercell are necessary to cancel the -1 st and 0 th diffraction orders, as illustrated in Fig. 3 (a). Figure 3 (d) depicts the frequency response of the metagrating performing anomalous reflection at angle of 50 o . The situation is more difficult in case of a small-angle anomalous reflection with the presence of many high diffraction orders and when the energy should be scattered only in the first one. Indeed, in the example of Fig. 3 (b) there are nine diffraction orders and the metagrating with nine polarization currents is used to cancel scattering in all of them except the first one corresponding to an anomalously reflected wave at angle of 12.5 o . Figure 3 (e) demonstrates the frequency dependence of the metagrating's performance efficiency. Clearly, metagratings are not restricted to anomalous reflection application and can be used for multichannel reflection. One can distribute the energy of an incident wave between all diffraction orders in a desirable manner. For instance, Figs. 3 (c) and (f) demonstrate the scenario when the metagrating having five polarization currents is used to split normally incident waves between the -2 nd and 1 st diffraction orders and cancel scattering in the other three diffraction orders.

In conclusion, it has been shown that a metagrating having the number of polarization line currents per super cell equal to the number of plane waves scattered in the far-field can be used for controlling the diffraction pattern. Namely, equations ( 3) and (4) allowing one to find currents realizing desirable transformations have been derived. Since there are no analytical formulas of reactive load impedance densities (5) for direct design, genetic algorithms have been implemented for that purpose. The diffraction orders control has been demonstrated by means of 3D full-wave simulations on the ex-amples of anomalous reflection and equal redistribution of the energy of the incident wave between two diffraction orders.

The validation results can be very interesting for the metamaterials community to perform highly efficient control of light scattering. It allows one to significantly decrease the number of used elements and simplify the design, which is very convenient for optical and infra-red frequency ranges. Our findings also may serve as a way for development of efficient tunable antennas in the microwave domain.

E (loc) x

(y 0q , -h), where E (loc) x (y 0q , -h) is the total electric field at the location of the q th line current [y 0q = (q -1)d, z = -h]. The self-action of the q th current in the zeroth supercell is given by -kηI q H

(2) 0 [kr 0 ]/4 [r 0 is the effective radius of the line current which depends on its concrete practical realization]. For instance, when a line current is imitated as a perfectly conducting strip of width w the radius r 0 = w/4 [3]. The electric field created by q th line current from all supercells [but the zeroth one] reads Summarizing three previous equations and keeping in mind the excitation field the load-impedance densities Z q can be calculated by the following formula

Z q = 1 I q   E (exc)
x (y 0q , -h) -kη 4 The cumbersome equation (B2) can be rewritten in the compact form by introducing mutual impedance densities

Z q I q = E (exc)
x (y 0q , -h) -Z in I q -N p=1 Z (m) qp I p , (B2)

where Z in = kηH

(2) 0 [kr 0 ]/4 is the input impedance density of wires and the mutual impedance densities Z (m) qp which account for interaction between the wires and between the wires and the grounded substrate are expressed as

Z (m) qp = kη 4 ∞ n=-∞ H (2

FIG. 1 .

 1 FIG. 1. (a) System under consideration: a periodic array of line currents Jnq = Iq exp[-jk sin[θ]nL]δ(y -ynq, z + h)x0 (blue circles) placed on PEC-backed dielectric substrate having permittivity εs, permeability µs and thickness h. The array is excited by a plane wave incident at angle θ and having TE polarization. (b) A line current implemented as a 1D array of loaded dipoles. (c) A PEC strip dipole of length B and width w loaded with lumped impedance Z. (d) A PEC strip dipole loaded with printed circuit capacitance having arms of length A.

FIG. 2 .

 2 FIG. 2. (a) Dependence of the excitation field [θ = 0] acting on a metagrating on the thickness of the substrate h when εs = 4.5 and µs = 1. (b) Absolute values of the R T E 2 vs. the thickness of the substrate when θ = 0, r = l = 1, N = 3 and L = λ/ sin[60 o ]. The rest of R T E m does not have poles under these parameters. λ is the operating vacuum wavelength.

FIG. 3 .

 3 FIG. 3. The top row of figures demonstrates schematics of simulated metagratings with (a) N = 3, (b) N = 9 and (c) N = 5, the green and white lobes depict excited and canceled diffraction orders respectively. Figures in the bottom row depict obtained from 3D full-wave simulations frequency responses of the metagratings corresponding to the figures in the top row. (a), (d) The example of anomalous reflection at angle of 50 o with the metagrating having N = 3, L = λ/ sin[50 o ], and loads [in η/λ, κ = 26.80 mm] Z1 = -14.5j, Z2 = -6.86j, Z3 = -4.43j. (b), (e) The example of small-angle anomalous reflection [of 12.5 o ] with metasurface having N = 9, L = 4λ/ sin[60 o ] and loads [in η/λ, κ = 27.43 mm] Z1 = -j7.62, Z2 = -j6.96, Z3 = -j6.19, Z4 = -j5.55, Z5 = -j5.18, Z6 = -j3.57, Z7 = -j3.02, Z8 = -j18.7, Z9 = -j10.1. (c), (f) The example when out of five only -2 nd and 1 st diffraction orders are equally excited with the metagrating having N = 5, L = 2λ/ sin[50 o ] and loads [in η/λ, κ = 27.18 mm] Z1 = -9.00j, Z2 = -5.88j, Z3 = -6.59j, Z4 = -3.03j, Z5 = -5.14j. The substrate is Arlon AD450 [εs = 4.5], h = 3 mm, B = λ/10 = 3 mm and w = 3mil ≈ 76.2µm.

II

  p e -jk sin[θ]nL H (2) 0 [k|(q -p)d -nL|].Finally, the influence of the grounded substrate is acp e jξm(p-q)d R T E m β m .

  ) 0 [k|(q -p)d -nL|]e -jk sin[θ]nL

Supplementary Material: Controlling Diffraction Patterns With Metagratings

A single line current J(r) = Iδ(y, z)x 0 radiates as [1]

is the second kind Hankel function. Corresponding magnetic field is found from the Maxwell equations

Consequently, radiation field of a periodic array of N electric line currents per super cell J nq (r) = I q exp[-jk sin[θ]nL]δ(y -y nq , z)x 0 [where y nq = nL + (q -1)d, n and q take integer values from -∞ to +∞ and from 1 to N respectively] is given by a series of Hankel functions

that can be calculated by means of the Poisson's formula stating that for a given "good" function f (w) the following identity takes place

Thus,

e -jξmy-jβm|z| β m .

(A6)

The effect of the grounded substrate on the field radiated by the array J nq (r) = I q exp[-jk sin[θ]nL]δ(yy nq , z + h)x 0 can be derived in the same manner as in [2]. After some algebra one would arrive at the following expressions for the electric field profile inside (0 > z > -h) and outside (z < -h) the substrate Appendix B: Calculation of load-impedance densities Necessary currents I q can be achieved by engineering the load-impedance densities Z q defined as I q Z q =


Fernando I Ireta Muñoz

Andrew I Comport
email: andrew.comport@cnrs.fr

Point-to-hyperplane ICP: Fusing Different Metric Measurements for Pose Estimation

Keywords: Point-to-hyperplane, Visual Odometry, RGB-D Pose estimation, Visual SLAM

The objective of this article is to provide a generalized framework of a novel method that investigates the problem of combining and fusing different types of measurements for pose estimation. The proposed method allows to jointly minimize the different metric errors as a single measurement vector in n-dimensions without requiring a scaling factor to tune their importance. This paper is an extended version of previous works that introduced the Point-to-hyperplane ICP approach. In this approach an increased convergence domain and a faster alignment was demonstrated by considering a 4-dimensional measurement vector (3D Euclidean points + Intensity). The method has the advantages of the classic Point-to-plane ICP method, but extends this to higher dimensions. For demonstration purposes, this paper will focus on a RGB-D sensor that provides color and depth measurements simultaneously and an optimal error in higher dimensions will be minimized from this. Results on both, simulated and real environments will be provided and the performance of the proposed method will be carried on real-time visual SLAM.

Introduction

One of the most common problems in view registration is estimating the pose that relates sets of measurements obtained by a moving sensor (or sensors). This problem has been widely studied by the computer vision and robotics communities and it is specially applied for 3D reconstruction, visual odometry and autonomous navigation tasks.

Depending on the type of sensor, different types of measurements of the environment can be registered through pose estimation. Classically, when more than one type of sensor is employed for pose estimation, the alignment between the extended measurements has been achieved by minimizing each sensor's error separately or by jointly optimizing over each type of measurement in a so-called hybrid -manner.

Nowadays, the availability of RGB-D sensors such as the Microsoft Kinect V1, V2 or Asus Xtion have provided the possibility to acquire color and depth information simultaneously at a considerably high framerate, which has been useful for real-time pose estimation. The metric information obtained by RGB-D sensors has been individually studied in the literature. One case is by using depth images, where geometric-based methods, such as the well known Iterative Closest Point (ICP) [START_REF] Besl | A method for registration of 3-d shapes[END_REF] and its variants, have demonstrated the ability to obtain robust alignments when enough geometric information is available and they can obtain fast alignment if the datasets are closely overlapping. Particular variants such as the Point-to-plane ICP strategy [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] and the Generalized-ICP [START_REF] Segal | Generalized-ICP[END_REF] have demonstrated to be the most effective and robust methods when combined with robust estimations approaches such as the M-estimators [START_REF] Huber | Robust statistics[END_REF]. On the other hand, color images have been used for pose estimation processes by performing photometric-based minimization. Strategies as direct approaches based on view synthesis [START_REF] Irani | About direct methods[END_REF] or feature-based strategies such as: SIFT (Scale-Invariant Feature Transform) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], SURF (Speeded Up Robust Features) [START_REF] Bay | Speeded-up robust features (surf)[END_REF], BRIEF (Binary Robust Independent Elementary Features) [START_REF] Calonder | Brief: Binary robust independent elementary features[END_REF] or ORB (Oriented FAST and Rotated BRIEF) [START_REF] Rublee | Orb: An efficient alternative to sift or surf[END_REF] have been widely used. In [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], the feature descriptors are obtaining by computing a histogram of local oriented gradients around a keypoint, however, it can be computationally expensive due to the high dimensionality of its descriptors. An improved version has been proposed in [START_REF] Bay | Speeded-up robust features (surf)[END_REF], which relies on local gradient histograms but the matching is accelerated by using integral images. However, methods [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Bay | Speeded-up robust features (surf)[END_REF] are highly discriminant and features of the same object under different illumination cannot be properly detected. A similar performance is shown in [START_REF] Calonder | Brief: Binary robust independent elementary features[END_REF], but the feature matches are improved by training a set of classification trees and by using the Hamming distance as the metric for matching instead of Euclidean distance. Finally, an efficient alternative has been presented in [START_REF] Rublee | Orb: An efficient alternative to sift or surf[END_REF], where an efficient computation of BRIEF features is performed.

It can be noted at this point that feature-based approaches first extract geometric information from the image before performing estimation on a geometric error. Therefore, for purposes of this article, feature-based approaches won't be detailed here since they can be considered as a sub-part of direct approaches.

The geometric-based and photometric-based approaches share much similarity and subsequently, the common pose estimation framework of both strategies involves the following non-linear Iteratively Re-weighted Least Squares (IRLS) process:

(1) Acquisition of the measurements at different times.

(2) Transform (warp) the measurements using the last pose estimate.

(3) Find the closest points between the sets of measurements. (4) Minimize the robust weighted error functions. [START_REF] Irani | About direct methods[END_REF] Estimate a new incremental update on the pose. ( 6) Perform all the steps from 2 until convergence.

Recently, several strategies have combined geometric and photometric-based methods together to obtain the main benefits of each via a so-called hybrid method (The main recent surveys are cited in [START_REF] Fang | Experimental evaluation of RGB-D visual odometry methods[END_REF]). The advantages of hybrid approaches in combining different measurements include increased efficiency, accuracy and robustness for pose estimation processes. However, the contribution of each measurement during the minimization process should be weighted by a tuning parameter λ, which scales the relative importance of each measurement (Figure 1). Various prominent hybrid methods proposed in the literature are those that simultaneously minimize the geometric and photometric error functions in real-time such as [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF][START_REF] Meilland | On unifying key-frame and voxel-based dense visual SLAM at large scales[END_REF][START_REF] Whelan | Elasticfusion: Dense slam without a pose graph[END_REF][START_REF] Morency | Stereo tracking using ICP and normal flow constraint[END_REF][START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF]. The aforementioned methods differ in how the tuning parameter is estimated and how the closest points are found.

The cited hybrid strategies in this paper do not necessarily consider the color and depth simultaneously when computing the closest points. All the methods perform the closest point searching separately for both color and depth except for [START_REF] Morency | Stereo tracking using ICP and normal flow constraint[END_REF], which estimates the closest points using a k d-tree (k -dimensional) in a 4-dimensional space (3D Euclidean points + intensity). Finding the closest points by considering the fused information increases the accuracy of finding the true nearest neighbours, however, this approach requires an efficient search in a higher dimensional space. Depending on how the closest points are found, this step can potentially be the most computationally expensive part of the pose estimation pipeline.

The choice of λ has a huge influence while estimating the pose. If the parameter is well determined, then it can speed-up the alignment and increase the convergence rate. Depending on how λ is chosen, a variety of hybrid strategies have been categorized into adaptive or non-adaptive methods in [START_REF] Muñoz | A proof that fusing measurements using point-to-hyperplane registration is invariant to relative scale[END_REF]. Basically, adaptive methods are those that determine the tuning parameter at each iteration of the minimization process and the non-adaptive methods estimate λ only once and its value is used for all the following iterations. For purposes of this article, methods that perform real-time tasks for pose estimation are selected for comparison.

Adaptive methods such as [START_REF] Meilland | On unifying key-frame and voxel-based dense visual SLAM at large scales[END_REF] and [START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF] estimate the tuning parameter by obtaining the ratio between the Median Absolute Deviation (MAD) of the photometric and geometric error functions. In [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF] the uncertainty between the metric measurements is compensated by computing the covariance matrix between both metric errors and in [START_REF] Morency | Stereo tracking using ICP and normal flow constraint[END_REF] the tuning parameter is obtained by a sigmoid function which increases the importance of the photometric error over the geometric error or viceversa. The strategy proposed in [START_REF] Whelan | Elasticfusion: Dense slam without a pose graph[END_REF] is classified as the non-adaptive since λ is chosen experimentally.

The aim of this article is to provide an extended framework of a previous work on fusing different metric measurements via the Point-to-hyperplane ICP approach [START_REF] Muñoz | Point-to-hyperplane RGB-D pose estimation: Fusing photometric and geometric measurements[END_REF]. The invariance to any tuning parameter will be proven mathematically, which will demonstrate that the Point-to-hyperplane ICP method is invariant to λ in hybrid pose estimation processes if the normals are estimated in higher dimensions. Particularly, here the method is applied for RGB-D pose estimation in a 4D and 6D space by fusing both geometric and photometric techniques based on Point-to-plane ICP and direct methods, respectively. With respect to previous work, this paper addresses the issue of computing the 4D normal when geometric or photometric information together are not available. Various real RGB-D sequences that allow to better compare texture vs structure will be compared for the proposed method alongside hybrid strategies that estimate a scale factor. This paper is organized as follows: In Section 2 an overview of the hybrid method briefly explains how the RGB-D pose estimation can be performed by jointly minimizing over the color and depth by using a direct method for the color and Point-to-plane ICP for the geometric 3D points. In this section it will be shown that the tuning parameter λ has a huge influence on these methods. In Section 3 the Point-to-hyperplane ICP method will be introduced and the invariance to the scale parameter will be demonstrated by minimizing the error as a single vector for n-dimensions. Finally, extended results for both, real and synthetic environments, will be shown.

Hybrid RGB-D pose estimation

Hybrid approaches have been useful when color or depth alone are not significant enough for obtaining a correct alignment between RGB-D frames, that have been acquired at different times and that are not in correspondence. A IRLS pose estimation process can be employed to minimize the geometric or photometric error functions separately, but hybrid methods estimate the unknown pose by iteratively minimizing the non-linear error functions simultaneously. Hybrid methods can converge faster than using geometric or photometric approaches individually, and attempt to retain the main benefits of each by weighting their respective contribution.

Joint error for pose estimation

The generated errors between two sets of extended measurements (color + depth) can be jointly minimized to estimate the pose since the color and depth pose estimation pipeline shares too much similarity. So-called hybrid methods have been introduced to minimize both error functions simultaneously, where a 3D Euclidean point P i ∈ R 3 is associated with an unique intensity I i by weighting each contribution with an uncertainty factor λ λ λ. Consider here two augmented point clouds obtained at different times. Let M * be the reference point cloud and M be the current point cloud measurements. The hybrid error function for the i -th joint measurement vector, e Gi and e Ii (geometric and photometric error, respectively) can be represented as:

e Hi = λ λ λ
e Gi e Ii = λ λ λ P * i -w(P i , x)

I I I * i -w(I I I i , x) ∈ R 6 (1)
where λ λ λ = diag(λ X , λ Y , λ Z , λ I ) = λ λ λ G 0 0 λ I are the tuning parameters which weigh the contribution of each error function. For 3D Euclidean points, λ λ λ G = diag(λ X , λ Y , λ Z ). The alignment between the measurement vector M * = [P * I * ] and M = [P I] is found by iteratively minimizing the error function [START_REF] Besl | A method for registration of 3-d shapes[END_REF]. This involves transforming the current dataset M with the estimated pose x with the transformation function represented here as: w(•). A 3D Euclidean point can be determined by using the depth-back-projection function as:

P i = K -1 p i Z i = X i Y i Z i ∈ R 3 ,
where Z i ∈ R + is the depth measurement for each pixel coordinate p i = [u i v i 1] ∈ R 3 of the depth image and K ∈ R 3×3 is the intrinsic camera calibration matrix as:

K =   f w f s θ c x 0 f h c y 0 0 1   R 3×3 (2)
where f is the focal distance, s θ is the skew angle of a pixel (which is usually set to 0), w and h is the width and height of the image, respectively and c x , c y are the coordinates of the center of the image.

The 6DOF (Degrees of freedom) pose parameter x can be decomposed into rotational and translational components and it will be defined here via the homogeneous transformation matrix

T(x) = R(x) t(x) 0 3 1 ∈ SE(3)
which is parametrized by the linear υ υ υ ∈ R 3 and angular velocity ω ω ω ∈ R 3 , respectively. The relationship between the velocity twist and the homogeneous pose matrix is given by the exponential map as T(x) = e [x] , with the operator [•] defined as:

[x] = [ω ω ω] × υ υ υ 0 0 ∈ SE(3) where [•] × is the skew symmetric operator.
Here the non-linear error defined in ( 1) is minimized iteratively using a Gauss-Newton approach to compute the unknown 6DOF pose parameters with increments given by:

x = -(J WJ) -1 J W λ λ λ G e G λ I e I ∈ R 6 (3)
where J = [ J G J I ] is the stacked Jacobian matrix obtained by deriving the stacked error functions, and the weight matrix

W = diag(ρ 1 , ρ 2 , • • • , ρ n )
contains the stacked weights associated with each set of coordinates obtained by M-estimation [START_REF] Huber | Robust statistics[END_REF]. Often, robust M-estimation is performed separately for each different measurement type. The cited hybrid strategies [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF][START_REF] Meilland | On unifying key-frame and voxel-based dense visual SLAM at large scales[END_REF][START_REF] Whelan | Elasticfusion: Dense slam without a pose graph[END_REF][START_REF] Morency | Stereo tracking using ICP and normal flow constraint[END_REF][START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF] perform the Point-to-Plane ICP algorithm [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] and a direct image-based method [START_REF] Irani | About direct methods[END_REF] whilst minimizing the error simultaneously. Generally, these approaches minimize an error function1 similar to:

e Hi = λ G N * i (P m i -P w i ) λ I (I I I m i -I I I w i ) ∈ R 6 (4)
where P w i = Π 3 TT (x) P * i ∈ R 3 is the warped 3D point and N * i ∈ R 3 is the surface normal for each 3D point P * i ∈ R 4 . In the photometric term, I I I w i = I I I ω( TT(x), P * i ) is the warped image through the geometric warping function ω(•) as:

p i w = KΠ 3 TT(x)M i * e 3 TT(x)M i * =   u w i v w i 1   =   c x + f x X w i /Z w i c y + f y X w i /Z w i 1   ∈ R 3 (5)
where The closest image intensity is found by interpolating the current intensity function at the warped pixel coordinates. Therefore, the corresponding intensities can be estimated as:

Π 3 = [1, 0] ∈ R 3×4 projects
I w i (p * i ) = I i (p w i ) ∈ Z + .
The 3D point correspondences and the matched intensities are defined as P m i and I m i , respectively. Finally, the pose estimation T(x) is computed at each iteration and is updated incrementally as: T ← TT(x) until convergence. Advanced˙Robotics˙Journal˙2017 13)) when (a) λ is not estimated, (b) when the intensities and 3D points are normalized between [0, 1] [18] and (c) when λ is estimated at each iteration [START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF]. An example of a cost function is shown in (d), which indicates that the number of iterations can be improved when a good choice of λ is made.

Uncertainty between depth and color

Since the geometric and photometric measurements and subsequently their uncertainties are not in the same units or order of magnitude, the contribution in the minimization process of each measurement should be weighted to compensate for the relative uncertainty between the different error functions. Consider as an example the case when the intensity is almost uniform in the scene but the geometric features are not (Figure 2(b) and 2(a)). The pose estimation function should give more importance to the geometric features since the errors generated in the photometric term are not significant enough to constrain all degrees of freedom. On the other hand, the opposite case can be found when rich texture can be registered from flat surfaces (Figure 2(d) and 2(c), the geometric information does not constrain all degrees of freedom for obtaining the alignment. An example of the influence of λ λ λ = diag(λ X , λ Y , λ Z , λ I ) in the minimization process is shown in Figure 3 where a Gaussian distribution has been fitted into the residuals.

As is shown in [START_REF] Besl | A method for registration of 3-d shapes[END_REF], each intensity is associated with its corresponding 3D Euclidean point through a matrix λ λ λ that scales the importance of the geometric points w.r.t. the intensities. As mentioned in the introduction, many methods have been proposed to choose this parameter ranging from manual tuning to more complex estimation approaches. Manually fixing λ λ λ is not optimal nor efficient for real-time applications, and estimating its adequate value can require extra computational cost. Various strategies, which obtain λ λ λ in different ways, have been cited in [START_REF] Muñoz | Point-to-hyperplane RGB-D pose estimation: Fusing photometric and geometric measurements[END_REF]. Three efficient real-time possibilities will be considered here and they will be compared in the results section. They include adaptive methods: such as the ratio of the Median Absolute Deviations (MAD) [START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF] or computing the covariance matrix for each measurement vector as in [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF], and non-adaptive methods: using the normalization of the metric measurements to scale the relative error distributions. In the first iteration, the closest points can be estimated by finding in a k d-tree with the 4D vector and linear interpolation in the following iterations. The method is invariant to a tuning parameter λ.

Point-to-hyperplane ICP

The aim of this paper is to present a method that performs RGB-D registration by minimizing over color and depth components simultaneously in such a way that it is invariant to the scale between both measurements. The n-D space generated by considering measurement vector of n-Dimensions, has additional degrees of freedom. The proposed method therefore consist in extending the classic Point-to-plane method [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] for 2D and 3D points, to higher dimensions. Therefore, the estimated n-dimensional normal (as well as the 3D normal) will be orthogonal to a surface in n-dimensions which spans both geometry and color. This n-dimensional surface will be referred in this paper as the hyperplane.

Based on the Point-to-plane method for 3D points, an error function in higher dimensions can be defined as follows:

e Hi = λN * i (M * i -w(M i , x)) (6)
where a tuning parameter λ = det(λ λ λ) is added to deal with the uncertainty between different measurements and the normal N * i is perpendicular to the formed hyperplane. For the purpose of this article, 4 dimensions will be considered (3D Euclidean points + intensity) for the experimentation (See Figure 4). The 4-vector is defined as

M i = [P i I i ] ∈ R 4 .
The normals N * are computed on the reference 4D measurements vector M * , which will be referred throughout this paper as the reference dataset. Therefore, the pose vector x can be estimated by iteratively minimizing the error function that projects the Point-to-point distance onto the normal direction as:

x = -(J WJ) -1 J We H (7)

Invariance to a tuning parameter λ

The invariance to any scale factor λ λ λ = diag(λ 1 , λ 2 , • • • , λ n ) in the Point-to-hyperplane ICP method was experimentally observed in [START_REF] Muñoz | Point-to-hyperplane RGB-D pose estimation: Fusing photometric and geometric measurements[END_REF] and a demonstration to this invariance is given in [START_REF] Muñoz | A proof that fusing measurements using point-to-hyperplane registration is invariant to relative scale[END_REF], where the mathematical proof of the following lemma is given:

Lemma. The integrated error e H in n-dimensions is invariant to the relative scale λ if it is minimized by a Point-to-hyperplane ICP method.

e Hi = N * i (M * i -w(M i , x)) = λN * i (M * i -w(M i , x)) (8)
The projection of the error onto the normal direction has the effect of canceling out the effect of λ = det(λ λ λ) between the geometric and photometric error since the direction of the normal is invariant to any tuning parameter. The mathematical proof given in [START_REF] Muñoz | A proof that fusing measurements using point-to-hyperplane registration is invariant to relative scale[END_REF], demonstrated that the error function is not influenced by any scale parameter if the n-dimensional normals are normalized.

In fact, for this paper it was observed that the invariance to the scale parameter is due to the fact that the elements of λ λ λ have no influence on the direction of the estimated normals. For purposes of simplicity, this invariance will be shown for the 2D case here. Consider two sets of measurements taken at different times as

M * = [X * Y * ] and M = [X Y]
where each measurement (X, Y) has a different scale (e.g. centimeters and millimeters). The measurements are scaled to the same order of magnitude by the tuning parameter λ λ λ = diag(λ X , λ Y ). The equation of the tangent line ax + by + c = 0 (where a, b are the coordinates of the normal) at an i -th reference point can be obtained from the definition of the slope of a line:

y -λ Y Y * i x -λ X X * i = λ Y (Y * i+1 -Y * i ) λ X (X * i+1 -X * i ) (9)
that can be written as:

λ Y (Y * i+1 -Y * i )x + λ X (X * i -X * i+1 )y + c = 0, where c = λ X λ Y (X * i+1 Y * i - Y * i+1 X * i ). The normal a, b = λ Y (Y * i+1 -Y * i ), λ X (X * i -X * i+1
) can be represented for simplicity as:

-→

N * i = det (λ λ λ) • λ λ λ -1 V * i where V * i = [Y * i+1 -Y * i X * i -X * i+1 ] .
For simplicity this last will be written as:

V * i = [N * 1i N * 2i ] The projection of the point-to-point error ----→ M * i M i = λ λ λ(M i -M * i ) onto the normal direction -→ N * i
defines the distance of a point to a line (Figure 5(a)). It is clearly seen that the error function can be computed as:

e Hi = -→ N * i λ λ λ (M * i -M i ). Replacing -→ N * i into this error function as: e Hi = det (λ λ λ) • λ λ λ -1 V * i λ λ λ (M * i -M i ) allows to rewrite it as e Hi = λ X λ Y N * 1i (X i -X * i ) + N * 2i (Y i -Y * i )
, where is demonstrated that the tuning parameter λ λ λ has no effect on the minimization process since it has no influence on the direction of the normal and it scales its magnitude only.

The invariance for the 3D case has been demonstrated in [START_REF] Muñoz | A proof that fusing measurements using point-to-hyperplane registration is invariant to relative scale[END_REF] and it has been extended to ndimensions. In this paper, a better presentation of the proof will be given. The normal in three dimensions is obtained by performing the cross product between two hybrid vectors, scaled by λ λ λ = diag(λ X , λ Y , λ I ) (2D geometric points + intensity) as:

N * i = det(λ λ λ) • λ λ λ -1 V ik × V il
, where V ik and V il are defined as the k -th and the l -th closest point to M * i that lies on the reference dataset as Therefore, a 3D normal at the i -th reference point is obtained as follows:

V ik = λ λ λ (M * k -M * i ) and V il = λ λ λ (M * l -M * i ),
N * i =   λ Y λ I V ik Y V il I -V ik I V il Y λ X λ I V ik I V il X -V ik X V il I λ X λ Y V ik X V il Y -V ik Y V il X   =   λ Y λ I N 1 λ X λ I N 2 λ X λ Y N 3   (10)
The error function generated between two sets of 3D hybrid measurements can be defined then as:

e Hi = det(λ λ λ) V ik × V il (M * i -M i ).
The estimation of the normal presented in (10) can be easily extended to higher dimensions. The n-dimensional normal is estimated by performing the n-dimensional cross product between the n -1 vectors such as:

N * i = det(λ λ λ) • λ λ λ -1 V 1 × V 2 × • • • × V n-1 ∈ R n×1 (11)
The n-dimensional error function [START_REF] Muñoz | A proof that fusing measurements using point-to-hyperplane registration is invariant to relative scale[END_REF] between two sets of measurements can be defined as:

e Hi = det(λ λ λ)λ λ λ -1 V 1 × V 2 × • • • × V n-1 λ λ λ(M * i -M i )
, which can be re-written as:

e Hi = det(λ λ λ) V 1 × V 2 × • • • × V n-1 (M * i -M i ) (12)
where det(λ λ λ)

= λ 1 λ 2 • • • λ n .
The aforementioned normals can be computed by performing an n-dimensional cross product but other strategies can be equally used. In the Generalized-ICP strategy, the PCA (Principal Component Analysis) is computed. The eigenvector associated with its lowest eigenvalue is considered as the normal. Recently, an alternative solution has been provided in [START_REF] Segal | Generalized-ICP[END_REF], where the normals are efficiently and accurately computed by performing the Prewitt operator on projected spherical coordinates onto a spherical range image. This approach, however, only applies to the 3D case. For the 4-dimensional case presented in this article, a PCA analysis was performed to estimate the 4D normal, which can be written as:

N * i = λ Y λ Z λ I N * 1i λ X λ Z λ I N * 2i λ X λ Y λ I N * 3i λ X λ Y λ Z N * 4i
, where N 1 , N 2 , N 3 and N 4 are the components of the normal. These components can be estimated by considering the lowest eigenvalues of the nearest 4D points to a central i -th 4D point (Figure 6). Therefore, equation ( 6) can be rewritten for the 4-dimensional space as follows:

e Hi = det(λ λ λ) N Xi (X * i -X w i ) + N Yi (Y * i -Y w i ) + N Zi (Z * i -Z w i ) + N Ii (I * i -I w i )
where det(λ λ λ) = λ X λ Y λ Z λ I .

Results

In order to evaluate the Point-to-hyperplane ICP method, some parameters considered for the experiments are established. All the experiments were performed on both, real and synthetic RGB-D grayscale images in MATLAB. Furthermore, visual SLAM in real-time was performed in C++.

All the experiments were validated on a PC with Ubuntu 14.04, Intel core i7-4770K and 16GB ram.

A multi-resolution pyramid was used to improve computational efficiency (resolution: 160 × 120 at the top), where a pose is estimated at the top of the pyramid and the estimated transformation is employed to initializate the transformation in the next level until reaching the base of the pyramid.

The minimization process can be stopped by two criteria: an established maximum number of iterations (200 iterations for the experiments performed in this paper) or if the norm of the pose parameter is less than 1 × 10 -6 in rotation and 1 × 10 -5 in translation. To reject outliers, the Huber influence function was employed in only one M-estimator (as opossed to [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF][START_REF] Meilland | On unifying key-frame and voxel-based dense visual SLAM at large scales[END_REF] where the M-estimation is performed separately for color and depth). The M-estimation allows the use of different minimization functions not necessarily corresponding to normally distributed data.

The Point-to-hyperplane ICP method is compared in this paper with variants (which differ in how the uncertainty factors λ λ λ = diag(λ X , λ Y , λ Z , λ I ) are estimated) of the error function proposed in [START_REF] Meilland | On unifying key-frame and voxel-based dense visual SLAM at large scales[END_REF]. For this strategy, the classic Point-to-plane [START_REF] Chen | Object modeling by registration of multiple range images[END_REF] approach is employed to minimize the geometric term and a direct method for the photometric term as:

e Hi = λ λ λ   RR(x)N * i P m i -Π 3 TT (x) P * i I i w( TT(x), P * i ) -I * i (p * i )   ∈ R 4 ( 13
)
where the first and second row correspond to the geometric and photometric error function, respectively.

P m i ∈ R 3 is the closest 3D Euclidean point in the current cloud, R ← RR(x) is the incremental update of rotations, N * i = [N xi N yi N zi ]
are the normals of the reference points and Π 3 = [1, 0] ∈ R 3×4 is the projection matrix. For the purposes of this paper, the photometric term is minimized by using the Second Order Minimization (ESM) method [START_REF] Benhimane | Real-time image-based tracking of planes using efficient second-order minimization[END_REF].

A strategy to locally find the closest points is needed for computing the normals. The nearest neighbours to a central pixel in the image are considered to find its associated 3D point. At least three 4D points should be considered to estimate the 4D normal. A 3 × 3 window was considered in this paper for the experiments, which estimate the resultant normal of the 8 nearest neigbours (Figure 6). For the real-time application, however, a 2 × 2 window is employed to speed up the performance. It should be mentioned here that the computational cost is linear to the size of this window, but the accuracy is increased.

For the comparisons, the performance of the Point-to-hyperplane ICP method is compared with three different strategies that compute a non-adaptive or adaptive λ λ λ: 1) The intensities are normalized λ I = I i /255) (non adaptive), 2) an adaptive λ λ λ as in [START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF], where the scale parameter is the ratio between the Median Absolute Deviation (MAD) of the errors λ G = M AD(e I )/M AD(e G ) and 3) the covariance matrix of the metric errors as λ λ λ = cov(e G , e I ). For this last strategy, the T-distribution was employed to reject outliers as in [START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF]. The minimization of the error presented in [START_REF] Whelan | Elasticfusion: Dense slam without a pose graph[END_REF] will also be compared with a λ λ λ = eye(1) (λ λ λ is not estimated) in order to demonstrate that the parameter λ λ λ can improve the hybrid methods if it is well estimated. Alternatively, the estimation of the closest points were also done by searching a k d-tree (Labeled in Figure 7(b) as NN4D). This strategy demonstrated a better performance while aligning the frames when they are not close enough, but increasing the computational cost. A similar performance was observed in the 1000 synthetic frames that were equally tested. The label NN4D indicates that the nearest neighbours were obtained by using a k d-tree in the first iteration. The benchmarck datasets [START_REF] Handa | A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM[END_REF] and [START_REF] Sturm | A benchmark for the evaluation of RGB-D SLAM systems[END_REF] were used.

Simulated environments

The performance of the aforementioned strategies were compared in a synthetic environment.

The motivation for using synthetic data is that the generated images provide a groundtruth for the evaluation since the transformation between the frames is known. For the comparisons, 1000 synthesized RGB-D frames were generated with random poses and Gaussian noise was added.

Averages in time and number of iterations are shown in Table 1 2 . The time shown, however, does not consider the computation of the normals or construction of the k d-tree in the reference image. Therefore, for the employed RGB-D image in this experiment, the normals obtained 9.33 seconds for a 3 × 3 window and the construction of the k d-tree took 0.0056 seconds in MATLAB. For comparisons in this experiment, the matching points obtained by the k d-tree were used by considering 4D and 6D points. It was observed that the searching of the closest points by using this strategy reduces the number of iterations and convergence time if it is performed in the first iteration only. The chances to find the true nearest neighbours increases when more dimensions are considered. This is useful when the overlapping area between RGB-D frames is not large enough. However, the searching of the closest points in the k d-tree require extra computational time. Therefore, for purposes of this paper, only 4 dimensions were considered (This balances the computational cost and accuracy while estimating the pose). Figure 7(a) shows an example of the estimated trajectories in the convergence domain by different strategies.

Real environments

The well known living room ICL-NUIM RGB-D [START_REF] Handa | A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM[END_REF], freiburg1 and freiburg3 TUM [START_REF] Sturm | A benchmark for the evaluation of RGB-D SLAM systems[END_REF] benchmark datasets were employed to perform visual odometry to compare different hybrid strategies. For this experiment, a frame-to-frame alignment was employed. The estimated poses were used to evaluate the ATE (Absolute Trajectory Error) and RPE (Relative Pose Error). Various examples of the ATE evaluation for the Point-to-hyperplane ICP method are shown in Figure 8, where it can be seen that the Point-to-hyperplane ICP method can obtain close solutions w.r.t. the groundtruth without employing extra strategies for pose refinement or loop closure methods. [START_REF] Handa | A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM[END_REF]. It can be seen that the Point-to-hyperplane method [START_REF] Muñoz | Point-to-hyperplane RGB-D pose estimation: Fusing photometric and geometric measurements[END_REF] improves hybrid methods that combine the direct approach and the geometric Point-to-plane approach. The numerical results are shown in Table 2, where the methods are listed as follows:

( From Table 2, it can be seen that the Point-to-hyperplane methods improve other methods while obtaining less computational cost and less number of iterations. It was observed that when the frame-to-frame alignment is employed, the Strategies 2 and 4 obtain about the same results as Strategies 1 and 3, respectively. Therefore, the results for these strategies are shown together in Table 3, where it can be noted that the adaptive λ methods obtained less ATE and RPE error for the freiburg3 sequences (benchmark structure vs texture), however the computational cost is high w.r.t. Point-to-hyperplane strategies. It can be noted that the Point-to-hyperplane ICP method obtained more robust results in challenging 360 degree scenarios.

The results obtained by performing visual odometry in the synthetic environment demonstrated the robustness of the Point-to-hyperplane ICP method when rich color and depth features can be associated. During the experiments in real scenarios, the Point-to-hyperplane ICP method obtained better estimations in challenging sequences with blurred images such as fr1/room and fr1/360 (closed loops), demonstrating the robustness of the method. It was observed in the experiments that adaptive methods can improve the accuracy of the pose estimation methods when rich geometric and photometric information is available (as the case of sequences fr/3), however they were not robust enough for closed loop sequences. On the other hand, non adaptive methods achieved faster convergence and they obtained better estimations for sequences as fr1/desk, fr1/desk2 (where the sequences contain several sweeps) and fr1/floor, where the geometric features are not significant enough but texture. An example of the performance of the Point-to-hyperplane method can be shown in Figure 10, where the 3D reconstruction of loop closed sequences is obtained by transforming the cloud of points by the estimated 6DOF pose parameter. The Point-to-hyperplane ICP method has been implemented successfully for real time applications by employing the ASUS Xtion sensor. The proposed method here has demonstrated its robustness in a long corridor which contain similar geometric and photometric features (See Figure 9). The depth map obtained by the sensor contain holes in the depth map which are not valid values. This can generate a problem in the estimation of the normal for the Point-to-hyperplane ICP. In previous works, the corresponding intensities with non-valid depth values were considered as outliers. Another solution is to interpolate the surrounding valid depth values or assign the value of the closest point. By doing this, the method can achieve the main advantages of the Point-to-plane ICP or the direct method if either the color or the 3D point are not available.

The implementation of the Point-to-hyperplane ICP method in real-time can be an alternative for performing visual navigation, 3D reconstruction and localization of robotics platforms. The method can be also used for combining other types of measurements obtained by different sensors.

Conclusion

In this article, extended results of previous works on the Point-to-hyperplane strategy were shown. Particularly, the method has been extended to higher dimensions (3D Euclidean points + intensity and 3D Euclidean points + 3 channels of color) and it was mathematically demonstrated that the joint error function projected onto the normal direction has the effect of canceling out the λ tuning parameter since it does not change the direction of the normal. The future aim is to exploit other types of measurements for estimating the pose.

Two strategies for obtaining the closest points were compared, k d-trees and bilinear interpolation. For the benchmark sequences presented, both strategies obtained about the same performance. The bilinear interpolation was employed for the real-time application where a 2 × 2 window, which has been used to improve the computational cost. The real-time visual SLAM is running under CPU, obtaining a mean computational time of 1236 ms for the estimation of the normals in 4 dimensions. As a future work, the proposed method here will be implemented using GPU and a refinement method will be employed by estimating the global pose for a keyframe-to-frame and keyframe-to-keyframe tracking.

The reconstructions can be refined by any post-processing algorithms. The post-processing refinement strategies were not introduced here, but strategies that perform global convergence can be considered [START_REF] Yang | Go-icp: Solving 3d registration efficiently and globally optimally[END_REF][START_REF] Campbell | GOGMA: globally-optimal gaussian mixture alignment[END_REF][START_REF] Straub | Efficient globally optimal point cloud alignment using bayesian nonparametric mixtures[END_REF]. The refinement of the Point-to-hyperplane method has been recently performed by estimating the global pose of a RGB-D frame w.r.t. the generated 3D model in [START_REF] Muñoz | Generalized global rgb-d registration: Collaborative local and global pose estimation by fusing color and depth[END_REF]. In the first column the groundtruth obtained by an external motion capture system is shown, the second column shown the result of the Point-to-hyperplane method. This difficult 360 degree sequence with motion blur shows that the proposed method can achieve more robust estimations.

Figure 1 .

 1 Figure 1. Classic hybrid approach diagram. The geometric and photometric error (e G and e I , respectively) are jointly minimized. A tuning parameter λ λ λ = (λ G , λ I ) weights the contribution of each measurement during the minimization process. The metric measurements are represented in a i -th single vector M i = P i I i ∈ R 4 , which contains the 3D Euclidean points P and its associated intensities I.

Figure 2 .

 2 Figure 2. Examples of RGB-D images: Texture vs Color. (a) and (b) are images taken from an environment where the geometric measurements are more significant than photometric measurements while (c) and (d) were taken from a scenario where the texture is more significant than geometry. The depth image is visualized as an intensity image, where black indicates a non valid depth measurement.

the 4 ×

 4 4 pose matrix onto the 3 × 4 space and e 3 = [0 0 1] extracts the depth component of a transformed 3D point.

Figure 3 .

 3 Figure 3. Influence of the scale coefficient on the error function residuals (Equation (13)) when (a) λ is not estimated, (b) when the intensities and 3D points are normalized between [0, 1] [18] and (c) when λ is estimated at each iteration[START_REF] Tykkälä | Direct Iterative Closest Point for Real-time Visual Odometry[END_REF]. An example of a cost function is shown in (d), which indicates that the number of iterations can be improved when a good choice of λ is made.

Figure 4 .

 4 Figure 4. Point-to-hyperplane approach diagram.In the first iteration, the closest points can be estimated by finding in a k d-tree with the 4D vector and linear interpolation in the following iterations. The method is invariant to a tuning parameter λ.

FebruaryFigure 5 .

 5 Figure 5. Principle of the Point-to-hyperplane applied in 2 dimensions Point-to-line and 3 dimensions Point-to-plane. It should be noted that the axis x and y are not in the same scale in (a). A tuning parameter λ is added in order to demonstrate the invariance of the Point-to-hyperplane method in the 2D case. The distance e H i is the result of projecting the vector ----→ M * i M i onto the normal direction N * i .

  respectively (See Figure5(b)).

Figure 6 .

 6 Figure 6. The neighbours of a central pixel (in grey) in a 3×3 window on a image are employed to associate its closest 3D points and estimate the normal to the connecting plane (or hyperplane). Any pixel on the image can be selected as a central pixel to compute the distance with its surrounding neighbours (except for the corners).

FebruaryFigure 7 .

 7 Figure 7. Example of the estimated trajectories between a current and a reference frame in the transformation space. (a) The green and red dot indicates the initial and final pose respectively. The Point-to-hyperplane method improves the other hybrid methods by obtaining more direct trajectories and a less number of iterations (b).A similar performance was observed in the 1000 synthetic frames that were equally tested. The label NN4D indicates that the nearest neighbours were obtained by using a k d-tree in the first iteration.

Figure 8 .

 8 Figure 8. Examples of the Absolute Trajectory Error evaluation obtained by the Point-to-Hyperplane method. The benchmarck datasets[START_REF] Handa | A benchmark for RGB-D visual odometry, 3D reconstruction and SLAM[END_REF] and[START_REF] Sturm | A benchmark for the evaluation of RGB-D SLAM systems[END_REF] were used.

1 )

 1 Hybrid + non-adaptive λ I = I i /255 (2) Hybrid + non-adaptive λ I = I i /255 (NN4D* 3 ) (3) Point-to-hyperplane (4) Point-to-hyperplane (NN4D) (5) Hybrid + adaptive λ G = M AD(e I )/M AD(e G ) [15] (6) Hybrid + adaptive λ λ λ = cov(e G , e I )[START_REF] Kerl | Dense visual SLAM for RGB-D cameras[END_REF]

  Advanced˙Robotics˙Journal˙2017

Figure 9 .

 9 Figure 9. 3D reconstruction of a long corridor by the Point-to-hyperplane method

Figure 10 .

 10 Figure 10. 3D reconstruction of sequences freiburg1: room, 360, plant and teddy (shown at each row, respectively). In the first column the groundtruth obtained by an external motion capture system is shown, the second column shown the result of the Point-to-hyperplane method. This difficult 360 degree sequence with motion blur shows that the proposed method can achieve more robust estimations.

Table 1 .

 1 Averages in time and in the number of iterations until convergence for 1000 synthesized Images at Random Poses. The legend NN4D or NN6D indicates that the closest points were estimated in the first iteration only by searching the nearest neighbours in the 4D or 6D k d-tree.

	Method	# Iterations Time (sec)
	Hybrid (λ λ λ = ones)	157.668	2.046
	Hybrid + non-adaptive λ I = I i /255	124.419	1.598
	Hybrid + non-adaptive λ I = I i /255 (NN4D)	116.609	1.563
	Hybrid + adaptive λ G = M AD(e I )/M AD(e G ) [15]	154.966	2.010
	Hybrid + adaptive λ λ λ = cov(e G , e I ) [11]	155.455	6.079
	Point-to-hyperplane (3D points + Intensity)	48.038	0.531
	Point-to-hyperplane (NN4D)	13.224	0.191
	Point-to-hyperplane (3D points + RGB)	96.79	2.1572
	Point-to-hyperplane (NN6D)	79.439	1.7978

Table 2 .

 2 Averages in Time (miliseconds), number of iterations, Relative Pose Error (RPE) and Absolute Trajectory Error (ATE) for the synthetic dataset

Table 3 .

 3 Average in Time (miliseconds), number of iterations, Relative Pose Error (RPE) and Absolute Trajectory Error (ATE) for the dataset freiburg1 and freiburg3[START_REF] Sturm | A benchmark for the evaluation of RGB-D SLAM systems[END_REF].

	Sequence	Method	ATE (m) RMSE MEAN RMSE RPE translational (m) RPE rotational (deg) MEAN RMSE MEAN	AVERAGE Time(sec) #Iterations
		1 & 2	0.068	0.064	0.041	0.035	2.485	2.068	0.227	18.41
	fr1/xyz	3 & 4 5	0.045 0.092	0.038 0.087	0.021 0.040	0.019 0.035	1.106 2.407	0.998 2.034	0.301 0.364	26.94 28.97
		6	0.086	0.080	0.041	0.036	2.421	2.080	0.946	24.51
		1 & 2	0.102	0.087	0.040	0.034	2.918	2.566	0.273	21.16
	fr1/rpy	3 & 4 5	0.035 0.129	0.032 0.111	0.038 0.045	0.032 0.037	2.820 3.034	2.652 2.643	0.445 0.288	36.79 23.15
		6	0.131	0.114	0.046	0.038	2.947	2.594	0.997	25.09
		1 & 2	0.353	0.332	0.111	0.090	3.917	3.550	0.295	22.95
	fr1/360	3 & 4 5	0.322 0.190	0.296 0.179	0.152 0.094	0.114 0.085	3.159 4.113	2.859 3.583	0.460 0.285	38.68 22.18
		6	0.268	0.245	0.191	0.149	5.210	4.539	0.962	24.55
		1 & 2	0.375	0.353	0.075	0.055	3.373	2.836	0.255	19.75
	fr1/room	3 & 4 5	0.152 0.323	0.131 0.286	0.056 0.063	0.047 0.051	2.673 3.250	2.329 2.743	0.375 0.308	33.36 23.34
		6	0.363	0.305	0.068	0.055	3.434	2.902	0.896	23.41
		1 & 2	0.064	0.060	0.043	0.036	2.738	2.403	0.236	19.78
	fr1/desk	3 & 4 5	0.071 0.069	0.067 0.065	0.044 0.044	0.036 0.036	2.310 2.580	2.028 2.233	0.408 0.300	34.98 23.79
		6	0.065	0.062	0.045	0.037	2.667	2.285	0.877	24.19
		1 & 2	0.083	0.079	0.058	0.049	3.816	3.133	0.243	20.00
	fr1/desk2	3 & 4 5	0.133 0.560	0.116 0.233	0.060 0.624	0.051 0.154	3.026 24.448	2.641 7.055	0.496 0.328	38.33 25.02
		6	0.409	0.188	0.463	0.129	16.883	5.887	0.940	25.20
		1 & 2	0.124	0.105	0.035	0.023	1.946	1.364	0.198	15.61
	fr1/floor	3 & 4 5	0.473 0.273	0.405 0.224	0.080 0.085	0.051 0.037	3.909 2.846	1.915 1.754	0.355 0.281	31.67 21.64
		6	2.050	1.765	0.384	0.089	21.096	4.932	0.873	22.75
		1 & 2	0.066	0.056	0.035	0.030	1.740	1.568	0.262	20.63
	fr1/plant	3 & 4 5	0.101 0.067	0.093 0.055	0.055 0.031	0.043 0.027	2.130 1.608	1.947 1.405	0.395 0.329	34.88 25.48
		6	0.065	0.054	0.033	0.028	1.623	1.427	0.914	23.52
		1 & 2	0.260	0.219	0.061	0.046	1.996	1.656	0.282	20.97
	fr1/teddy	3 & 4 5	0.169 0.303	0.158 0.267	0.070 0.086	0.056 0.049	2.287 2.432	1.954 1.734	0.424 0.322	36.86 23.74
		6	0.339	0.299	0.100	0.055	2.681	1.821	0.926	24.12
		1 & 2	0.136	0.133	0.028	0.025	0.936	0.856	0.201	16.71
	fr3/s t far	3 & 4 5	0.401 0.044	0.361 0.040	0.067 0.024	0.061 0.021	1.561 0.707	1.428 0.639	0.172 0.274	15.52 22.32
		6	0.044	0.042	0.021	0.018	0.649	0.590	0.617	16.32
		1 & 2	0.152	0.143	0.027	0.021	1.568	1.176	0.161	13.59
	fr3/s t near	3 & 4 5	0.185 0.056	0.157 0.052	0.045 0.018	0.036 0.016	1.862 0.965	1.525 0.848	0.240 0.345	21.32 27.43
		6	0.066	0.063	0.019	0.016	0.997	0.856	0.813	21.25

Note that in (4) the parameter λ G = det(λ G λ G λ G ) is a scalar for the geometric point-to-plane error

In order to better present the performance of the strategies, the best obtained results are displayed in bold in all tables.

The legend NN4D means that the closest points were estimated by a k d-tree in the first iteration only

Funding

This work is supported by the European H2020 project: COMANOID, Université Côte d'Azur, CNRS, I3S, France and CONACYT (Consejo Nacional de Ciencia y Tecnología) under grant [265807], Mexico.