Franziska Schmidt

Sogol Kharrazi

Sigurdur Erlingsson

Carl Van Geem

Xavier Cocu

Bernard Jacob

FALCON II. : Input for a european PBS definition : review of vehicle legislations and infrastructure design criteria

Keywords: Performance Based Standard (PBS), Smart Infrastructure Access Policy (SIAP), vehicle policy, infrastructure design criteria

published or not. The documents may come

FALCON II. : INPUT FOR A EUROPEAN PBS DEFINITION: REVIEW OF VEHICLE LEGISLATIONS AND INFRASTRUCTURE DESIGN CRITERIA

Vehicle Policy

As part of the FALCON project, relevant European legislation for commercial heavy vehicles, and the corresponding regulations implemented in the countries involved in the FALCON project, namely Sweden, Norway, Netherlands, Germany, UK, France and Belgium were reviewed and compared. It should be noted that there are two types of European legislations: regulations and directives. Regulations have general application and are applicable in all member states, while directives set out general rules to be transferred into national law by each country as they deem appropriate.

Here a brief comparison of the length and weight limits in the studied European countries is presented, more information about the vehicle policy can be found in the deliverable 3.3 of the FALCON project. The length limit of motor vehicles in the EU is regulated in the R (EU) No 1230/2012 which is also applied in the studied countries [6]. However, in Norway it does not apply to timber transport, and in Sweden it is only applicable to modular vehicles. The length of vehicle combinations in Europe is regulated in the Directive 96/53/EC, with a limit of 16.5 m for articulated vehicles and 18.75 m for road trains. However, article 4 of the directive gives each member country the possibility to use longer vehicle combinations for national transport, if they are based on the modular system. A modular combination is a vehicle combination that consists of vehicle units defined in Annex I of the directive [7].

In Belgium, UK, France, Germany and The Netherlands, the European length limits apply. But Sweden has a length limit of 24m which is also the limit for timber transport in Norway, otherwise limits of 17.5 m and 19.5 are used for articulated vehicles and road trains respectively in Norway. European Modular System (EMS) combinations, which are 25.25 m combinations of EMS units are allowed in Sweden, Norway, Belgium, the Netherlands, and in 13 of the 16 provinces (Länder) in Germany.

The single axle load limits are very similar in the studied countries complying with the EU limits for international traffic stated in the Directive 96/53/EC: 11.5 t for driving axles and 10 t for other axles; however, France has a higher axle load limit of 13 t. The load limits for a bogie are also comparable and mostly complying with the Directive 96/53, but the reference axle distances for setting the bogie load limit are slightly different for some countries. For instance, in Norway 0.8 m and in France 0.9 m is used as the reference axle distance, below which the lowest load limit is applied, while in other countries 1 m is used which is the same as the EU regulations for international traffic. It is a similar case with triple axles loads, i.e. the load limits are comparable, but the reference axle distances are not uniform. Again, France allows higher load limits, e.g. up to 31.5 t on a tridem instead of 27 t, and Norway has the lowest load limit for an axle distance below 1 m.

The weight limit for a motor vehicle depends on its number of axles in all the considered countries and is quite similar to the European limits for the international traffic stated in the Dir 96/53/EC , the Netherlands is an exemption with higher limits. For regulation of the weight limits of trailers and semi-trailers, different approaches are used in each country. Commonly the weight limits are regulated based on features such as the axle distances, number of axles and the vehicle type. For instance, in Sweden the weight limit depends on the axle distance between the foremost and rearmost axles in the vehicle/vehicle combination, while in the Netherlands, the axle load limits and the total weight limit of the vehicle combination determine the weight limits on the constituent units, see Table 1.

Infrastructure Catalogue

Pavements, pavement structures

The infrastructure catalogue developed within the FALCON project is the basis for the development of both the PBS and the SIAP (Smart Infrastructure Access Programme): indeed, on the one hand, the infrastructure catalogue distinguishes the various roads or road networks with various access policies; on the other, this catalogue shows the infrastructure elements for which the design criteria must be determined.

A series of assumptions have been made to establish this infrastructure catalogue: Only design values are considered regarding axel loads and tyre pressures, physical parameters such as layer thicknesses, material characteristics and ambient climate. Therefore, the considered infrastructure is supposed to be in a design (meaning nominal) state. Moreover, a second assumption assumes that only current design codes are considered. This catalogue deals with three infrastructure types: pavements, bridges and tunnels. Geometrical and mechanical (load capacity) aspects are considered.

Pavements are multi-layer structures built on top of the subgrade soil. Their main structural function is to support the axle loads, and spread them downwards to the subgrade (natural ground) avoiding overstressing of all layers, that is reducing the induced stresses and strains to tolerable levels during all climatic conditions. A typical pavement structure consists of:  a surface course, which can be divided into a wearing course (or top layer), a binder course and a bound road base.  Unbound base layers, usually divided in two layers: a base course and a subbase, both build of materials with good mechanical resistance.  a pavement foundation: the subgrade (native soil), sometimes topped with a capping layer. The main functions of the surface course is to provide good pavement surface characteristics (evenness, skid resistance), to ensure a good rolling quality and appropriate safety conditions for road users, and to protect the road base from the wear due to traffic loads, and from the penetration of water (impermeability). The main functions of the base layers are to withstand the loads induced by traffic (fatigue resistance), and to distribute the stresses on the pavement foundation. The capping layer protects the subgrade during the works phase, improves the homogeneity and bearing capacity of the natural subgrade, and protects it from frost actions. European pavements can be categorized according to the materials used in the different layers, namely flexible pavements (with bituminous bound layers or asphalt concrete) or rigid pavements (with Portland cement concrete layers). Flexible pavement flex under traffic loading giving a localized deformation bowl under the load that is further distributed downwards over HVTT15: Review of vehicle legislations and infrastructure design criteria HVTT15: Review of vehicle legislations and infrastructure design criteria gradually increasing area that is proportional to the stiffness's of respectively layer. Rigid pavement structure's load transmission is different from flexible structures. It relies on the rigid slab action that spreads the loading over a large area. Further are pavement structures with both bituminous and concrete bound layers usually referred to as semi-rigid structures, see Figure 1. The various layers may have different thicknesses depending on the country considered. Thin flexible pavements (see figure 1a) consist of a relatively thin bituminous surface course (or asphalt concrete layer) resting on one or more layers of unbound granular materials. They are usually the most economical, but due to the low stiffness of the unbound granular base layers, they are only suitable for low to medium traffic levels. They are frequently designed to last for 20 years. The other types of pavements in Figure 1 are more suited for heavier traffic levels on main arterials roads and motorways. Concrete pavements are usually designed to last longer (30 to 40 years) than the tick-bituminous or semi-rigid pavements (20 to 30 years), they are usually more expensive to build but they need less maintenance.

Two damage phenomena are commonly studied in flexible pavement design: That is fatigue cracking (bottom-up cracking) starting at the bottom of the lowest bituminous bound layer and rutting that is accumulation of plastic deformation in the pavement structure. Fatigue is the phenomenon of damage induced by repeated applications of small traffic load actions. Fatigue laws express the relation between the number of applications of loads and tensile strain or stress leading to failure. Only bituminous pavements are subject to rutting. Rutting is the manifestation of accumulation of contribution of permanent deformation in the different layers of the pavement structure. Thus, all layer can contribute to the rutting development. Rutting design calculations are usually associated with the repeated induced vertical strain at the top of the subgrade material. From extensive rutting tests performed in different European countries with different tyres, tyre configurations, axle loads, inflation pressures, etc. (e;g. [START_REF]Effects of Wide Single Tyres and Dual Tyres[END_REF]) a tyre configuration factor (TCF) was defined. The TCF value relates the pavement wear of a given tyre to the pavement wear of a reference tyre. Within different axle categories (steered, driven or towed axle), there is a wide range of TCF values. For rigid pavements the main steps in the structural thickness design are associated with calculations of critical stresses at three locations that are related to crack initiations due to traffic loading. These three locations are at the interior (center), the edge and corner of the proposed concrete slab. Usually the tensile stress at the top of the slab when loaded at the corner of the slab is the most critical one.

Road managers design pavement structures using a wide range of parameters such as the expected transport needs, the available materials, and climate conditions. Hence, the pavement structures vary from country to country and even along the road network managed by one road authority leading to a wide variety of pavement structures. For FALCON, the following main factors need to be considered in order to establish a representative "catalogue", or a representative library of pavement structures:  The type of pavement structure: thin flexible, thick bituminous, …, or concrete pavement,  The level of traffic, expressed by the number of heavy vehicles (HV) per day and ESALs (Equivalent Single Axle Loads),  The chosen level of service, which can be expressed by a factor of safety, or a risk of failure,  The bearing capacity of the subgrade,  The mechanical characteristics of the pavement materials.  The ambient climate Bridges can be described by influence lines, i.e. transfer functions from unit loads to the considered effects (mainly bending moments and shear forces at supports). Therefore, the main characteristics are the number of supports and the various span lengths.

Bridges

Bridges are classified using various criteria: depending of the type of structure (suspended, cable-stayed, …), of the materials (reinforced or prestressed concrete, steel, composite, …), of the number of spans (one, two, …), of the supports and the degrees of freedom (isostatic, hyperstatic), … But when comparing the effect of one given type of vehicle on a bridge with the effect of another vehicle, it is sufficient to calculate the convolution of the characteristics of the vehicles (axle loads, distances between axles) with a structural information of the structure called "influence line". The influence line of a given effect (for example bending moment, shear stress, …) is a function which gives this effect when a unit force is applied to the structure. Therefore, for the bridge catalogue of the FALCON a set of influence lines have been chosen, namely:

 One-span bridge, bending moment at midspan and shear at support, span lengths of 10~m, 20 m, 35 m, 50 m and 100 m,  Two-span bridge, bending moment at midspan, at central support and shear on central support, both span lengths of 5 m, 10 m, 17.5 m, 25 m and 50 m (see Figure 2).

Road geometry

The geometry of the road is an important point when talking about dimensions of trucks. Indeed, trucks have to be compliant, in a static and a dynamic way, with the existing infrastructure, meaning the roundabout, the right turn, left turn, the ramps,etc. Road geometry is usually defined following several parameters: -the horizontal alignment, consisting of straight sections connected by circular horizontal curves (themselves defined by their radius and length) and transition curve (horizontal curve with a varying radius to connect other elements); -the longitudinal profile, characterised by grades (road slopes) connected by parabolic vertical curves (sags or crests) used to provide a gradual change from one road slope to another for a smooth vehicle navigation; -the cross section, i.e. the number of lanes, their widths and cross slopes, as well as the presence or absence of various features like shoulders, curbs, sidewalks, drains and ditches -other various elements like the clearance under bridges and in tunnels, etc.

Road geometry is strongly related to vehicle dynamics. The stability of vehicles depends on the observance of rules linking vehicle speed, radius of curves and cross slope. Road alignment also influences sight distance, which is an important factor for safety. While studies of contributing factors to road accidents show that human factors predominate, roadway factors are the second most common category of elements that influence the road safety performance.

For road design and depending on the road category (typically motorways, express roads, interurban arterials and ordinary roads), there are target or minimal acceptable values of geometrical parameters pertaining the curve radius, transition curve length, longitudinal profile, width of traffic lanes (the normal traffic lane width being 3,50 m) and shoulder, cross slope. Nominal values of some geometrical characteristics are provided in Table 2.

The geometrical characteristics of highway exit/access lanes, at-grade intersections and roundabout are also regulated to allow smooth left-or right-turn maneuvers, and crossing of roundabouts, see Figure 3.

Infrastructure Design Criteria

For the infrastructure elements included in the infrastructure catalogue, European design criteria have been compared. It has been observed that while the National design criteria are quite similar, differences exist.

Pavements, pavement structures

Most countries use mechanistic-empirical pavement design methods, which are similar in their principle. They are based on two main steps:

 A calculation of the stress-strain response of the pavement based on a reference load (generally defined as the "equivalent standard axle load", or ESAL), using a multi-layer linear elastic pavement model.  The application of several pavement design criteria, which allow to calculate the number of standard axle loads (ESALS) which can be supported by the pavement before failure (also called the pavement life), in function of the maximum level of stress or strain calculated in each pavement layer.

In the design process the accumulated traffic loads during service life is converted into a number of Equivalent Standard Axle Loads (ESAL) that the pavement must support. The bearing capacity of the subgrade is expressed by its elastic modulus and the mechanical properties of the pavement materials comprise the elastic modulus, Poisson ratio and fatigue properties. The ambient climate, and in particular the temperature (single value, or several climatic periods) is considered for the bituminous materials, which exhibit temperature-dependent behaviour. Sometimes a factor of safety (or risk coefficient) is used to adjust the number of loads to failure. Based on the above-mentioned parameters the required layer thicknesses are calculated so the pavement can last for the required designed life.

Several of these design parameters are country-specific:  The design lifetime of the structure is generally 20 years in Europe, but some countries may use up to 30 years (Belgium).  The traffic volume and loads vary from one country to another, depending on the local economy and the geographical location. Moreover, the definition of ESAL (Equivalent Standard Axle Load) also slightly differs: generally, ESAL is given in equivalent 100 kN-axle loads (for a given configuration of tyres with given air pressure), but the French ESAL is expressed in equivalent 130 kN-axle loads.  The fatigue and rutting criteria in the pavement structure are not expressed with the same formula and at the same location.  Climatic parameters are also different from one country to another: obviously, the targeted temperature varies by countries, but also by region within a given country: e.g. in Sweden, this temperature is defined according to five climatic regions and six seasons per year. Frost-thaw can therefore be taken into account.  Materials used in the pavement layers in different countries have different characteristics and thus the commonly used values for their elastic modulus and Poisson ratio as well as their fatigue properties are different.

Thus, there is no unique and standardized approach in Europe to design pavements. However, there is a general agreement about the main steps in the design framework that is the design approach and design criteria's. The design tool and the input parameters values for the mechanistic-empirical pavement design can differ for different regions or countries throughout Europe. However, this common ground in pavement design can be exploited in the development of a vehicle policy framework and a PBS Bridge design is standardized by the Eurocodes. The Eurocode 1 specifies the actions on structures, such as wind loads, traffic loads [START_REF] En | Eurocode 1: Actions on structures -Part 2: Traffic loads on bridges[END_REF], etc. These actions are described by load models, to be applied for the bridge (and other structures) design. Several traffic load models exist for calculation of extreme effects. Specific fatigue load models (five) are proposed for lifetime assessment, the simplest one being the most conservative and the most detailed being the less conservative. Some differences of the load model intensities are allowed by countries, by the application of α-factors.

Bridges

Bridge design is standardized in Europe through codes called Eurocodes: for example, Eurocode 2 is the code dedicated to the design of bridges and civil engineering structures in concrete (plain, reinforced, prestressed). Eurocode 3 is dedicated to steel construction. The traffic actions to take into account when designing bridges are given in Eurocode 1 [START_REF] En | Eurocode 1: Actions on structures -Part 2: Traffic loads on bridges[END_REF]. More precisely load model 1 ("LM1") of Eurocode 1 defines the actions for extreme load calculations, whereas load model 3 ("LM3") are those for calculation of fatigue damage. LM1 is composed of a uniformly distributed load, whose numerical values can be different between the traffic lane (slow or fast lane), and punctual vertical corresponding to axle loads. These values are nationally specific through the α-factors (see Figure 4 and So, while standardizing the design and the building of bridges in Europe, national coefficients make it possible to adapt the load models to national requirements. It should be mentioned here that it is not possible to compare at a glance the various national load models (meaning that the effect on a given structure would depend on this infrastructure).

Road geometry

For road geometry, the design criteria are country-specific and the target values for the parameters (as explained in Section 2) are quite similar from one country to another. Therefore, in the FALCON project, only the cases of one main road and one secondary road are taken into account.

The design criteria can then be applied to vehicles included in the vehicle policy framework in order to determine the PBS.

Conclusions and Discussion

The work, presented here and in deliverables D3.2, D3.3 and D3.4 of the FALCON project, is the background for the development of Performance Based Standards, and then the development of SIAP. Indeed, the information on vehicle policy and different types of infrastructure allow determining which impact of vehicles on infrastructure is sustainable. This will give more flexibility for vehicle design and usability, and therefore for innovation in vehicle industry, without inducing increased damage on the existing road infrastructure, which is ageing and whose limited maintenance budget is decreasing.

Figure 1 :

 1 Figure 1: Typical pavement structures.

Figure 2 :

 2 Figure 2: Bending moment at midspan ((red curve) and on support (blue curve), for two-span structure with both span lengths equal to 50 m.

HVTT15:

  Review of vehicle legislations and infrastructure design criteria HVTT15: Review of vehicle legislations and infrastructure design criteria

Figure 4 :

 4 Figure 4: Load model 1 of Eurocode 1: on each lane (lanes 1, 2 and 3), an uniformly load and several punctual loads are applied. .

Table 1 . Vehicle weight limits (ton)

 1

		EU International	Sweden	Norway	Netherlands	Germany	France	UK	Belgium Flanders & Wallonia
	Motor	18/25(26) 1 /32	18/25(26) 1 /31(32) 1	19/26/26-32	21.5/28-	18/25(26) 1 /32	19/26/32	18/25(26) 1 /30(32) 1	19/26/32
	vehicle	2/3/4+ axles	2/3/4+ axles	2/3/4+ axles	31.5/34(37) 1	2/3/4+ axles	2/3/4+ axles	2/3/4 axles	2/3/4 axles
					2/3/4+ axles
	Trailer	18/24	GVW/GCW table	10/18,20/24,27	Depends on the	18/24	19/26	18/24	10/18/24
	Semitrailer	2/3 axles	for axle distance	1/2/3 axles	axle distance and	2/3 axles	for 2/3 axles	for 2/3 axles	1/2/3 axles
				ST or CT	number of axles,	Trailer			Trailer
				20/28/30 1/2/3 axles FT or DY-ST	see Erreur ! Source du renvoi introuvable.				22-44 Semitrailer
	Vehicle	36/40	64	50	50	28/36/40(44) 3
	combination	4/5 axles Road train	GVW/GCW table for axle distance	GCW table for axle distance	60 EMS	3/4/5 axles Road train
		36(38) 2 /40(44) 3		60		28/
		4/5 axles Articulated vehicle		EMS & timber

36(38) 2 /40(44) 3

  If the semitrailer axle distance is bigger than 1.8m and the driving axle is fitted with twin tyres and air suspension3 If carrying a 45-feet ISO container, 42t for if the motor vehicle has two axles and 44t for if the motor vehicle has three axles4 If the single axle load does not exit 12t 5 With air suspension CT=Centre Axle Trailer, FT=Full trailer, ST=Semitrailer, TK=Truck It should also be noted that in the R (EU) No 1230/2012 and the Dir 96/53/EC, there are extra criteria that indirectly impose restrictions on the dimensions and load distribution of the vehicle to ensure maneuverability and traction (EC 2012). Examples of such criteria are the swept area in a roundabout, ratio of the load on steer or drive axles, and engine power based on the vehicle weight.

		38/40(44) 4	26/36/40	29/35
		4/5 axles	3/4/5 axles	TK2-CT1/2+
		Road train	Road train	36/42(44) 5
		38/40(44) 4	26/36(38) 2 /40(44) 3	TK3-CT1/2+
	3/4/5 axles Articulated vehicle	4/5 axles Articulated vehicle	3/4/5 axles Articulated vehicle	39/44 4/5 axles
				Other road trains
				29/39/43(44) 5
				3/4/5+ axles
				Articulated vehicle
				60
				EMS

1

If driving axle is fitted with twin tyres and a) air suspension (or equivalent) or b) drive axle load does not exceed 9.5 t 2

Table 3 )

 3 .

Table 3 : Comparison of α-factors in several European countries (extracted from deliverable D3.4 of FALCON project).

 3

	Bridges	Norway	Netherlands Germany	France	UK	Belgium
				𝛼 𝑄1
				= 0.8	LM1:
	α-factor	Alpha-factors (LM1 of Eurocode 1): 𝛼 𝑄1 = 0.6 All the other alpha-factors= 1.0.	αq1 = 1,15 and for i >1 αqi = 1,40	𝛼 𝑄2 = 0.8 The other alpha-factors 𝛼 𝑄𝑖 equal to 0.	𝛼 𝑄1 = 1 𝛼 𝑄𝑖 = 1, ∀𝑖 ≥ 2 𝛼 𝑞1 = 1 𝛼 𝑞1 = 1.2, ∀𝑖 ≥ 2	LM1: ∀𝑖, 𝛼 𝑄1 = 1 𝛼 𝑞1 = 0.61 𝛼 𝑞𝑖 = 2.2, ∀𝑖 ≥ 2	For new bridges: European class 1 (alpha-factors = 1) For existing bridges: Belgian class 2 (alpha-factors in general = 0,8).
				∀𝑖, 𝛼 𝑞𝑖	𝛼 𝑞𝑟 = 1.2
				= 1


Nathan Crombez
email: nathan.crombez@utbm.fr

El Mustapha Mouaddib
email: mouaddib@u-picardie.fr

Franc Guillaume Caron
email: guillaume.caron@u-picardie.fr

¸ois Chaumette

Franc ¸ois Chaumette
email: fran-cois.chaumette@inria.fr

Visual Servoing with Photometric Gaussian Mixtures as Dense Feature

Keywords: Visual servoing, Photometric Gaussian Mixture, large convergence domain

The direct use of the entire photometric image information as dense feature for visual servoing brings several advantages. First, it does not require any feature detection, matching or tracking process. Thanks to the redundancy of visual information, the precision at convergence is really accurate. However, the corresponding highly nonlinear cost function reduces the convergence domain. In this paper, we propose a visual servoing based on the analytical formulation of Gaussian mixtures to enlarge the convergence domain. Pixels are represented by 2D Gaussian functions that denotes a "power of attraction". In addition to the control of the camera velocities during the servoing, we also optimize the Gaussian spreads allowing the camera to precisely converge to a desired pose even from a far initial one. Simulations show that our approach outperform the state of the art and real experiments show the effectiveness, robustness and accuracy of our approach.

I. INTRODUCTION

V ISUAL servoing is a closed-loop control method for dynamic systems, such as robots, which uses visual information as feedback. Visual information is usually obtained by a digital camera directly set into motion by the system i.e. the eye-in-hand configuration. In another way, the camera may be stationary, external to the system in motion. Visual perception can be made with one or several cameras. In our case, we are particularly interested in eye-in-hand configuration with one camera. More precisely, the control law is based on visual information acquired by a camera in order to move a robot effector to a desired pose, implicitly defined by an image.

The visual information extracted from images are usually geometric features such as image points, straight lines, shapes and even 3D poses [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF]. Exploiting these features requires complex image processing like their robust extraction, their matching and their real-time spatiotemporal tracking in the images acquired during the servoing.

To avoid these disadvantages, it has been proposed to directly use all image intensities as visual features (luminance feature) [3][4]. This approach known as photometric visual servoing only requires as image processing the computation of the image spatial gradient. Visual servoing based on the luminance feature has shown very accurate positioning even with approximated depths, partial occlusions, specular and low-textured environments. However, the highly nonlinear photometric cost function induces a relatively small convergence domain. Because of that, the visual difference between the initial image and the desired one should not be too large.

Several methods have been proposed to improve the pure photometric approach. In [START_REF] Delabarre | Visual Servoing using the Sum of Conditional Variance[END_REF], the authors propose to adapt the desired image at each iteration of the visual servoing according to the illumination of the current image. The sum of conditional variance between the desired and the adapted current image is computed to achieve direct visual servoing. A robust M-estimator has been directly integrated to the method in order to reject the pixels of the current image which are too different from the desired one [START_REF] Comport | A visual servoing control law that is robust to image outliers[END_REF]. Another approach [START_REF] Dame | Improving mutual information based visual servoing[END_REF] compares the current and the desired images by computing their mutual information. The mutual information is a similarity measure that represents the quantity of visual information shared by two images. Using this similarity measure, the control law is naturally robust to partial occlusions and illumination variations. Another improvement is to represent the images by intensity histograms [START_REF] Bateux | Direct visual servoing based on multiple intensity histograms[END_REF]. The use of these compact global descriptors increases robustness to noise and illumination changes. The control law minimizes the Matusita distance between the intensity histogram of the desired image and the histogram of the current image. However, even with all these improvements, the convergence domain remains relatively narrow.

Recently, photometric moments have been introduced as visual feature [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF]. Pixel intensities of the whole image are considered to compute a set of photometric moments which represent essential characteristics of the image. The structure of the interaction matrix presents nice decoupling properties. This method has shown interesting results even for large displacement. However, the modeling of the interaction matrix is made under a constant image border hypothesis which leads to failures when parts of the scene enter or leave the camera field-of-view. A recent improvement has been proposed to counter this issue [START_REF] Bakthavatchalam | An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance[END_REF]. A spatial weighting function is included into the photometric moments formulation. However, this improvement alters the invariance properties of moments which disturbs the control of the camera rotation around the two axes orthogonal to the camera optical axis. Kernel-based Visual Servoing (KBVS) [START_REF] Kallem | Kernel-based visual servoing[END_REF][12] also enables the independent control of the rotational and the translational motions. The KBVS approach proposes to use kernels as image operators. Gaussian kernels and spatial Fourier transform are used as visual features to control four degrees of freedom (dof) (the three translations and the rotation around the optical axis).

As photometric moments, KBVS methods do not resolve the difficulties regarding the control of the rotations around the two orthogonal axes to the optical camera axis. They also encounter difficulties in dealing with appearance and disappearance of parts of the scene. Wavelets-based visual servoing has also been recently proposed [START_REF] Ourak | Wavelets-based 6 dof visual servoing[END_REF]. The authors have developed a visual servoing control law using the decomposition of the images by the multiple resolution wavelet transform (MWT). Shearlet, a natural extension of wavelet, has also been introduced as visual feature for visual servoing [START_REF] Duflot | Toward Ultrasound-based Visual Servoing using Shearlet Coefficients[END_REF]. The interaction matrix has been developed for shearlet transform coefficients but its computation is numerical and estimated offline. Similarly to wavelets, shearlets are very efficient in denoising and to represent anisotropic images features. Wavelets and shearlets have demonstrated interesting results in terms of robustness, stability and accuracy for six dof. However the convergence domain remains relatively tight for both methods.

In this paper, we propose a new approach that keeps the same advantages as the pure photometric method and considerably increases the convergence domain. For this purpose, each image pixel is represented as a Gaussian function w.r.t. its luminance information. This can be seen as an enlargement of the attraction power of each pixel. The combination of every pixel Gaussian function forms the Photometric Gaussian Mixture that represents the image.

Gaussian mixtures as visual servoing features have also been studied in [START_REF] Hafez | Visual servoing based on gaussian mixture models[END_REF]. In the latter work, image point features are modeled as a Gaussian mixture. Thanks to the representation of the geometric features as Gaussian distributions the points matching and tracking are not necessary. However, the points detection is still a crucial step and points extracted from the images during the servoing must be exactly the same as the points extracted from the desired image. In other words, the desired points should always be in the camera field-of-view and the points detector must have a 100% repeatability. Authors chose to set the same variance for each Gaussian function which induces ambiguities. Similarly, but in a different context, Gaussian mixtures have also been exploited in [START_REF] Boughorbel | A new method for the registration of three-dimensional point-sets: The gaussian fields framework[END_REF] for automatic 3D point clouds registration. Usually, 3D registration are realized in two steps: first, a coarse initialization using feature extraction and matching methods and then a refinement using Iterative Closest Point (ICP) based algorithms. The size of the region of convergence is significantly extended by the use of the Gaussian mixtures and the usual first step of initialization can be avoided. Moreover, this method handles local convergence problems of standard ICP algorithms.

The goal of our visual servoing is to minimize the difference between the desired Gaussian mixture and the Gaussian mixture computed from the current image varying over time. The interaction matrix is a key concept of visual servoing. It expresses the relationship between the image variations and the camera displacements. The photometric Gaussian mixtures of the images have an analytical formulation, thus enabling the determination of the analytical form of this interaction matrix. Our method does not need feature detection and has been extensively evaluated in experiments with a real camera on an actual industrial robot.

While we already published this idea in [START_REF] Crombez | Photometric gaussian mixtures based visual servoing[END_REF], two important contributions are proposed in this paper: a new model of photometric Gaussian has been designed (Section II-A) and a new modeling of the interaction matrix using the Green's theorem has been developed (Section III-A). Comparisons between the two photometric Gaussian models (Section II-C) and an evaluation of their robustness to noise (Section V-C) have been conducted. Both highlight the benefits of the new model. In addition, we present a strategy to initialize and to control the Gaussian extension (Section III-C). New simulations in complex virtual 3D scenes (Section V-B), comparisons with state-of-the-art approaches (Section V-D) and experimentations in real conditions (Section V-E) have also been carried out.

The remainder of this paper is organized as follows. First, Section II presents the formulation of Gaussian mixtures of 2D images. The influence of the Gaussian parameter is also studied. After that, Section III describes the development of the interaction matrix for the Photometric Gaussian Mixture feature. Section IV explains the control scheme to minimize our cost function. Then Section V shows simulation results, comparisons and real experiments. Finally, conclusions and future works are presented in Section VI.

II. PHOTOMETRIC GAUSSIAN MIXTURE AS VISUAL

FEATURE

A digital image described in a 2D discrete space is derived from an analog signal in a 2D continuous space by a sampling process. This digitization can be modelled as the convolution of the continuous signal by a Dirac comb. Therefore, the obtained digital image can be seen as a comb of pulses whose pulses height depends on the pixels intensity. As mentioned above, considering that image pixels have a power of attraction, pixels represented as Dirac pulses (Fig. 1) have a zero attraction everywhere except to their own position. That is why, instead of considering a pixel of an image as a pulse we propose to represent it as a 2D Gaussian. Figure 1 illustrates, for an one-dimensional case, one of the advantages of this new signal representation: the difference between the desired and the current signal varies continuously, while it remains constant for Dirac pulses. This power of attraction enlarges the convergence domain of the visual servoing.

A. Photometric Gaussian Mixture

2D Gaussian is a suitable characterization of the power of attraction that we want to assign to pixels. Indeed, the further away from a pixel location the smaller the attraction of this pixel. A Gaussian function is defined on [-∞, +∞], therefore a pixel has an influence on every other image pixel. Moreover, a Gaussian function is differentiable which is substantial to express analytically the interaction matrix (Section III). Considering two uncorrelated variables x ∈ R, y ∈ R which compose the vector x = (x, y), the two-dimensional Gaussian function is the distribution function given by:

f (x, y) = A exp - (x -x) 2 2σ 2 x + (y -y) 2 2σ 2 y ( 1
)
where A is the amplitude coefficient, x is the expected value and σ = (σ x , σ y ) is the variance. In the following, we do not use this terminology because our use of the Gaussian is not statistical. To avoid ambiguities, we name x as the Gaussian center and σ as the Gaussian extensions along x and y axes. An image is composed by N × M pixels, each pixel has a location u = (u, v) and a luminance I(u). We note I(r) an image acquired from a camera pose r represented by a vector (t x , t y , t z , θ x , θ y , θ z ). I(r) is the stacking of every I(u) where u belongs to the discrete set of the acquired image coordinates grid U. a) Gaussian center x: we want the highest attraction of a pixel to be located at its own position. We therefore consider each image pixel u as a Gaussian function centred around its location. Thus, we have: x = u. b) Gaussian amplitude A: in order to keep pixel distinctiveness, under the temporal luminance constancy hypothesis, pixels must attract the other pixels that have the same intensity, Gaussians amplitude are thus related to pixels intensity:

A = I(u).
c) Gaussian extension σ: considering the Gaussian extension as a power of attraction, there is no reason to give more priority to a pixel than another or to give more priority to an image axis than the others. Thus, every pixel has an equal extension λ g along the u and v image axes:

σ u = σ v = λ g .
Following this parametrization, the photometric Gaussian function g is:

g(I, u g , u, λ g ) = I(u)E(u g -u) (2)
where u g are the Gaussian function coordinates and

E(u g -u) = exp - (u g -u) 2 + (v g -v) 2 2λ 2 g , (3)
for compactness. A Gaussian mixture is defined as the sum combination of a finite number of Gaussian density functions. We denote the Photometric Gaussian Mixture associated with an image I computed with an extension parameter λ g by G(I, λ g ) (Fig. 2). More precisely, a spatial sampling of the Photometric Gaussian Mixture associated with the image I at u g is expressed as:

G(I, u g , λ g ) = u∈U g(I, u g , u, λ g ) = u∈U I(u)E(u g -u) (4)
Fig. 2: A Photometric Gaussian Mixture G(I) (red) of a onedimensional image composed by 3 pixels I(u 0 ), I(u 1 ) and I(u 2 ) (respectively the blue, green and yellow pulses). G(I) is the sum combination of every Gaussian function related to image pixels.

Note that another interpretation of G(I, u g , λ g ) consists to consider the intensity of every pixel u around u g with a weight proportional to the distance between u and u g . G is the stacking of every G. Fig. 3 shows several Gaussian mixtures G(I, λ g ) of a same image I for different extension parameters λ g .

B. Role of the Gaussian extension

Visual servoing based on the pure photometric feature (pixels luminance of the entire images) is designed to minimize the difference between a desired image I(r * ) and images acquired during the servoing I(r). Vectors r * and r respectively represent the desired and the current camera poses. Because of the limited convergence domain, the displacement between r * and r must be short enough to guarantee a large photometric overlapping between the initial image and the desired one.

With a large extension parameter λ g , the power of attraction of every pixel is more important. Thus, the Gaussian mixture representation (e.g. Fig. 3f) offers more chance to have overlapping areas between the representation of the desired image and the representation of the current image. On the contrary, it is interesting to note that for a small extension parameter, the Gaussian mixture (e.g. Fig. 3b) and the original image (Fig. 3a) are similar. This can be verified by:

lim λg-→0 G(I, u g , λ g ) = I(u)| u=ug (5)
and consequently, The comparison of cost functions shape obtained for different λ g confirms these observations. Let us consider that Fig. 3a shows the desired image. Fig. 4a shows the cost function obtained by computing the differences between the desired image I(r * ) and images I(r) obtained around the desired pose w.r.t. translations along u and v image axes. Fig. 4(bf) show cost functions obtained by computing the differences between the desired Gaussian mixture G(I(r * ), λ g ) and the Gaussian mixture G(I(r), λ g ) for different extension parameter λ g values. We can see that the higher λ g is, the larger the convex domain. On the contrary, when λ g is low, the cost function (Fig. 4b) has the same tight shape than the photometric one (Fig. 4a). It appears thus appropriate to start the visual servoing with a high λ g to ensure the convergence and then decrease it for obtaining a high accuracy. The control of this supplementary degree of freedom simultaneously to the classical ones (i.e. the six dof of the camera) is discussed in Section III-C.

lim λg-→0 G(I, λ g ) = I. ( 6
) (a) (b) (c) (d) (e) (f)

C. Photometric Gaussian models comparison

In [START_REF] Crombez | Photometric gaussian mixtures based visual servoing[END_REF], we had proposed the following photometric Gaussian expression (Model 1):

g(I, u g , u, λ g ) = exp - (u g -u) 2 + (v g -v i ) 2 2I 2 (u)λ 2 g ( 7
)
where the extension parameter is related to the pixel intensity:

σ u = σ v = λ g I(u)
and where every pixel has the same Gaussian amplitude: A = 1.

Due to the direct relation between the extension parameter and the pixels intensity, Model 1 is less robust to image noise than the new proposed model. Indeed, noise may generate low and high intensity pixels distributed over the images that affect the shape of the Gaussian mixture. To highlight that problem, Gaussian mixtures have been computed from one-dimensional images for ease of reading. The first row of Fig. 5 shows three one-dimensional images. The first one is noiseless, and noisy pulses that follow normal distributions with different variances have been added on the two others. The second and the third rows show the Gaussian mixtures computed using respectively Model 1 [START_REF] Crombez | Photometric gaussian mixtures based visual servoing[END_REF] and the new model presented in this paper (Section II-A). Se can see that a small noise has a significant influence on the shape of the Gaussian mixture when the Model 1 is considered. The noisy pixels produce narrow Gaussians and the resulting Gaussian mixtures are less convex and do not provide the expected representation of the initial images. On the contrary, the new model ( 2) always produces a smooth Gaussian mixture.

These dissimilarities between the two photometric Gaussian models induce different photometric Gaussian mixtures and thus different visual servoing behaviors. Experimental results (Section V-C) validate these observations and illustrate the higher performances of the new model.

III. MODELING OF THE INTERACTION MATRIX

The key of visual servoing is the interaction matrix L s which links the time variation of visual features s to the camera velocity v = (υ, ω) where υ and ω are respectively the linear and angular components [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF]:

ṡ = L s v (8)
In order to use Photometric Gaussian Mixtures as visual features, we have to develop the interaction matrix for s = G (Eq. 4). The Photometric Gaussian Mixture based Visual Servoing objective is to regulate to 0 the following error:

= G(I(r), λ g ) -G(I(r * ), λ g * ) = G -G * (9)
where G(I(r * ), λ g * ) is the Gaussian mixture associated with the desired image computed with a fixed extension parameter λ g * and G(I(r), λ g ) is the Gaussian mixture associated with the image acquired during the servoing computed with the extension parameter λ g .

We present in this paper two methods to perform the modeling of L G . As in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], the first approach is based on the Green's theorem which permits to avoid the computation of image gradients. The second is based on the assumption that the 3D scene itself is a Gaussian mixture viewed by a camera. Comparison results (Section V-A) show that this assumption offers good approximations of the interaction matrix.

A. Method #1: Green's Theorem Modeling (GTM)

In ( 4), the only term that is time dependent is the image intensity I(u) that varies due to the camera displacements. Consequently:

Ġ(u g ) = u İ(u)E(u g -u) , (10)
Considering, as usual, that the temporal luminance constancy hypothesis is ensured [START_REF] Collewet | Visual servoing set free from image processing[END_REF][9] [START_REF] Bakthavatchalam | An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance[END_REF], which means the optical flow constraint (OFC) [START_REF] Horn | Determining optical flow[END_REF] is valid, we have [START_REF] Collewet | Visual servoing set free from image processing[END_REF]:

İ(u) = - ∂I(u) ∂u ∂u ∂r ṙ = -∇I T L u v (11)
where L u is the interaction matrix related to an image pixel point. Knowing the intrinsic parameters of the perspective projection model: α u , α v for the horizontal and vertical scale factors of the camera photosensitive matrix and u 0 , v 0 for the principal point coordinates, L u can be decomposed as:

L u = L ux L x = α u 0 0 α v -1 Z 0 x Z xy -(1 + x 2 ) y 0 -1 Z y Z 1 + y 2 -xy -x (12)
where L x is the so called interaction matrix related to a point x = (x, y), expressed in the normalized image plane of the perspective projection model [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF].

Injecting [START_REF] Kallem | Kernel-based visual servoing[END_REF] in [START_REF] Bakthavatchalam | An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance[END_REF], we obtain:

Ġ(u g ) = - u ∇I T E(u g -u)L u v ( 13
)
from which we deduce by analogy with ( 8):

L G = - u ∇I T E(u g -u)L u = [L Gvx L Gvy L Gvz L Gwx L Gwy L Gwz ] (14)
Each component of the interaction matrix L G involves image gradients ∇I = [ δI δu , δI δv ] T . With regards to the pure photometric feature [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF], these image gradients are an approximation of the actual image derivatives computed using image processing methods (derivative filters) along horizontal and vertical image axes. Concerning the photometric moments [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], authors proposed to avoid the image gradients computation using the Green's theorem. This simplification step can also be employed for the Photometric Gaussian Mixture feature.

For example, let us consider L Gvx the component of L G corresponding to the camera translational velocity in x. We suppose that the scene is planar and that the depth of all the points is equal to a constant value Z 1 :

L Gvx = - u ∂I ∂u E ∂I ∂v E α u 0 0 α v -1 Z 0 = α u Z u ∂I ∂u E (15)
We note Q = EI and P = 0, then:

∂Q ∂u = ∂E ∂u I + E ∂I ∂u , ∂P ∂v = 0 (16)
The Green's theorem gives us:

u ∂Q ∂u - ∂P ∂v = ∂u P + ∂v Q, (17)
1 One can note that, as for the photometric moments [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF], it is also possible to model the interaction matrix considering that

1 Z = Ax + By + C. This part is left for future work.
and given that ∂P ∂v = 0 from equation ( 16):

u ∂Q ∂u = ∂v Q (18)
Substituting ( 16) in [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF], we obtain:

u ∂E ∂u I + u E ∂I ∂u = ∂v EI (19)
Therefore, we have:

u ∂I ∂u E = ∂v EI - u ∂E ∂u I (20)
As in [START_REF] Bakthavatchalam | Photometric moments: New promising candidates for visual servoing[END_REF] we consider being under the zero border assumption. Indeed if the pixel intensities I(u) lying on the border of the image are all zero, the term ∂v EI in (20) is equal to zero. Then, replacing (20) in [START_REF] Hafez | Visual servoing based on gaussian mixture models[END_REF], we finally get:

L Gvx = - α u Z u I ∂E ∂u (21)
where ∂E ∂u is easy to express from (3):

∇E T = δE δu δE δv =   (ug-u) λ 2 g E (vg-v) λ 2 g E   (22)
The interaction matrix of each component can be developed following the same way. We obtain:

L Gvx = - α u Z u ∂E ∂u I L Gvy = - α v Z u ∂E ∂v I L Gvz = 1 Z α u u ∂E ∂u Ix + α v u ∂E ∂v Iy + 2 u EI L Gwx =α u u ∂E ∂u Ixy + α v u I ∂E ∂v (1 + y 2 ) + 3 u EIy L Gwy = -α u u ∂E ∂u I(1 + x 2 ) -α v u ∂E ∂v Ixy -3 u EIx L Gwz =α u u ∂E ∂u Iy -α v u ∂E ∂v Ix (23)
It is important to note that the image gradients are not needed anymore to compute the interaction matrix. We propose in the next subsection a second method to model the interaction matrix that is easier and faster to compute.

B. Method #2: "Photometric Gaussian Consistency" (PGC)

Photometric Gaussian Mixtures are representations of images -images which are the results of the projection of a 3D scene on a 2D plane. To model the interaction matrix, we consider here that the Gaussian mixtures are direct representations of the 3D scene itself. In other words, we consider that the scene is already a Gaussian mixture. This assumption can be seen as an approximation of the projection of the 3D Gaussians on the 2D Gaussian mixtures. Thus, the temporal luminance consistency is adapted to the temporal photometric Gaussian consistency:

g(u g + ∆u g , t + ∆t) = g(u g , t), (24)
where u, I and λ g are omitted in function g (25) for compactness. A first order Taylor development of (24) gives:

∂g ∂u g ug + ∂g ∂t = ∇g T ug + ġ = 0 (25)
with ∇g T the spatial gradient of g(u g , t) and ġ its time derivation, thus leading to a relationship similar to the optical flow constraint [START_REF] Kallem | Kernel-based visual servoing[END_REF]. Assuming this temporal photometric Gaussian consistency and following a similar reasoning, ( 4) is now decomposed as:

Ġ = u ġ = - u ∇g T ug = - u ∂(I(u)E(u g -u)) ∂u g ug = - u I(u) ∂E(u g -u) ∂u g L ug v (26)
Here, L ug is out of the summation because it does not depend on u while, in the Green's theorem based modeling, L u is inside the summation [START_REF] Duflot | Toward Ultrasound-based Visual Servoing using Shearlet Coefficients[END_REF]. The computation is consequently faster with (26). As in the previous modeling (Section III-A), the expression of each component of the interaction matrix L G does not contain image gradients but Gaussian derivatives ∇E g :

∇E g T = δE δug δE δvg =   - (ug-u) λ 2 g E - (vg-v) λ 2 g E   = -∇E T (27)
Then, an analytic formulation of the interaction matrix of each component can be developed:

L Gvx =α u K u /Z L Gvy =α v K v /Z L Gvz = -(α u K u x g + α v K v y g ) /Z L Gwx = -α u K u x g y g -α v K v (1 + y 2 g ) L Gwy =α u K u (1 + x 2 g ) + α v K v (x g y g ) L Gwz = -α u K u y g + α v K v x g ( 28
)
where

K u = u I ∂E ∂ug and K v = u I ∂E ∂vg .
It is interesting to note that even if the two modeling approaches are different, the expressions of the interaction matrix are close (see (23)). More precisely, the components L Gvx and L Gvy are strictly identical since ∇E T = -∇E g T while the other ones are more simple. Furthermore K u and K v are computed once for every u g while in the Green's method a lot of terms have to be computed.

C. Extension parameter modeling

The influences of the extension parameter observed in Section II-B show that it is interesting to start the servoing using a high extension parameter and to finish it with a small one. Indeed, this respectively ensures to enlarge the convergence domain while keeping the convergence accuracy at least similar to the pure photometric case. The extension parameter λ g is optimized during the visual servoing as well as the camera velocities v. To this end, we compute the derivative of the Photometric Gaussian Mixture with respect to the parameter λ g :

∇Λ G(ug) = u ∂g(I, u g , u, λ g ) ∂λ g (29) = u I(u) ∂E(u g , u) ∂λ g (30) = u I(u) (u g -u) 2 λ 3 g exp - (u g -u) 2 2λ 2 g ( 31
)
IV. IMPLEMENTATION AND VALIDATION

We extend the classical control law [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF] to:

v λ = -µ L + G λ (G -G * ) (32)
where v λ = (υ, ω, δλ g ) is composed of the three linear and three angular camera velocities, and the extension parameter increment. The convergence can be improved by tuning the gain µ. A Gauss-Newton scheme is used here to validate the use of the Photometric Gaussian Mixture feature.

Other optimization algorithms like the Efficient Second order Minimization [START_REF] Malis | Improving vision-based control using efficient second-order minimization techniques[END_REF] or the Levenberg-Marquardt [START_REF] Hammouda | Real-time visual servoing based on new global visual features[END_REF] could be used. L + G λ is the pseudo inverse of the extended interaction matrix L G λ related to the Gaussian mixture G. The stacking of every L G interaction matrices (14) related to Gaussian mixture values at every location u g and the equation (31), gives the global interaction matrix:

L G λ =     . . . L G(ug) Λ G(ug) . . .     , (33)
where L G(ug) is computed with (23) or (28), depending on the chosen model, and Λ G(ug) from (31). For both methods, the computation of L G λ requires the scene depth Z from the camera [START_REF] Swensen | Empirical characterization of convergence properties for kernel-based visual servoing[END_REF]. In this paper, we consider the depth of each point as constant and equal to the distance between the scene and the camera at the desired pose ( Z = Z * ). That is why the pseudoinverse of the interaction matrix is noted L + G λ . This control law is locally asymptotically stable when L + G λ L G λ > 0. The choice of initial and desired λ g must ensure a large convergence and a high accuracy. These two aspects are respectively satisfied by initializing λ g with a large value and by ending the servoing with a very small value. We propose a two-steps extension parameter evolution strategy:

• Step 1 to enlarge the convergence domain: we initialize λ g with a large value λ gi and the desired extension parameter is initialized following λ g * = λ gi /2.

• Step 2 to accurately finish the visual servoing: both λ g and λ g * are set with a same small value.

The transition from the first to the second step is performed when |λ g -λ g * | < S where S is an empirically fixed value.

More details about the choice of the latter parameters values are provided in Section V-F.

V. EXPERIMENTAL RESULTS

The two proposed interaction matrix modeling methods, i.e. the Green Theorem Modeling (GTM) and the Photometric Gaussian Consistency (PGC), are compared in Section V-A. Then, challenging experiments that highlight the contributions of using Photometric Gaussian Mixtures as visual features are presented in Section V-B. Section V-C studies the robustness of the proposed method when some noise is added to the images. In Section V-D, our method is opposed with three state-of-the-art approaches. Section V-E presents experiments conducted with a real robot in 3D environments. Finally, the initialization of the required parameters are discussed in Section V-F.

A. Comparisons between GTM and PGC

Experiment #V1 (Fig. 6): The first simulated experiment demonstrates the concept of the proposed approaches for 2 dof (v x and v y ). A single untextured object is present in the camera field-of-view. At the desired camera pose, the object is projected on the top right corner of the image (Fig. 6d). At the initial camera pose, the object is projected on the bottom left corner of the acquired image (Fig. 6b). The projections of the object on the desired and initial images have not any photometric overlapping area. The initial extension parameter has been chosen such that λ g is high enough to induce a sufficient overlapping area between the two Gaussian mixtures. In the first step of the extension parameter evolution strategy, the desired Gaussian mixture is computed using a fixed λ g * = 0.5 * λ gi where λ gi = 90.0. The transition to the second step is performed when |λ gi -λ g * | < 0.1. Both extension parameters λ g * and λ g are set to 1.0 (as explained previously in Section IV).

As seen on Fig. 6, the two modeling methods (GTM and PGC) give strictly the same results because the components L Gvx and L Gvy are strictly identical. In the beginning, the error does not decrease exponentially because of the large error between the desired and the initial images. But after a few iterations, it decreases exponentially (Fig. 6h and Fig. 6l). Despite the large displacement between the desired and the initial poses, the camera has converged to the desired one as well as the Gaussian expansion parameter (Fig. 6i and6m). This example shows how the Photometric Gaussian Mixture drastically enlarges the convergence domain of the photometric visual servoing while maintaining the final accuracy. In this situation, the pure photometric feature [START_REF] Collewet | Visual servoing set free from image processing[END_REF] could not at all drives the camera to the desired pose.

Experiment #V2 (Fig. 7): This experiment has been conducted using a complex scene (textured plane) containing several shapes and controlling six dof. Note that some parts of the scene which are present in the desired image (Fig. 7a) are not visible in the initial one (Fig. 7b) and vice versa. The image and its associated Gaussian mixture at the end of the servoing (Fig. 7d) illustrate that the camera has reached its desired pose. However, due to the visual appearances and disappearances of Translational velocities (m/s) parts of the scene during the servoing, the residual error does not follow a perfect exponential decrease (Fig. 7h and Fig. 7l).

Iterations v X v Y (n) -5 0 5 -2 0 2 4 6 X(m) Y(m) (o)
The Gaussian expansion parameter has converged to the value Translational velocities (m/s) Rotational velocities (rad/s) Translational velocities (m/s) Rotational velocities (rad/s) used for the desired Gaussian mixture (Fig. 7i and Fig. 7m).

Iterations v X v Y v Z w X w Y w Z (j)
Iterations v X v Y v Z w X w Y w Z (n) -3 -2 -1 0 0 1 2 3
The initial and final difference images (Fig. 7e and Fig. 7g) highlight the large gap at the beginning of the task and the accuracy of the convergence at its end. Even if the GTM modeling is mathematically more rigorous than the PGC one, the obtained visual servoing behaviors are very close. PGC can actually be seen as a very good approximation of GTM. We use the PGC modeling for the following experimentations because of its smaller processing time.

B. Visual servoing in complex and 3D virtual environments

Experiment #V3 (Fig. 8): In addition to a complex textured plane, we consider a very challenging positioning task. As can be seen on the initial image of differences (Fig. 8b and8d), the desired and the initial images do not share a lot of overlapped photometric information. In addition, the textured plane is partially outside the camera field-of-view in the initial image (Fig. 8). The displacements on (t x , t y , t z , θ x , θ y , θ z ) 2 is given by (10.31m, -4.39m, 0.1m, -4.58 • , 0 • , -13.75 • ). Despite this very large displacement and the high difference between the initial and the desired images, the Photometric Gaussian Mixture-based Visual Servoing has successfully controlled the camera motion to precisely reach the desired pose. The difference image at convergence is null. Of course, because of the very small initial overlapped photometric information, the trajectory taken by the camera to converge to the desired pose (Fig. 8k) is not straight.

Experiment #V4 (Fig. 9): In this simulation, a 3D scene is considered. We suppose that the depth of the scene is unknown. We used the same value ( Z = 50m) as approximation for every pixel. In simulation, this depth is available and could be used in Equations ( 28) and ( 32), but we choose to "ignore" it in order to have relevant comparison with real scenes for which the depth is unknown. The displacement between the desired and the initial camera poses is (26.95m, -5.02m, -11.14m, 14.40 • , -27.54 • , -5.88 • ). We can observe in Fig. 9e that image differences are also very large. Fig. 9c corresponds to the iteration just before the step transition. The camera converges perfectly to the desired pose (Fig. 9g) with a final pose error equal to (1.6mm, 2.6mm, 4.1mm, 0.02 • , 0.01 • , 0.02 • ). We insist that the transformation between the initial and desired poses, and particularly the orientations around the two orthogonal axes to the optical camera axis, are very important and make this case very challenging. The visual differences between the initial and the desired images also reflect the challenging nature of this experiment.

C. Evaluation of the robustness to noise

Fig. 10 and Table I describe an experiment that evaluates the noise robustness of the proposed approach where the two models (( 2) and ( 7)) of photometric Gaussian (Section II-C) 2 In the following, all the differences between desired and initial image poses are given in this order: (tx, ty, tz, θx, θy, θz), with tx, ty and tz are in meters, θx, θy and θz are in degrees while positioning errors at convergence are given in millimeters. are compared. The displacement between the desired and the initial camera poses is always the same for each experiment: (6.94m, 7.71m, 0.76m, 41.56 • , 33.19 • , 60.34 • ). First, visual servoing is carried out without any image noise (Fig. 10a). Then, a Gaussian noise is added on both the desired and the current images (Fig. 10b-e). More precisely, the noise added on the desired image is static and the noise added on the current images varies at each iteration. Between every experiment, the intensity of the noise is enhanced increasing its standard deviation σ.

Iterations v X v Y v Z w X w Y w Z (j)
Table I shows the final errors at convergence for the photometric Gaussian model proposed in [START_REF] Crombez | Photometric gaussian mixtures based visual servoing[END_REF] and the new model proposed in this paper. As we can see, photometric Gaussian mixtures based visual servoing is particularly robust against image noise. As expected from the theoretical expressions (Section II-C), the proposed new model of photometric TABLE I: Noise robustness evaluation: final error (position and orientation) of the estimated pose regarding the used photometric Gaussian models and the standard deviation σ of the Gaussian noise (Fig. 10).

Iterations v X v Y v Z w X w Y w Z (j)

New model

Model Gaussian is more robust to noise than Model 1. Of course, the accuracy at convergence decreases as the noise intensity increases but it remains very good regarding the excessively high noises that are added.

D. Comparisons with state-of-the-art methods

In this section we oppose our method (PGC) to two state-ofthe-art methods: the pure luminance (PL) [START_REF] Collewet | Photometric visual servoing[END_REF] and the weighted photometric moments (WPM) [START_REF] Bakthavatchalam | An Improved Modelling Scheme for Photometric Moments with Inclusion of Spatial Weights for Visual Servoing with Partial Appearance/Disappearance[END_REF]. The PL method is based on a Levenberg-Marquart optimization scheme, while the two other methods use Gauss-Newton algorithm. Simulations are conducted to ensure that the conditions of experimentation are exactly the same. Fig. 11 shows the obtained camera 3D trajectories using the three methods in the case of a small translation facing a planar scene. We can observe that the WPM approaches gives a quasi-straight 3D trajectory, while those obtained with the PL and the PGC methods are more twisted.

The convergence efficiency of the PGC (Model 1 and the new one), PL and WPM methods are then compared for the difficult case of a 3D scene. The goal of this experiment is to reach a same desired pose starting from 20 random initial camera poses. Fig. 12a shows the desired image and Fig. 12 shows the initial images generated from these 20 random poses. The size of the images is 200 × 150 pixels for every method. The two photometric Gaussian models for the PGC approach are also compared for several initial extension parameters λ gi .

We consider that a camera has successfully converged to the desired pose when the final error is less than (5.00mm, 5.00mm, 5.00mm, 0.1 II provides the successful and failed convergences of this experiment. Note that each method has the same parametrization throughout the experiments. The new model of photometric Gaussian is slightly better than the Model 1 proposed in [START_REF] Crombez | Photometric gaussian mixtures based visual servoing[END_REF]. When the extension parameter of the new model is initialized with a high enough value, the PGC approach converges for almost every initial pose (and even all the 20 poses when λ gi = 25). Even if the results obtained with the Model 1 are good too, it is more difficult to identify a unique initial value of λ gi that works for every initial pose. The PL approach has only drove the camera to the desired pose in 3 out of the 20 initial poses and numerous iterations were required. The WPM method diverges completely for 11 poses, which is not surprising due to the large discrepancies between the initial and desired images. For the 9 poses marked with an orange check mark in Table II, it drives the camera next to the desired pose with a non negligible final visual alignment. This is probably due to the fact that a 3D scene is considered, which implies that parts of the scene appear and disappear during the camera motion, inducing ambiguities in the values of the photometric moments involved. However, the final visual alignments are low enough to complete the convergence to the desired pose by switching for instance to the PL method.

• , 0.1 • , 0.1 • ). Table
Finally, Fig. 13 reports some challenging cases with different scenes where large parts of the desired image are absent in the initial image. In these situations, only the PGC method drives successfully the camera to the desired poses. For all these cases, the extension parameter has been initialized with the same value λ gi = 5.

To conclude this comparison, the PL is very accurate, fast and does not need any parameter tuning, but the convergence domain is limited. The WPM enlarges the convergence domain, it is fast and guarantees a more direct 3D trajectory, when it converges but it is inefficient in case of 3D scenes when parts appear or disappear near the center of the image. Our PGC method enlarges drastically the convergence domain and is accurate, but it is more time consuming. As a comparison, considering 100 × 100 pixels images, one iteration of the PL method is done in less than 20ms. One iteration of our parallelized version of the PGC method takes 100ms on a Intel Core I7 2.3GHz with a GeForce GT 630M GPU running Linux. This allows a frequency for the servo loop around 10Hz. II.

E. Real experiments

Three real experiments using a 6 axis industrial robot (Stäubli TX60) with a perspective camera mounted on its endeffector are presented. Figs. 14, 16 and 18 show this robot and the experimental environments. The intrinsic parameters of the camera have been estimated. The depth Z is unknown and is supposed constant for every pixel all along the motion of the camera.

Experiment #R5 (Fig. 15): The goal of this experiment is to demonstrate that the proposed visual servoing works well even under common lighting conditions with large TABLE II: PGC, PL and WPM methods comparion. Successful () and failed () convergences for 20 random camera initial poses (Fig. 12). An orange mark (H) means that the camera has converged next to the desired pose with a non negligible final visual alignment. For the two PGC approaches, several initial extension parameters have been used. differences in the images. As it can be seen in Fig. 14, the scene contains several 3D objects. The mean distance between the scene and the camera is about 0.5 meter at the desired pose. We used this distance as depth Ẑ for every pixel to compute the interaction matrix. This scene is quasi-planar with a relative depth of approximatively 0.02m. The initial displacement is composed by translations and rotations (-0.043m, -0.300m, 0.018m, 20.80 • , -7.84 • , 5.97 • ). Several occlusions have been voluntary introduced during the visual servoing (Figs. 15h-15k). These disturbances affect the behavior of the control as it can be seen on the residual error, velocities and extension parameters curves around the iterations 500 and 1100. At convergence, the final image of differences (Fig. 15g) is not absolutely null due to a global illumination change between the beginning and the end of the experiment, but the final pose error is very small (1.48mm, 0.97mm, 0.54mm, -0.13 • , -0.21 • , 0.06 • ).

Experiment #R6 (Fig. 16 and Fig. 17): As previously, this experiment has been conducted under common lighting conditions with a large displacement between the initial and the desired camera poses (0.0m, 0.265m, 0.040m, -0.89 • , 0.05 • , 27.04 • ). The visual difference between the initial and desired images is also important. The 3D scene contains various shapes of objects (monitors, robots, specular surfaces, ...) and different colors as it can be seen in Fig. 16. The relative depth of the 3D objects present in the scene is around one meter. Moreover, scene occlusions and lighting changes have been introduced during the experiment between the iterations 50 and 150. For example, the Fig. 17d shows the intrusion of an object in the camera field-of-view. This explains why the curves (Figs. 17e, 17f and 17g) are shaky. Because of our robot singularities and its limited working space, the initial displacement in terms of image differences is not as impressive as for the experiments conducted in virtual environments. However, this experiment shows that the visual servoing based on photometric Gaussian Mixtures works well even in real conditions. The final error at convergence is (-1.39mm, 6.67mm, -1.25mm, -0.9 • , 0.05 • , 0.95 • ). This relative accuracy is due to variations in light conditions.

Experiment #R7 (Fig. 18 and Fig. 19): The last experiment is achieved in a real environment (windows, tables, mobile robots, ...) with a large relative depth as it can be observed in Fig. 18. The displacement between the initial and the desired poses is (0.194m, 0.0584m, 0.062m, 0 • , 0 • , 14.10 • ).

(c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o)
However, as usual we used the same depth Ẑ = 2m for every pixel. As before, scene occlusions and lighting changes have been introduced during the process. The two columns of Fig. 19c show respectively external views of the environment and the image acquired by the robot camera at four moments of the visual servoing process. More precisely, the first row corresponds to the start, the second shows an obstruction of the light source, the third highlights the intrusion of a person in the camera field-of-view, occluding the top-right part of the scene, and the fourth row corresponds to the desired state. The introduced perturbations explain why the residual error curve (Fig. 19i) is not smooth. Despite these perturbations, the visual servoing converges and remains stable thanks to the redundancy of the used information. The final error is (7.8mm, 3.5mm, 0.5mm, 0.01 • , 0.52 • , 0.39 • ).

The convergence is a bit less accurate than for the previous experiments. This is due to the higher relative depth of the scene and to the natural light coming from the windows (Figure 19).

F. Discussion

The previous results show that our method ensures the convergence even if there is almost no overlapping in the desired and initial images within realistic 3D environments. However, as explained in Section IV two parameters are involved in our method: the initial value λ gi of the extension parameter and the switching threshold. In general, the choice of these parameters depends on the images and on the difference between the desired and initial images. As it is very difficult to establish the exact relation to estimate these values, we propose an empiric strategy:

• λ gi must confer to the Gaussian mixture a huge power of attraction (Fig. 3f). • λ g * must guarantee that the desired Gaussian mixture overlaps the current one. The first aspect can be justified by the relation that is observed between the size of the Gaussian Mixture (close to a unimodal Gaussian for a large λ gi ) and the size of the image. In the presented experiments, V 2, V 3 and V 4 have all the same λ gi because the size of the images (200 × 150 pixels) is the same and these images are textured. In contrary, for experiment V 1, the value of λ gi is large because there is no texture in the image thus the Gaussian centers are far away. This is an unusual case, which explains why λ gi has a high value. For the real experiments, the values of λ gi are lower than those used in virtual scenes because the images size has been reduced to speed up the calculation time for the robot control. The link between the size of the image and λ gi appears in Table III where the reported values of λ gi have permitted the convergence of the visual servoing regarding three sizes of the image (Fig. 20). However, a too large λ gi leads either to divergence or to an imprecise control of the rotations because the Gaussian mixture is close to a unimodal Gaussian. To show the influence of λ gi , we present in Fig. 20, 3D trajectories produced by the visual servoing for an easy case (small difference between desired and initial images). We can observe that the quality of trajectories depends on the chosen value of λ gi , but the visual servoing converges for a large interval of this parameter: λ gi ∈ [2, 20] for a 40 × 40 image size. So, a value between 2 and the half image size can be assigned to λ gi . The second parameter is the switching threshold. It is not useful when the scene is very simple. For example, in Fig. 6, the scene contains only one object and a Gaussian mixture with a large value of λ gi is sufficient to cancel the difference between the desired and the current images. When the scene is complex, the value of the switching threshold S has consequences on the behavior of the visual servoing.

If it is too small, the algorithm switches too late. Then, λ g stays large and the visual servoing may diverge because the orientations are poorly controlled in this case. If this threshold is too high, switching happens too early and the visual servoing does not converge because the current image is too far from the desired image (local minimum). However, for all our different experiments, we set a same threshold equal to 0.1.

VI. CONCLUSION

We introduced in this paper Photometric Gaussian Mixtures as visual features for dense visual servoing. Limitations of the convergence domain regarding the pure photometric feature have encouraged this research. Usual dense visual servoing methods fail if there is not enough shared photometric areas between the desired and the initial images. Our basic idea was to assign a power of attraction to each image pixel. For that, instead of using images as intensity pulses, we consider every pixel as a Gaussian function. The combination of every Gaussian creates a Photometric Gaussian Mixture which is a representation of the image. Thanks to this representation, even if there is almost no overlapping between the desired and the initial images, the convergence domain is sufficiently enlarged to drive the camera to the desired pose. Beyond the power of attraction concept, Photometric Gaussian Mixtures are also tunable. Indeed, in addition to the camera velocities, the extension parameter is also optimized during the servoing. The variation of the Gaussian extension allows us to enlarge the convergence domain and to ensure a convergence as accurate as with the pure photometric feature.

To validate the Photometric Gaussian Mixtures as visual features, a Gauss-Newton control law has been used to minimize the proposed cost-function. The interaction matrix -key of the visual servoing -has been modeled for the proposed Photometric Gaussian Mixtures. We have presented and compared two modeling approaches. As for the photometric moments, the first one uses the Green's theorem to avoid the computation of the image gradients. The second approach is based on 3D assumptions and can be seen as a good approximation of the first modeling. Both simulation and real experiments have been led that confirm the validity of the two modelings, over performing the previous dense visual servoing approaches.

Fig. 1 :

 1 Fig. 1: Image pixels represented as pulses (red) and their associated Gaussian functions (blue)

Fig. 3 :

 3 Fig. 3: Influence of λ g : (a) Grayscale image I, (b-f) Representations of the Gaussian mixture G(I, λ g ) for respectively λ g = {0.1, 5.0, 10.0, 20.0, 25.0}

Fig. 4 :Fig. 5 :

 45 Fig. 4: Cost functions shapes comparison: (a) Photometric dense features, (b-f) Gaussian mixtures for respectively λ g = {0.1, 5.0, 10.0, 20.0, 25.0}.

Fig. 6 :

 6 Fig. 6: Experiment #V1: Comparison between GTM (from h to k) and PGC (from l to o). (a) Desired image and its Gaussian mixture. (b) Initial image and its Gaussian mixture, (c) Image and its Gaussian mixture just before switching, (d) Final image and its Gaussian mixture. (e)-(g) Image of difference, (h) and (l) Residual error, (i) and (m) Extension parameters, (j) and (n) Velocities, (k) and (o) Trajectories.

Fig. 7 :

 7 Fig. 7: Experiment #V2: Comparison between GTM (from h to k) and PGC (from l to o). (a) Desired image and its Gaussian mixture. (b) Initial image and its Gaussian mixture, (c) Image and its Gaussian mixture just before switching, (d) Final image and its Gaussian mixture. (e)-(g) Image of difference, (h) and (l) Residual error, (i) and (m) Extension parameters, (j) and Velocities, (k) and (o) Trajectories.

Fig. 8 :

 8 Fig. 8: Experiment #V3: Complex scene and large difference. (a) Desired image and its Gaussian mixture, (b) Initial image and Gaussian mixture, (c) Image and Gaussian mixture just before switching, (d) Final image and Gaussian mixture. (e)-(g) Image of difference, (h) Residual error, (i) Extension parameters, (j) Velocities, (k) Trajectory.

  velocities (m/s) Rotational velocities (rad/s)

Fig. 9 :

 9 Fig. 9: Experiment #V4: Virtual 3D scene and large difference. (a) Desired image and its Gaussian mixture, (b) Initial image and its Gaussian mixture, (c) Image and its Gaussian mixture just before switching, (d) Final image and its Gaussian mixture. (e)-(g) Image of difference, (h) Residual error, (i) Extension parameters, (j) Velocities, (k) Trajectory.

Fig. 10 :

 10 Fig. 10: Noise robustness evaluation: Desired images (first row) and Initial images (second row) increasing the standard deviation σ of the Gaussian noise. Results of this evaluation are shown in TableI.

  Fig. 10: Noise robustness evaluation: Desired images (first row) and Initial images (second row) increasing the standard deviation σ of the Gaussian noise. Results of this evaluation are shown in TableI.

Fig. 11 :

 11 Fig. 11: (a) Desired image, (b) Initial image and (c) Cameras 3D trajectory obtained using the PL, the WPM and the PGC methods.

Fig. 12 :

 12 Fig. 12: PGC, PL and WPM methods comparion. (a) Desired image and (b) Initial images generated from 20 random camera poses. The results of this experiment are shown in TableII.

Fig. 13 :

 13 Fig.13: Other examples with large difference between initial and desired images, and 3D trajectories for the PGC method.

Fig. 14 :

 14 Fig. 14: Experiment #R5: Experimental environment and the robot (Staubli TX60).

Fig. 15 :

 15 Fig. 15: Experiment #R5: Real 3D scene. (a) Initial image and Gaussian mixture, (b) Image and Gaussian mixture just before switching, (c) Image and Gaussian mixture just after switching, (d) Desired image and its Gaussian mixture, (e)-(g) Image of difference, (h)-(k) Examples of images with occlusions, (l) Residual error, (m) Extension parameters, (n) Velocities, (o) Trajectory.

Fig. 16 :Fig. 17 :

 1617 Fig. 16: Experiment #6: Experimental environment and the view of the scene from the camera.

Fig. 18 :

 18 Fig. 18: Experiment #7: Experimental environment (with distances between some objects and the robot) and the view of the scene from the camera.

Fig. 19 :

 19 Fig. 19: Experiment #R7: Real 3D scene with voluntary disturbances. External views of the environment and the image acquired by the robot camera at four moments of the experiment: (a, b) Initialization, (c, d) High light variation, (e, f) Important occlusion and (g, h) Desired pose. (i) Residual error, (j) 3D Trajectory.

Fig. 20 :

 20 Fig. 20: Influence of the initial extension parameter: (a) desired, (b) initial and (c) difference images ((40 × 40) pixels) and (d) camera trajectories. Legend: SL for straight line between initial and final camera pose; 2, 6, ..., 20 are the values of λ gi .

TABLE III :

 III λ gi interval related to the image size.

	Image size	λ gi
	(40 × 40) ≈ 2 .. 20
	(60 × 60) ≈ 2 .. 30
	(80 × 80) ≈ 2 .. 40


Patrick Laug

Géraldine Pichot

Romain Le Goc

Caroline Darcel

Philippe Davy

Laug

Patrick

Automatic meshing of Discrete Fracture Networks

Keywords: DFN, mesh generation, flow simulation

Introduction

DFN models may contain a large number of fractures, whose sizes span over several order of magnitudes. In addition, fracture intersections also have a wide range of sizes. Their influence is reinforced by the presence, in natural fractured media, of specific configurations like T-terminations that should be properly handled. The rock matrix is considered impervious here. The challenge is to build a mesh of the fractures with necessary quality for subsequent flow simulations. We require that the generated mesh:

• includes the fracture intersections and the highly critical T-terminations, whatever their size and position in the fracture plane, • contains a relevant number of elements to accurately model the flow properties but which remains reasonable for the flow simulation to be feasible.

Presented work

We will present BLSURF FRAC software, which first builds a geometric model including fracture intersections. In a second phase, it generates a mesh of the geometric model by calling an associated planar mesher [START_REF] Borouchaki | Parametric surface meshing using a combined advancing-front -generalized-Delaunay approach[END_REF]. The general scheme is the following : 1. Compute the intersections of the disks and save topological information; 2. Split curves and find common extremities (for intersections and disks boundaries); 3. Build valid curve discretizations; 4. Generate a mesh of each planar disk, using a user-selected planar mesher, and create the final tridimensional mesh.

The first three steps are handled by BLSURF FRAC software. The main difficulty is to build valid curve discretizations (step 3). To do so, BLSURF FRAC implements automatic corrections. As planar meshers (step 4), we propose to consider the two followings, BAMG and BL2D.

As benchmark test cases, we extend those proposed in [START_REF] Laug | Realistic geometric modeling of fracture networks[END_REF][START_REF] De Dreuzy | Synthetic benchmark for modeling flow in 3D fractured media[END_REF] to DFNs generated with the UFM framework [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF][START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF]. They are large scale DFNs where the fracture size distribution matches the observations and where fractures are organized so that large fractures inhibit the smaller ones, creating T-termination configurations.

As illustration, we present the results on such a DFN, called DFNL50. Figure 1 shows this DFN, made of 3477 fractures, together with the meshes obtained with BAMG and BL2D respectively, using a constant mesh step. BL2D also has the abibility to adaptively refine the mesh, for more accurate and less expensive flow computations. As an example, figure 2 illustrates, for the same network, DFNL50, a mesh refinement around the intersections, with a significant reduction of the number of triangles, by comparison with the number of triangles obtained with a constant mesh step generation (as shown on Figure 1). Finally, we will also present the results of steady-state single phase flow solutions on different benchmark test cases. Figure 3 shows the flow solution, for a permeameter test case, on the network DFNL50 meshed with BLSURF FRAC combined with BL2D.

Figure 1 :

 1 Figure 1: (Left) DFNL50 Test case: 3477 fractures, 4334 intersections. Zoom on the mesh generated with: (middle) BAMG: #triangles: 1,513,802; (right) BL2D: #triangles: 1,450,993.

Figure 2 :

 2 Figure 2: DFNL50 Test case -Adaptive mesh refinement with BLSURF -# triangles: 591,117.

Figure 3 :

 3 Figure 3: (Left) Mean hydraulic head computed on the network DFNL50. (Right) Zoom on the conforming mesh generated by BLSURF, #triangles: 5,924,391.


Guillaume Lapouge

Jocelyne Troccaz

Philippe Poignet

Needle localization for needle steering under 3D ultrasound feedback

Introduction

Tools Conclusion References

[1] G. Lapouge, J. Troccaz and P. Poignet, Multi-rate unscented Kalman filtering for pose and curvature estimation in 3D ultrasound-guided needle steering, Control Engineering Practice, vol. 80, pp. 116-124, 2018 [2] N. Hungr, M. Baumann, J.-A. Long, and J. Troccaz, A 3-D ultrasound robotic prostate brachytherapy system with prostate motion tracking, IEEE Transactions on Robotics, vol. 28, no. 6, pp. 1382Robotics, vol. 28, no. 6, pp. -1397Robotics, vol. 28, no. 6, pp. , 2012 [3] R. J. Webster, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, Nonholonomic Modeling of Needle Steering, The International Journal of Robotics Research, vol. 25, no. 5-6, pp. 509-525, 2006 [4]

  H. Younes, S. Voros, and J. Troccaz, Automatic needle localization in 3D ultrasound images for brachytherapy, International Symposium on Biomedical Imaging (ISBI), IEEE, pp.1203-1207, 2018 • PROSPER robot for brachytherapy procedure [2] • 6 degrees of freedom, needle insertion and rotation module Needle localization for needle steering under 3D ultrasound feedback • The tip pose estimation is accurate, robust to needle disappearance and 3D ultrasound imaging artefacts • An adapted segmentation method, inherited from [4], uses the resulting ROI for needle tip segmentation in the 3D US volumes • The good estimation of the tip behavior and its uncertainties could benefit to adapted control laws and path planning methods Robotic device Ultrasound (US) imaging • The needle tip path is modeled by a bicycle kinematic model [3] • β cut is the cutting angle of the needle, as deduced from data fitting 3D kinematic bicycle model Multi-rate unscented Kalman filter Methods • 3D B-mode US volume @ 1Hz • 3D end-fire probe 4DEC-9/10 with Ultrasonix Sonix RP US system Motivation and objectives • In needle steering, estimating the needle pose is a critical problem • In 3D ultrasound volumes, fine needle localization is difficult and requires a combination of estimation and image processing to be successful. Indeed, 3D ultrasound imaging suffers from noise, artifacts and works at a low frequency • We propose a needle tip pose estimation method in the context of 3D robotic needle steering under 3D ultrasound feedback, based on multi-rate, multisensor fusion [1] • This estimation feeds a segmentation algorithm for robust needle detection 3D US imaging in B-mode Pre-operative shear wave imaging (SWE) • SWE estimates tissues Young's modulus from their response to a shear wave US stimulation • The needle curvature is proportional to the tissue stiffness Echogenic Needle Coating Bi-layered Agar phantom and associated SWE The red part is stiffer than the blue one • Needle-related artefacts appear such as comet tail artefacts, reverberation etc. • The needle tip pose and curvature are estimated by a multi-rate unscented Kalman filter • All available measurements are taken into account asynchronously • The filter updates the estimation error covariance online to account for poor 3D US quality and expected changes of tissue stiffness Results • 51 insertions of 8 cm at 1.5 mm.s -1 • Validation on phantoms and ex vivo tissue sample