Arthur Francisco
email: arthur.francisco@univ-poitiers.fr

Noël Brunetière

Gildas Merceron

Damaged Digital Surfaces Also Deserve Realistic Healing

Keywords: inpainting, surface healing, dental micro-wear

published or not. The documents may come

Damaged digital surfaces also deserve realistic healing

Introduction

Image inpainting consists in filling in damaged areas or missing value areas in an image. It is a well-known process in the field of image processing. The aim of image inpainting was originally to restore damaged photographs by removing pits and scratches. The main difficulty in filling in spaces is to obtain a realistic result; ideally speaking the correction should not be apparent. When the damaged areas are small enough, simple tools like Laplace's correction succeed in healing the image. The repair can be compared to a skin stretched over the hole border. Healing a picture in such a way can be satisfactory in low color gradient areas, but as soon as the restored zone crosses a contrasted feature, the result becomes too apparent. It is the key point that has initiated the present study: real worn surfaces -from mechanical components for instanceexhibit most of the time roughness, making Laplace's restoration not fitted for image healing needs; at least, at first sight. It is the only available tool in digital surface software, though.

There are many inpainting methods devoted to picture restoration: some are common-sense processes but other involve modern techniques like AI, artificial intelligence, with impressive results [START_REF] Liu | Image Inpainting for Irregular Holes Using Partial Convolutions[END_REF]. Let's put aside the latter, because they take advantage of complex neural networks that are trained to recreate special features like realistic eyes on a face, windows on a building, etc. By contrast, the need in digital surface restoration is just filling in areas by reproducing the neighboring texture. It is a more affordable task, for which widely used methods like point-and exemplar-based methods are well suited. There is an abundant literature on inpainting technology, presenting the historical needs for that kind of tools in image and video processing. Instead of collecting an exhaustive list of papers on the subject, one may refer to synthetic literature reviews as in Harrison thesis [START_REF] Francis | Image Texture Tools[END_REF], Gong et al. book [START_REF] Gong | Advanced Image and Video Processing Using MATLAB[END_REF], chapter six and Barbu's book [START_REF] Barbu | Novel Diffusion-Based Models for Image Restoration and Interpolation[END_REF] for mathematical-based methods. As mentioned above, we choose to restrict the investigation field to simple algorithms, as the aim of the present work is to prove that affordable alternatives to Laplace's healing are worth to be used for digital surfaces.

Bertalmio et al. [START_REF] Bertalmio | Image Inpainting[END_REF] fill in a hole, smoothly propagating the isophote lines -lines of equal color level -inward, using the surrounding existing points. The outside isophote are determined and extended inside the hole using the discretized Laplacian. Despite the major drawback of not reproducing the image texture -the short wavelengthsthis paper has been a starting point for many other works. A few years later, Bertalmio et al. [START_REF] Bertalmio | Simultaneous structure and texture image inpainting[END_REF] update, and signigicantly improve, the technique using image decomposition: structure and texture. The structure is the macro features, and the texture includes the short wavelength features. The structure restoration is performed following their first work, and the texture creation is ensured according Efros and Leung's seminal paper [?]. Efros and Leung's algorithm works well for textured images, where texture refers to pattern -texels -rather than various short wavelength features. Note that noisy compounds are preserved as well because unlike the preceding approach, the new points are created by copying points of the surroundings. For this purpose, a few known neighboring points of a centered point to be created are compared to the existing ones outside the hole: when a matching set of points is found, the central point is copied. The process goes on until there is no remaining point to be filled. A refined and efficient method is proposed by Harrison [START_REF] Harrison | A non-hierarchical procedure for re-synthesis of complex textures[END_REF] and improved later [START_REF] Francis | Image Texture Tools[END_REF]. Bertalmio's work is based on Partial Differential Equations, Efros' and Harrison's works are based on pixel surroundings, and defined here as point-based. Another way of reproducing information inside a hole, taking into account similar zones outside, is to copy several contiguous pixels at a time. As for point-based techniques, a metric is defined to compare an inside kernel to an outside target kernel; when they match, the inside pixels are replaced by a copy of the outside pixels. Criminisi et al. [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] use this strategy, along with isophote determination, to propagate patterns and preserve high frequencies. The approach, called exemplar-based, also produces very convincing results in removing objects in a picture.

The last two approaches -point-and exemplar-based methods -seem well indicated for a digital surface restoration. Where picture information is made of three color layers, digital surface information is heights. It makes no big difference in the algorithm implementation. However, unlike image inpainting goal, which is only visual, digital surface healing must not significantly alter the subsequent analyzes to carry out. In the present work, we introduce a new approach which aims to preserve the surface anisotropy as well as the height statistics. It is referred as global method. To test the efficiency of the different inpainting methods, three surfaces are especially selected. They are chosen among dental surfaces, because they are sometimes subject to artifacts: in most cases, the artifacts are due to sharp slopes but it can also be due to remaining bubbles in the mold. The surfaces have different anisotropy ratios, and different height statistics. Many parameters are calculated on both original surfaces and arbitrarily damaged then repaired surfaces. If the results remain close for a particular method, it is considered as efficient; the visual result criterion alone is no more sufficient.

Material

As stated before, dental surfaces exhibit a great variability in roughness levels, in feature distribution, etc. To accurately catch the micro-texture of teeth, the dental surfaces are molded with a silicon polyvinylsiloxane elastomer. Then, the mold surfaces are digitized with a confocal, white light profilometer (Leica Microsystem DCM8; 100× magnification lens, Numerical Aperture = 0.9; Working Distance = 0.9 mm, Lateral Resolution up to 140 nm, Vertical Resolution up to 2 nm).

Three surfaces are chosen for their anisotropy differences. The surfaces have their 8 th degree least square polynomial subtracted, so that the long wavelength texture is removed [START_REF] Francisco | Enamel surface topography analysis for diet discrimination. A methodology to enhance and select discriminative parameters[END_REF], [START_REF] Francisco | Gathering and Analyzing Surface Parameters for Diet Identification Purposes[END_REF].

Anisotropy1 surface

Anisotropy1 is a dental surface from a semi-wild cattle Bos taurus from Camargue, Rhône delta, France. Because it feeds on high amounts of herbaceous dicots, a strong anisotropy is observed, see figure 1.

Anisotropy3 surface

Anisotropy3 is a dental surface from a red deer Cervus elaphus from Bia lowieża, Poland which feeds on grass, forbs and fruits. It can be seen on figure 3 that there is no preferred texture direction, thus the surface is globally isotropic.

Methods

First of all, Laplace's restoration is so often used that we shall explain why. The principle is to minimize the surface of the patch given the border heights.

Laplacian filling

When 1D data are to be filled between two points with no particular assumption, linear interpolation is used. However, regarding 2D data inside a closed curve -a "hole" border -linear interpolation is no longer applicable, see figure 4 for an illustration. The lesser constraint becomes the so-called "minimal surface"; to put it simple each point minimizes the surface created by its neighboring points. There is no reason to linearly interpolate [AB] rather than [CD]: if a central point is to be drawn, it is at a fraction of the distance between (AB) and (CD) A mathematical explanation for the use of the Laplacian can be found in the Electronic Supplementary Material (ESM0). Another way of seeing things is to consider that each point inside the hole is the mean of its closest neighboring points. Let Ω be uniformly discretized and f i,j be the surface height at point (i, j) on the grid. For the sake of simplicity, the grid step is supposed to be the same along x and y: h x = h y = h.

The simplest way of expressing that f i,j is the mean of the closest points is:

f i,j = 1 4 (f i,j-1 + f i-1,j + f i,j+1 + f i+1,j )
or,

f i,j-1 + f i-1,j + f i,j+1 + f i+1,j -4f i,j = 0 (1)
Equation 1 is also the discretized Laplacian applied to f . Indeed, discretizing ∂ 2 f /∂x 2 using centered finite differences yields:

∂ 2 f ∂x 2 i,j = 1 h f i+1,j -f i,j h - f i,j -f i-1,j h (2)
and then,

(∆f ) i,j = 1 h 2 (f i,j-1 + f i-1,j + f i,j+1 + f i+1,j -4f i,j )
The most valuable qualities of Laplacian healing are its processing speed and robustness. The later means that no "surprising" data are introduced, whatever the hole shape, because the patch is extremely smooth: it is a soap film lying on the hole boundary, figure 5.

Figure 5: The Laplacian operator fills in Ω from the hole frontier Γ. The shape that is obtained is the same as the one of a soap film lying on the frontier The main drawback of the Laplacian healing is obviously the other side of the minimal surface principle. When the hole is large compared to the apparent wavelength, the patch becomes too much apparent, and may alter the surface parameters, such as the height statistics. Because it does not take into account the hole vicinity, the more complex the surface is, the poorer the result. However, Laplace's interpolation is well suited for small holes, even if they account for an important ratio of missing data, provided that they are well scattered over the surface.

Global approach inpainting

The goal is to propose a novel technique to fill in surface holes taking into account the surface overall anisotropy and the height statistics of the hole vicinity. It is an ad hoc technique that aims to reproduce the global characteristics of the surface inside the hole.

The main idea is to:

Step 1: fill in the hole area using Laplacian interpolation

Step 2: remove the fifth order least square polynomial of the surface masking the hole area

Step 3: determine a convenient digital filter H, and then apply it to the surface heights

Step 4: determine the height profile of the crown and generate as many heights as needed and fill in the hole area

Step 5: reintroduce the polynomial As a result the anisotropy inside the hole area should reflect the overall anisotropy, and the inside heights should exhibit similar statistics as the outside of the hole area. The whole procedure is illustrated in ESM(ESM0).

Step 1: Laplace healing. Step 1 is necessary to compute the surface height AutoCorrelation Function, ACF, because the surface patch is smooth enough to not disrupt the anisotropy characteristics.

Step 2: textural removal. Step 2 must be tuned regarding the surface correlation lengths: the polynomial degree to be removed should be adjusted to the surface wavelengths. A fifth -or less -order polynomial P 5 (x, y) leads to satisfactory results with the surfaces used in the present work.

Step 3: surface reconstruction. The principle developed by Hu and Tonder [START_REF] Hu | Simulation of 3-D random rough surface by 2-D digital filter and fourier analysis[END_REF] is recalled hereafter, with their own notation.

Let η be the initial heights of the surface, z the desired final heights with the right anisotropy and h the digital filter so that z = h * η.

Let us suppose that η is random white noise. Applying a Fast Fourier Transform, FFT, to the preceding relationship, one easily obtains

Z = H × A. Z 2 = H 2 A 2 , Z 2 being the F F T of the ACF of z, Z 2 = F F T (ACF (z)) -Wiener-Khinchin theorem.
A 2 is a constant because ACF (η) = 0 everywhere except in (0, 0). h being symmetric regarding x and y, the following relationship holds:

h = F F T -1 F F T (ACF (z)) (3)
Hence, providing the height autocorrelation -the surface anisotropy -the digital filter h is determined. The strong assumption here is to consider the surface heights as a random distribution, and more specifically as white noise. What are the consequences? Actually, when applying the digital filter h, the resulting heights do not exhibit the exact anisotropy as planned, but an approximated one. However this is not a big deal because the aim of the maneuver is to generate a realistic anisotropy inside the hole. Hence having an approximate one is enough.

Step 4: change the heights inside the hole area. After applying the digital filter to the surface heights, the resulting anisotropy is close to the overall anisotropy but the heights are far from being realistic. The idea as exposed by Francisco and Brunetière [START_REF] Francisco | A hybrid method for fast and efficient rough surface generation[END_REF] is to keep the filtered height order but to replace them by more realistic ones. The later are taken in a crown around the hole; if there are not enough heights then, the crown height profile is used to generate the missing heights as illustrated in figure 6.

Figure 6: The crown ∆Ω heights may not be enough to fill Ω. Therefore the height profile is determined in order to generate additional heights. S = Ω ∪ Ω, hole border:∂Ω

Step 5: reintroduce the macro texture. The last step is the addition of the least square polynomial P 5 (x, y). The surface heights have been altered outside Ω, so the final surface is the modified Ω heights combined with the outside hole Ω initial heights.

Point-based inpainting

For each new point inside Ω, the neighboring existing points are examined to find similar zones in Ω. When a good candidate is found, the point is copied from Ω into Ω. The procedure, based on Efros' work [?] and Harrison's work [START_REF] Harrison | A non-hierarchical procedure for re-synthesis of complex textures[END_REF], is iteratively executed as follows:

Step 1: Initialization

Step 2: Get Ω and ∂Ω, the unfilled points if any, or terminate

Step 3: For a maximum of ten of unfilled points

• find the best zone to search for correspondences in Ω • find the best patch in the search zone and copy the candidate into Ω

Step 4: return to Step 2

Step 1: initialization. Before the procedure starts, the anisotropy ratio is determined. In the present work, it is defined as the ratio Stri of the highest correlation length over the smallest. The correlation length is defined as the length where the ACF is 50% of its maximal value.

If we suppose that what influences a point is its neighbors ±3 points -see the figure 7 -and that the closest points influence more the central point that the far ones, then the weighting function can be chosen Gaussian with a standard deviation equal to 6.

Figure 7: The procedure looks for a similar template in the search zone ∆Ω in the direction of the anisotropy. The template points are not given the same weight: it is a 2D Gaussian function which big axis is the anisotropy direction. Once a similar template is found, the central point height is copied into ∂Ω to fill the hole. The dashed zone is the zone already filled in.

Step 2: get or update Ω. Each unknown point of Ω must be determined. However, as the process is a point-by-point inpainting process, only the border ∂Ω points are needed. Depending on how the ∂Ω point list is determined, the final result can be different. That is why the points are assigned a random order.

Step 3: find a similar template in ∆Ω and fill ∂Ω. If the whole border ∂Ω is filled in one single step, the final result exhibits a kind of connection line as shown on figure 8. Thus, a maximum of ten points are filled in a single step. The surface may exhibit scratches that should be preserved while inpainting Ω. Therefore the first place to look for similarities in ∆Ω is in the direction of the anisotropy. The distance D -in other words, the difference -between a ∂Ω-centered 7 × 7 kernel and a ∆Ω 7 × 7 kernel, is the weighted sum of the height absolute differences. All distances D are calculated, and the ∆Ω template that yields the lowest D is determined. The step ends with the copy of the new height -the template central point -into the unknown point. Step 4: loop to fill in the remaining points.

Exemplar-based inpainting

The procedure is based on Criminisi's work [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], the general principle is the same as the point-based inpainting. The differences are mainly the copy of the whole template found in ∆Ω and the order of filling in ∂Ω. To avoid the patches from being too marked on the final result, leaving a tiled appearance on the repaired hole, the process is carried out four times with smaller and smaller templates.

Step 1: Get Ω and ∂Ω, the unfilled points if any, or terminate. For each point a priority is assigned.

Step 2: For the point of maximum priority, fill in the empty points of the template with the source template points.

Step 3: return to Step 1, three times, with a smaller template width

The only step that is worth detailing is the first one. Features like scratches make the height gradient ∇h locally high, it is therefore a means to detect them. In order to extend the scratches from ∂Ω into Ω in the right direction, ∇∂Ω is determined, and if ∇∂Ω.∇h = 0 at a given point P , it means that the hole frontier is perpendicular to the scratch: P must be firstly filled in. The figure 9 shows a situation where points have already been filled in -dashed zone -and a priority p is assigned to point P : p = ∇∂Ω.∇ ⊥ h. Large values of p on ∂Ω means that the frontier is locally perpendicular to a marked feature.

Figure 9: ∂Ω is not randomly filled. To preserve the texture, the height gradient ∇h is determined on ∂Ω. At point P , if |∇h| is important and if the normal to ∂Ω is colinear to ∇ ⊥ h, it means that P is in the continuity of a relief element, like a scratch, and then is assigned a higher priority.

Surface parameters

The parameters that are tested on the Anisotropy1-2-3 surfaces are detailed by Francisco et al. [START_REF] Francisco | Enamel surface topography analysis for diet discrimination. A methodology to enhance and select discriminative parameters[END_REF], [START_REF] Francisco | Gathering and Analyzing Surface Parameters for Diet Identification Purposes[END_REF]. The parameters can be classified as spatial or anisotropy parameters, height parameters, topology parameters and fractal parameters. Why investigating the effect of surface healing on many different parameters? Provided the fact that the healed surface is to be analyzed, the recreated zones should alter as less as possible the results of the same sane surface: the new patch should meet the anisotropy conditions, the height conditions and the micro-texture conditions of the sane surface.

3.5.1. 3 spatial parameters. Spatial parameters are related to the height ordering. According the way that heights are organized, longer or shorter wavelengths appear. The anisotropy of a surface is related to the presence of preferred texture directions, then, a means to assess the anisotropy is to determined the AutoCorrelation Function, ACF. A detailed presentation is given by Francisco et al. [START_REF] Francisco | Enamel surface topography analysis for diet discrimination. A methodology to enhance and select discriminative parameters[END_REF].

On a normalized ACF, an horizontal cutting plane at height 0 < s < 1 defines an ellipsis-like line. Defining the minor and major ellipsis radii as Rmax and Sal respectively, the higher Stri=Rmax/Sal the more anisotropic surface. In the present work s = 0.5 and because Stri, Sal and Rmax are linked, Sal is left aside.

Another way of catching anisotropy is to study the ratio r.sl of the maximum of the ACF slope to the minimum: it quantifies the ACF relative decay speed in nearly perpendicular directions.

3 height parameters. Height parameters mainly refer to height statistics.

High order statistical moments are sensitive to outliers so, the height skewness Ssk and kurtosis Sku are chosen; if the healing process introduces abnormal heights or an abnormal height distribution, Ssk and Sku are affected.

An additional height parameter is analyzed, the height mean Sm. Although it is not a relevant parameter for the whole surface, when the surface is sampled, it becomes a wear marker. For instance, when dealing with dental surfaces, it reveals to be a diet discriminative parameter associated to locally stepped surfaces [START_REF] Francisco | Gathering and Analyzing Surface Parameters for Diet Identification Purposes[END_REF].

3.5.3. 2 topological parameters. Topological parameters are defined as parameters determined by morphological operations, such as binarization, dilation, erosion, etc., on surfaces.

Sh is the percentage of quasi-horizontal faces -face normal within a four degrees cone. It is sensitive to the surface long wavelength relief.

The ratio area Sk1 is the fraction of the surface heights above 85% of the height amplitude, it is sensitive to high elevation zones.

3 fractal parameters.

Fractal parameters are defined as parameters that take into account the surface area at different scales.

The relative area Sdar -developed area/projected area -strongly depends on the fractal dimension of the surface: the greater fractal dimension, the higher relative area.

asfc is called "surface complexity", and it is defined as a thousand times minus the steepest slope of the function log(Sdar) = f (log(se)), with se the element size at a given grid scale, see [START_REF] Francisco | Gathering and Analyzing Surface Parameters for Diet Identification Purposes[END_REF] for further details.

indf is a fractal index calculated with the box-counting method. It is the slope of the function log(nb) = f (log(se)), where nb is the number of boxes of size se under the surface at a given scale.

3.5.5.

Surface sampling. As the figures 1,2,3 suggest, the surfaces are not homogeneous. The key idea, as exposed by Francisco et al. [START_REF] Francisco | Enamel surface topography analysis for diet discrimination. A methodology to enhance and select discriminative parameters[END_REF], [START_REF] Francisco | Gathering and Analyzing Surface Parameters for Diet Identification Purposes[END_REF], is to uniformly sample the surface. Different tries have been carried out to determine an optimal number of samples, but experience has proved that it is very parameter-dependent. Indeed, when one tries to separate several diet categories upon a set of parameters, the different numbers of samples make the groups overlap more or less: in the end, it is all about compromise, the samples should not be too small or too large. In the present work, a 16 × 16 image slicing is chosen with 256 × 256 pixels sub-surfaces.

Several statistics are built from the sub-surfaces parameters:

• min.05 the parameter mean of the first 5% values, it is a kind of robust minimum,

• fst.25 the parameter first quartile,

• med the parameter median value,

• lst.25 the parameter third quartile, and

• max.05 the parameter mean of the last 5% values, it is a kind of robust maximum.

Hence, to each parameter calculated on the whole surface -suffixed with glob -5 derived parameters are added.

Analysis procedure

Two different but complementary studies are carried out to answer the questions: what is the effect of the damaged area size on the healing quality? What is the effect of the damaged area location on the healing quality? 3.6.1. Healed area size influence. A 1024 px × 1024 px sub-surface is extracted from each of Anisotropy1-2-3 surfaces. The study could as well be carried out on the 1551 px × 1551 px original surfaces but extracting a sub-surface allows for centering particular features -the deep scratch of Anisotropy1, figure 11, for instance; the impact of healing is thus better tested because the damaged zone is centered. It should be noted that the anistropy ratio is barely altered.

Six hole sizes are tested for healing. The holes are ellipsis-shaped, as on image processing software when defect area are selected. The table 1 details the hole sizes, as well as Da, the ratio of the damaged area on the whole surface area. We investigate here to what extend the defects should be treated with one healing method rather than another. It is expected that whatever the parameter involved, the values on global-, point-or exemplar-healed surfaces, will be closer to the values on the original surface, than for Laplace's infilling. For common sense reasons, the healed surface results should also be better for smaller damaged areas.

a (px) b (px)
Table 2 To ensure that the four healing methods yield very different results, the damaged area should be large enough. However too large an area would be unrealistic. The authors have chosen to randomly scatter four non-overlapping holes for an amount of 2.6 %, figure 10. As it will be seen with the first analysis, Da = 3% is enough for ranking the healing methods. However, a single 3% defect is not accepted in common use, because even if the ratio seems low enough, it is considered too much intrusive. It is the reason why we choose to break the "large" defect into four parts because an experienced user is more likely to accept several smaller defects rather than a large one. In addition, four elements let more location possibilities than a single one.

Then, the healing processes are carried out and the parameters are calculated for comparison with the original surface. The hole locations are not equivalent regarding the healing results, that is the reason why the calculations are performed many times. In order to compare the results using statistical tests, the amount of calculations must be high enough: this way, the tests will be more robust against departure from normality of the result distributions. n = 512 -the number of runs carried out for each case -is found to be a good compromise between statistical significance and computing time.

When performing a statistical test, an hypothesis H 0 is assumed. As an example, to test the normality of a set of data, one can use the well-known Shapiro-Wilk test; H 0 is "the data are normally distributed". For the moment, let us put aside the test assumption requirements. Some results, that will be presented in the result section, illustrate the whole set of n values for a given parameter, on a given surface, for each inpainting method, as well as the parameter value on the original surface. Let us suppose that an inpainting method has a mean value closer to the original one than another inpainting method, but that it has a greater variance in the mean time. Which one of the two is the best choice to fill in the holes? This is where the statistical tests can help. In particular, the one-sample Student t-test explores the possibility for a set of data to have a given value µ 0 as mean value. In the present study, the set of data is the n runs of four holes randomly scattered on a surface, and µ 0 is the parameter value calculated on the original surface.

One of the outputs of the statistical test is the well-known p-value, which is the probability to find sets of data more extreme that the one studied, assuming that H 0 is true, i.e. that the set mean value is µ 0 . In other words, the p-value quantifies how well the data fit H 0 . The p-value is a controversial subject because it is sometimes interpreted as how much the theory H 0 is truthfully; however, it is used here to compare the data distributions to a given mean value µ 0 , assuming H 0 .

As stated before, a statistical test requires more or less assumptions. The Student t-test assumes that the data are independent, normally distributed with no outlier. In a more practical way, as the t-test is fairly robust against departures from normality, especially in large samples [START_REF] Dalgaard | Introductory Statistics with R. Statistics and Computing[END_REF], we seek for near-normality and the 2% extreme values are excluded from the sample. It is to be noted that arbitrarily removing data from the data set is not recommended by some authors, as Quinn [START_REF] Quinn | Experimental Design and Data Analysis for Biologists[END_REF] p.68: [...] dropping observations just because they are messy or reduce the chance of getting a significant result is unethical, to say the least. The authors of the present study agree this statement; however here, potential outlier data are not chosen but rather systematically dropped.

For maximal chances to get near-normality, the data can be first monotonously transformed.

The Box-Cox transformation is a mathematical function that is automatically tuned to make the data closer to normality as possible. Even if it is of low interest because the data set is large -almost 500 values -the trick is anyway applied. In some cases the result distributions are too much skewed to just ignore it.

To end with the analysis protocol, as some results exhibit multi-modality, an additional statistical test is used. The one-sample Wilcoxon rank test is a non-parametric test, sometimes described as the alternative to the Student t-test when the data cannot be assumed to be normally distributed. This time the data median value is compared to the original value. The data being ranked, there is no need for normality but just symmetry. Once again with so large a data set, reasonable departures from perfect symmetry are allowed.

It is expected here that one of global-, point-or exemplar-healing turns to be a statistically much better choice than Laplace's because the results are visually far more realistic. Laplace-based method seems visually yield very poor results, while the other methods create realistic features. The results are quite impressive but to what extend the reconstructed area can be trusted? To put it simple, can we accurately replace any real area by a numeric patch? The answer is intuitively no, and it can be numerically brought to light with surface parameter comparisons. 2. Hereafter -figure 14 -are some representative cases of the parameters calculated on the whole surface Anisotropy1 : each plot synthesizes the results of a single parameter calculated on a whole particular surface, for the different hole sizes and the different methods. We choose to illustrate the results on the whole surface rather than a particular statistics (min.05, fst.25, med, lst.25, max.05 ) to help the understanding. In ESM(ESM1) almost all remaining parameters are presented on the same surface, Anisotropy1.

Results and discussion

As it clearly appears on the figure 14, the largest damaged areas -Da 12%yield the poorest results. It was nonetheless foreseeable that an artificial process could not reproduce such an important patch, despite visually convincing results.

The anisotropy factor Stri is not more than 3.5% different from the theoretical value if Da 6%. Laplace's operator shows rather good results, which could appear surprising at first sight, but as it can be seen on figure 11, it introduces no local parasitic anisotropy. Therefore, the overall anisotropy remains well caught by the autocorrelation function.

The mean height Sm seems to suffer more from surface healing, but a look at the scale reveals very small values, the variation remaining within a few nanometers. Once discarded the Da 12% restoration cases, the healing processes do not alter much Sm, Laplace's being the best choice. It is simply explained by the fact that the theoretical patch has a mean value very close to the whole surface; thus, Laplace's smooth patch barely modifies it.

The flatness ration Sh is well preserved with point-and exemplar-based methods, the error being around 0.8% when Da 6%. This is due to the method main characteristic: the reproduction of the neighboring relief. As concerns Laplace's method efficiency here, the reason is different. The patch being already almost 74% horizontal, Laplace's operator doesn't alter much Sh. In ESM(ESM2), the reader can find:

• the results of the figure 14 cleared from the two cases Da 12%,

• the results for the two other surfaces Anisotropy2-3 showing lower global flatness, and poorer results for Laplace's method.

It is out of the question to choose in practice the healing method for a given analysis because:

• there are often different kinds of parameters associated to a surface study,

• the methods are always improvable, except for Laplace's which is a particular equation implementation,

• the quality of the results depends obviously on the surface features, on their homogeneity, etc.

Here however, we can rank the methods upon their ability to recreate parts of a surface by averaging the results:

• on the three surfaces Anisotropy1-2-3,

• with all Da values,

• of all sub-parameters.

The reader is referred to ESM(ESM3) for a short explanation on the ranking process.

Doing so, the ranks shown on figure 15 highlight the advantages in using a pointbased technique to heal damaged surfaces. As warned before, restoring images with Da 12% is not realistic but, if the two related cases -12% and 24% -are excluded from the results, the graphs remain approximately the same. Additional results, that detail the synthetic results, are also provided in ESM(ESM3).

Defect location influence

When surfaces are not perfectly homogeneous, the damaged area locations are not equivalent from a healing point of view. In addition, as exposed in the method section, several moderate size damaged areas are more likely to be accepted than a single important one. The reasons for multiple defects to occur are

• the real surface is not clean enough: there are remaining debris on the surface • the microscope light is locally badly reflected: when the surface exhibits very steeped faces, the reflected rays are not caught

• the surface mold has trapped air bubbles, etc.

The figure 16 synthesizes the results of the 512 runs carried out on the surface Anisotropy1 for the same parameters as in the previous section, and a damaged area Da = 2.6%. See ESM(ESM4) for Anisotropy2-3. The first observation to be made is the pretty good symmetry of the data which safely allows one-sample Wilcoxon rank tests. Even if the data do not follow a normal law, thanks to the sample number and the low skew coefficient it is also expected Shapiro-Wilk tests to be applicable. Depending on the location of the damaged areas, Laplace's operator can significantly alter the anisotropy factor Stri : the range of values is at least twice than for the other methods. As a consequence, there are many points with a 10% error, even though the mean value remains close to the theoretical value. This statistical study gives clearly advantage to the exemplar-based method which contradicts the synthetic results of the figure 15. But, having a closer look on the more detailed results, figure 14 with Da = 3%, the exemplar-based method can provide the theoretical value. Actually, there is no real contradiction here: the mean rank results take into account the two other anisotropy parameters as well as their sub-parameters, i.e. fst.25, min.05, etc. Moreover, in the previous section, the damaged area was centered so the conclusions must remain cautious. The global method, which aimed to be particularly suited for anisotropic surfaces, yields at the end a bit disappointing results.

Concerning Sm results -the height mean -Laplace's patches approximate well the real areas. The same explanations as in the previous section hold. The two other parameters, the topological parameter Sh and the fractal parameter asfc, have less close values to the right values with Laplace's method than with other methods.

Figure 17 proves the meaningfulness of the statistical tests: despite some minor differences, the parametric test (Shapiro-Wilk) and non-parametric test (Wilcoxon) give the same results. The exemplar-based method -except for the height parameters -and the point-based method are the best way to heal the damaged surfaces. Once again, this statement does not contradict with the previous section observations because the situation is different. First, in the present situation Da 3% and it is summed over four areas, second, figure 17 is based on statistics; therefore many different cases are averaged, unlike the previous single cases. The reader may also refer to ESM(ESM5) for less synthetic results. When averaging the results on the surfaces Anisotropy1-2-3 and the two tests, the inpainting methods can be ranked, as proposed on figure 18. The conclusion is straightforward because the graphs show the best performance for point-and exemplarbased methods. If we take into account the previous section conclusions, it can be safely stated that the visual feelings about the different process quality are correlated with the quantitative analyses. Let us however remain careful about the superiority of the methods. The rank scales are not the same for the different parameter categories. The ranks are somewhat close to each other in the height graph, which means that Laplace yields not so much different results than the other methods.

Conclusion

The concluding remarks are the following:

• From a visual point of view, Laplace's results are the worst

• Despite very convincing visual aspects, except for Laplace, the healing processes significantly alter the results when the damaged area exceeds 3%

• Developing a method to specifically handle anisotropic surfaces turns not to be useful, compared to common methods like point-and exemplar-based methods

• Two categories of parameters are particularly sensitive to Laplace's inpainting: topology and fractal

• Statistics show that point-and exemplar-based methods produce better results than Laplace's method.

Point-and exemplar-based methods have been developed twenty years ago and are still improvable. They are fast, but for digital surface healing, speed is not a key point, accuracy is. We recognize that for scattered unread points -e.g. because of high local slopes -Laplace remains well suited. But, we think that it is worth testing more advanced methods like wavelet decomposition methods, or more recent ones like AI-based methods, when large damaged area occur. Otherwise, for moderate damaged areas, point-and exemplar-based methods remain efficient tools.

Figure 1 :

 1 Figure 1: Bos taurus, anisotropy: high

  Figure 2 exhibits longer wavelengths than Anisotropy1, but the anisotropy factor is smaller.

Figure 2 :

 2 Figure 2: Cephalophus silvicultor, anisotropy: moderate

Figure 3 :

 3 Figure 3: Cervus elaphus, anisotropy: low

Figure 4 :

 4 Figure 4: The points A, B, C and D are not coplanar, (AB) and (CD) do not intersect.There is no reason to linearly interpolate [AB] rather than [CD]: if a central point is to be drawn, it is at a fraction of the distance between (AB) and (CD)

Figure 8 :

 8 Figure 8: If each new ∂Ω is filled in one single step, the iterative borders evolve towards the hole center. When the center is reached, a gap is likely to appear between the bottom line and the top line.

Figure 10 :

 10 Figure 10: 1551 × 1551 surface with four elliptic holes covering 2.6 % of the surface area. The same surface ratio is presented in a single elliptic shape in the bottom left. A 256 × 256 sample is also represented.

4. 1 .

 1 Healed area size influence 4.1.1. Visual results. The following figures -from 11 to 13 -show the infilling results on the three surfaces Anisotropy1-2-3, with a central hole; the hole relative area is Da = 12 %.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: The original surface (O) -Anisotropy1 -is shown in central position. An elliptic area is erased on the surface center, and four infilling methods are tested: Laplace (L) top left, for which the infilled area clearly appears, Global-based (G) top right, Point-based (P) bottom left and Exemplar-based (E) bottom right

Figure 14 :

 14 Figure 14: Four different category parameters calculated on the surface Anisotropy1 are presented. Top left, the anisotropy factor Stri ; top right, the mean height Sm; bottom left, the horizontality ratio Sh; bottom right the area scale factor asfc. The healing methods are: (E) the exemplar-based method; (G) the global method; (L) Laplace's method; (P) the point-based method. The different hole relative sizes are expressed as percentages.

Figure 15 :

 15 Figure 15: The four ranked categories are presented. Top left, the anisotropy category; top right, the height category; bottom left, the topology category; bottom right the fractal category. The healing methods are: (E) the exemplar-based method; (G) the global method; (L) Laplace's method; (P) the point-based method. The lower rank, the better performance.

Figure 16 :

 16 Figure 16: Four different category parameters calculated on the surface Anisotropy1 are presented (Da = 2.6%). Each circle is a random location of the holes. Top left, the anisotropy factor Stri ; top right, the mean height Sm; bottom left, the horizontality ratio Sh; bottom right the area scale factor asfc. The healing methods are: (E) the exemplar-based method; (G) the global method; (L) Laplace's method; (P) the point-based method.

Figure 17 :

 17 Figure 17: The four ranked categories are presented. Top left, the anisotropy category; top right, the height category; bottom left, the topology category; bottom right the fractal category. The healing methods are: (E) the exemplar-based method; (G) the global method; (L) Laplace's method; (P) the point-based method. Two statistical tests are used: parametric (SH) and non-parametric (WI)

Figure 18 :

 18 Figure 18: The four ranked categories are presented. Top left, the anisotropy category; top right, the height category; bottom left, the topology category; bottom right the fractal category. The healing methods are: (E) the exemplar-based method; (G) the global method; (L) Laplace's method; (P) the point-based method.

Table 1 :

 1

			relative area Da (%)
	71	35	0.75
	100	50	1.50
	141	71	3.00
	200	100	6.00
	283	142	12.00
	400	200	24.00

Six hole sizes are tested to analyze the defect size influence on the healing results. a and b are the semi-major and semi-minor axes respectively.

Table 2 :

 2 recaps the whole parameter set; it represents more than 4750 cases. Recap of the different case combinations that are computed 3.6.2. Statistical analysis.

	Surface	Da Healing method Parameters	Sampling statistics
	Anisotropy1 0.75	(E)xemplar	Stri		glob
	Anisotropy2 1.50	(G)lobal	Rmax	Spatial	min.05
	Anisotropy3 3.00	(L)aplacian	r.sl		fst.25
		6.00	(P)oint-based Ssk		med
		12.00		Sku	Height	lst.25
		24.00		Sm		max.05
				Sh Sk1	Topology
				Sdar
				asfc	Fractal
				indf


Adel Mesbah

Jai Prakash

Sébiastien Lebègue

Jessica C Beard

Christos D Malliakas

James A Ibers

Ag5U(PS4)3: A Transition-Metal Actinide Phosphochalcogenide

The structure of Ag5U(PS4)3 is unique as in the literature there are no other structures of the type MAnPQ (M = transition metal, An = actinide, Q = S, Se, or Te). The compound has been synthesized at 1123 K by standard solid-state methods and its single-crystal X-ray structure has been determined at 100(2) K. Ag5U(PS4)3 crystallizes in a remarkable new structure type in space group P3221 of the trigonal system with three formula units in a hexagonal cell of dimensions a = b = 9.6635(2) Å, c = 17.1834(4) Å, and γ = 120. In the structure each U atom is coordinated to eight S atoms in a bicapped trigonal prismatic manner. Each P atom is tetrahedrally coordinated to four S atoms. Two of the three unique Ag atoms are connected to four S atoms in a distorted tetrahedral manner whereas the third unique Ag atom forms a Ag2S6 species. The overall structure consists of U polyhedra connected to each other via PS4 tetrahedra through edge-sharing in a zig-zag fashion along the c axis to form infinite layers.

PS4 groups and the Ag atoms pack these layers.

From DFT calculations Ag5U(PS4)3 is found to have finite spin polarization in the crystal cell: the magnetic moments of two of the U atoms are parallel, whereas the magnetic moment of the third U atom is antiparallel. l

INTRODUCTION

The search for new actinide chalcogenides (An = Th, U, Np; Q = S, Se, Te) has been widely conducted through the use of solid-state methods involving a flux method or direct combination of elements. [START_REF] Mesbah | Overview of the crystal chemistry of the actinide chalcogenides: incorporation of the alkaline-earth elements[END_REF][START_REF] Manos | The Chemistry of the Actinide and Transactinide Elements[END_REF][START_REF] Bugaris | Syntheses and characterization of some solid-state actinide (Th,U,Np) compounds[END_REF][START_REF] Narducci | Ternary and Quaternary Uranium and Thorium Chalcogenides[END_REF][START_REF] Koscielski | The Structural Chemistry of Quaternary Chalcogenides of the Type AMM'Q3[END_REF] The phosphochalcogens offer structural flexibility that often leads to a variety of exciting new structures. The greater diversity of stoichiometries among the phosphochalcogens arises from their ability to form Q-Q bonds that leads to polymeric substructures such as PQ3, P2Q6, P2Q10, P4Q13, P6Q12, and PQ6. [START_REF] Chondroudis | Rb4Sn2Ag4(P2Se6)3: First Example of a Quinary Selenophosphate and an Unusual Sn-Ag s 2 -d 10 Interaction[END_REF][START_REF] Gieck | Supramolecular" Solid-State Chemistry: Interpenetrating Diamond-Type Frameworks of U 4+ Ions Linked by S,S'-Bidentate P2S6 2-Molecular Rods in UP4S12[END_REF][START_REF] Hess | Synthesis and Characterization of Six New Quaternary Actinide Thiophosphate Compounds: Cs8U5(P3S10)2(PS4)6, K10Th3(P2S7)4(PS4)2, and A5An(PS4)3[END_REF][START_REF] Wu | Rb3Ti3(P3S13)(PS4)3 and Cs2Ti2(P2S8)(PS4 ) : Two Polar Titanium Thiophosphates with Complex One-Dimensional Tunnels[END_REF][START_REF] Cody | Ionothermal Synthesis of Four New Nickel Thiophosphate Anions: [Ni(P2S8)2] 2[END_REF] Known phosphosulfides include ternaries such as UP1-xSx, [START_REF] Baskin | Melting Point Studies of the System UP-US[END_REF] UPS, [START_REF] Kaczorowski | Crystal Structure, Magnetic and Electrical Transport Properties of UPS Single Crystals[END_REF] ThP2S6, [START_REF] Simon | Darstellung und Kristallstruktur von ZrP2S6 und ThP2S6[END_REF] UP2S6, 14 U(P2S6)2, 14 UP2S7, [START_REF] Neuhausen | Comprehensive Uranium Thiophosphate Chemnistry: Framework Compounds Based on Pseudotetrahedrally Coordinated Central Metal Atoms[END_REF][START_REF] Babo | Structure, and Spectroscopy of Two Ternary Uranium(IV) Thiphosphates: UP2S9 and UP2S7 Containing P2S9 2and P2S7 2-Ligands[END_REF] UP2S9, [START_REF] Babo | Structure, and Spectroscopy of Two Ternary Uranium(IV) Thiphosphates: UP2S9 and UP2S7 Containing P2S9 2and P2S7 2-Ligands[END_REF] and U3(PS4)4; [START_REF] Neuhausen | Comprehensive Uranium Thiophosphate Chemnistry: Framework Compounds Based on Pseudotetrahedrally Coordinated Central Metal Atoms[END_REF] quaternaries such as A11U7(PS4)13 (A = K, Rb), [START_REF] Gieck | Interlocking Inorganic Screw Helices: Synthesis, Structure, and Magnetism of the Novel Framework Uranium Orthothiophosphates A11U7(PS4)13 (A = K,Rb)[END_REF] CsLiU(PS4)2, 17 Cs8U5(P3S10)2(PS4)6, 8 A5An(PS4)3 (A = K, Rb, Cs; An = U, Th), [START_REF] Hess | Synthesis and Characterization of Six New Quaternary Actinide Thiophosphate Compounds: Cs8U5(P3S10)2(PS4)6, K10Th3(P2S7)4(PS4)2, and A5An(PS4)3[END_REF] AkAn(PS4)2 (Ak = Sr, Ba; An = Th, U), [START_REF] Mesbah | Crystal Structures, Optical and Theoretical Studies of the Actinide Thiophosphates SrU(PS4)2, BaU(PS4)2, and SrTh(PS4)2[END_REF] and Cs4Th2P6S18; [START_REF] Chan | Synthesis and Characterization of Two Quaternary Thorium Chalcophosphates: Cs4Th2P6S18 and Rb7Th2P6Se21[END_REF] and the quintaries or greater such as A6U3Sb2P8S32, [START_REF] Babo | A6U3Sb2P8S32 (A = Rb, Cs): Quinary Uranium (IV)Thiophosphates Containing the [Sb(PS4)3] 6-Anion[END_REF] Cs5Na6[U(PS4)4](PS4), [START_REF] Klepov | Complex Topologies from Simple Building Blocks: Uranium(IV) Thiophosphates[END_REF] Rb5Na3[U(PS4)4], [START_REF] Klepov | Complex Topologies from Simple Building Blocks: Uranium(IV) Thiophosphates[END_REF] CsNa[U(PS4)2], [START_REF] Klepov | Complex Topologies from Simple Building Blocks: Uranium(IV) Thiophosphates[END_REF] Cs1.67Na0.52I0. 19 [START_REF] Klepov | Complex Topologies from Simple Building Blocks: Uranium(IV) Thiophosphates[END_REF] But there are also examples of phosphoselenides and less so of phosphotellurides. These include AnPSe (An = U, Th), [START_REF] Oh | The beta-polymorph of uranium phosphide selenide[END_REF][START_REF] Hulliger | New Ternary Thorium and Uranium Compounds MYX[END_REF] K2U(P3Se9), 24

Rb4U2(P2Se6)2.5(Se2), [START_REF] Bellott | Synthesis and Structure of the Rubidium Uranium Selenophosphate Rb4U2P5Se17[END_REF] Cs2An2(P2Se9)(Se2)2 (An = U, Th), [START_REF] Briggs Piccoli | Structure, and Properties of Cs4Th4P4Se26: A Quaternary Thorium Selenophosphate Containing the (P2Se9) 6-Anion[END_REF][START_REF] Mesbah | Synthesis and Crystal Structure of Cs2U2(P2Se9)(Se2)2[END_REF] Rb7Th2(PSe4)3(P2Se6)1.5, 19

A2Th(P3Se9) (A = K, Rb), 28 Cs4Th2P5Se17 28 for the selenides and only AnPTe (An = U, Th) for tellurides. [START_REF] Stolze | ThPTe and U2PTe2O: actinide pnictide chalcogenides with diphosphide anions[END_REF] In the present paper, we report the synthesis, crystal structure, and electronic structure of the compound Ag5U(PS4)3. Its structure is the first example of an actinide transition-metal phosphochalcogenide.

EXPERIMENTAL METHODS

Synthesis and Analyses. Caution! Depleted U is an α-emitting radioisotope and as

such is considered a health risk. Experiments using radioisotopes require appropriate infrastructure and personnel trained in the handling of radioactive materials.

The following reactants were used as supplied: Ba (Johnson Matthey, 99.5%), P2S5

(Aldrich, 99%), S (Mallinckrodt 99.6%), and Ag (Aldrich, 99.99%). Depleted U powder was obtained after hydridization of U chips (IBI Laboratories) in a modification [START_REF] Bugaris | structures, and magnetic and optical properties of the compounds [Hg3Te2[END_REF] Structure Determination. Single-crystal X-Ray diffraction data for Ag5U(PS4)3 were collected at 100(2) K using a Bruker APEX2 Kappa diffractometer equipped with the CuKα (λ = 1.54178 Å) radiation. By the use of the algorithm COSMO in APEX2 [START_REF]5-1 Data Collection and Processing Software[END_REF] the data collection comprised a combination of ω and φ scans with steps of 0.3º and counting time of 10 s/frame.

Recorded data were indexed, refined, and integrated by SAINT in the APEX2 package. [START_REF]5-1 Data Collection and Processing Software[END_REF] The crystal-to-detector-distance was 40 mm. Face-indexed absorption, incident beam, and decay corrections were performed with the use of the program SADABS. [START_REF] Sheldrick | [END_REF] Precession images constructed from the data set showed no indication of super cells or modulation. The structure was difficult to solve. Ultimately, the program XT [START_REF] Sheldrick | SHELXT-Integrated space-group and crystal-structure determination[END_REF] provided a solution that could be refined with the use of programs in the SHELXTL 2014 package. [START_REF] Sheldrick | [END_REF][START_REF] Sheldrick | A short history of SHELX[END_REF] The asymmetric unit contains one U, three Ag, two P, and 6 S atoms, with atoms U1 and P2 having .2. symmetry. The anisotropic refinement of the atoms suggests some disorder in the Ag2 atoms. The occupancy of the Ag1 site was not constrained by symmetry and ultimately refined to 0.503 (3). Consequently, this occupancy was set to 0.5 to give the formula Ag5UP3S12, which is consistent with the EDX results. Further details are provided in Table 1and in Supporting Information.

l Theoretical Calculations. These have been conducted with density functional theory [START_REF] Kohn | Self-Consistent Equations Including Exchange and Correlation Effects[END_REF][START_REF] Hohenberg | Inhomogeneous Electron Gas[END_REF] using the projector augmented wave method [START_REF] Blöchl | Projector augmented-wave method[END_REF] as implemented in the VASP (Vienna ab Initio Simulation Package) [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] code. The generalized gradient approximation of Perdew, Burke and Ernzerhof [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] in its spin polarized version was used as the exchange and correlation potential. The parameters of the cell and the positions of the atoms were taken from the experiment, but because the experimental structure contains one Ag site that is half occupied we have created several models consistent with the chemical formula Ag15U3P9S36 in the following way: all the equivalent positions of this site were generated and then models were produced after populating half the Ag1 atoms while keeping the remaining unoccupied. This can be done in different ways, and we have retained the one with the lowest total energy for our analysis. The criterion for the convergence of the total energy was set to 10 -6 eV To obtain a sufficient convergence, a 2 × 2 × 1 mesh was used for the Brillouin zone. The default value was used for the cut-off of the plane wave part of the wave function.

RESULTS

Synthesis. The reaction of Ba, Ag, U, P2S5, and S at 1123 K resulted in the formation of black blocks of Ag5U(PS4)3 in a yield of about 50 wt% yield based on the U content. Side products included UP2S6 and BaS. The isolated crystals of Ag5U(PS4)3 were stable in oil for a few weeks. Unfortunately, by the time we knew the composition and the nature of the structure, the isolated crystals had decomposed.

Structure.

A search of the ICSD data base failed to reveal any structures of the type MAnPQ (M = transition metal, An = actinide, Q = S, Se, or Te); however, it did display the isostructural compounds AgLn(P2Se6) (Ln = rare earth = Sc, Er, Tm). [START_REF] Pfitzner | Synthesis and Structure Determination of AgScP2Se6, AgErP2Se6 and AgTmP2Se6[END_REF] The structure of these isostructural compounds is layered and very different from that of Ag5U(PS4)3. Ag5U(PS4)3 crystallizes in a remarkable new structure type in space group P3221 of the trigonal system with three formula units in a hexagonal cell of dimensions a = b= 9.6635(2) Å, c = 17.1834(4) Å, and γ = 120. The asymmetric unit contains one U, three Ag, two P, and six S atoms. All are in general positions except atoms U1 and P2, which have .2. site symmetry.

A projection down the [110] direction of the Ag5U(PS4)3 crystal structure is shown in Figure 1. Each U1 atom is coordinated to eight S atoms in a bicapped trigonal prismatic manner.

Each P atom is tetrahedrally coordinated to four S atoms. Each of the Ag1 and Ag3 atoms is connected to four S atoms in a distorted tetrahedral manner. In contrast, the Ag2 atoms exhibit a Ag2-Ag2 short distance of 2.774(5) Å and form a Ag2S6 species. The Ag-S network is presented in Figure 2. The overall structure consists of U polyhedra connected to each other via P2S4 tetrahedra through edge-sharing in a zig-zag fashion along the c axis to form infinite layers. The P1S4 and the Ag atoms pack these layers.

The U1-S distances in the structure of Ag5U(PS4)3 range from 2.784(5) to 2.851(5) Å, these values are typical for U 4+ coordinated to eight S atoms, as in the structures of K5U(PS4)3 [START_REF] Hess | Synthesis and Characterization of Six New Quaternary Actinide Thiophosphate Compounds: Cs8U5(P3S10)2(PS4)6, K10Th3(P2S7)4(PS4)2, and A5An(PS4)3[END_REF] (2.802(2) to 3.021(2) Å); CsLiU(PS4)2 [START_REF] Neuhausen | Modular Metal Chalcogenide Chemistry: Secondary Building Blocks as a Basis of the Silicate-Type Framework Structure of CsLiU(PS4)2[END_REF] (2.768(1) to 2.884(1) Å); K11U7(PS4)13 [START_REF] Gieck | Interlocking Inorganic Screw Helices: Synthesis, Structure, and Magnetism of the Novel Framework Uranium Orthothiophosphates A11U7(PS4)13 (A = K,Rb)[END_REF] [START_REF] Hu | I) Multiple Salts Assembled with Phosphonate, Perfluorocarboxylate, and the Multinuclear Silver-Ethynediide Synthon[END_REF] .

The crystal structure of the Ag5U(PS4)3 has no S-S short distances and thus the assignment of formal oxidation states is straightforward: 5 × Ag + , 1 × U 4+ , 3 × P 5+ , 12 × S 2-.

DFT Calculations. Ag5U(PS4)3 is found to have finite spin polarization in the crystal cell: the magnetic moments of two of the U atoms are parallel, whereas the magnetic moment of the third U atom is antiparallel. This is seen in Figure 3 where the total density of states (top plot) is asymmetric with respect to spin. The partial density of states for each inequivalent atom in the cell is also presented in Figure 3. Most of the magnetic moment is carried by the U atoms l which induces a small magnetization on the neighboring atoms. Also, the states at the Fermi level mostly arise from the U. From DFT calculations Ag5U(PS4)3 is found to have finite spin polarization in the crystal cell: the magnetic moments of two of the U atoms are parallel, whereas the magnetic moment of the third U atom is antiparallel.

CONCLUSIONS

ASSOCIATED CONTENT

Supporting Information

The crystallographic CIF file for Ag5U(PS4)3 has been deposited with CCDC (https://www.ccdc.cam.ac.uk) as entry CCDC1866631. c Rw(Fo 2 ) = {Σ [w(Fo 2 -Fc 2 ) 2 ] / ΣwFo 4 } 1/2 . For Fo 2 < 0, w -1 = σ 2 (Fo 2 ); for Fo 2 0, w -1 = σ 2 (Fo 2 ) + (qFo 2 ) 2 where q = 0.0321.

Table 2. Interatomic Lengths (Å) in Ag5U(PS4)3.

Atom-Atom Distance (Å)

Atom-Atom Distance (Å)

Ag5U(PS4) 3

 3 has a unique structure as there appear to be no other known structures of the type MAnPQ (M = transition metal, An = actinide, Q = S, Se, or Te). Ag5U(PS4)3 crystallizes in a remarkable new structure type in space group P3221 of the trigonal system with three formula units in a hexagonal cell of dimensions a = b= 9.6635(2) Å, c = 17.1834(4) Å, and γ = 120. The asymmetric unit contains one U, three Ag, two P, and six S atoms. All are in general positions except atoms U1 and P2, which have .2. site symmetry. In the structure each U atom is coordinated to eight S atoms in a bicapped trigonal prismatic manner. Each P atom is tetrahedrally coordinated to four S atoms. Two of the three unique Ag atoms are connected to four S atoms in a distorted tetrahedral manner whereas the third Ag atom forms a Ag2S6 species.The overall structure consists of U polyhedra connected to each other via PS4 tetrahedra through edge-sharing in a zig-zag fashion along the c axis to form infinite layers. PS4 groups and the Ag atoms pack these layers. The crystal structure of the Ag5U(PS4)3 has no S-S short distances and thus the assignment of formal oxidation states is straightforward: 5 × Ag + , 1 × U 4+ , 3 × P 5+ , 12 × S 2-.

Figure 1 .

 1 Figure 1. A projection down the [110] direction of the Ag5U(PS4)3 crystal structure.

Figure 2 .

 2 Figure 2. The Ag-S network of the Ag5U(PS4)3 crystal structure.

Figure 3 .

 3 Figure 3. Total (upper plot) and partial density of states (lower plots) of Ag5U(PS4)3.

Table 1 . Crystallographic Data and Structure Refinement Details for Ag5U(PS4)3.

 1

	Compound	Ag5U(PS4)3
	Space group	P3221
	fw (g.mol -1 )	1255.01
	a (Å)	9.6635(2)
	c (Å)	17.1834(4)
	V (Å 3 )	1389.66(7)
	Z	3
	T (K)	100(2)
	λ (Å)	1.54178
	ρ (g cm -3 )	4.499
	μ (mm -1 )	80.8
	R(F) b	0.04
	Rw(Fo 2 ) c	0.12

b R(F) = ΣFoFc/ ΣFo for Fo 2 > 2σ(Fo 2 ).

ACKNOWLEDGMENTS

Use was made of the IMSERC X-ray Facility at Northwestern University, which has received support from the Soft and Hybrid Nanotechnology Experimental (SHyNE) Resource (NSF ECCS-1542205); the State of Illinois, and International Institute for Nanotechnology (IIN). We have benefited from discussions with Prof. George Sheldrick concerning the crystal structure.

l

AUTHOR INFORMATION

*E-mail: ibers@chem.northwestern.edu

Notes

The authors declare no competing financial interest.


Guillaume Picard

Roland Lenain

Jean Laneurit

Benoit Thuilot

Christophe Cariou

A predictive control framework for edge following: Application to two types of mobile robots

In this paper, the accurate control of mobile robots is investigated in the framework of generic edge following. It proposes a new predictive control approach, based on the minimization of the lateral error along a distance of prediction. This permits to consider a distance of convergence independently from the yaw dynamics and the longitudinal velocity, allowing to achieve harsh maneuvers, with respect to potentially kinetically unachievable paths. The proposed algorithm is generic and permits to address the control of different kinds of robots (skid-steering, car-like, etc.) in a common framework and to consider independently the speed regulation. The control proposed here allows an accurate and reactive path tracking, even if the environment is complex and narrow. The efficiency of the approach is investigated through full-scale experiments in various conditions.

I. INTRODUCTION

The rising popularity of autonomous driving [START_REF] Marmoiton | Toward Smart Autonomous Cars[END_REF] shows the huge progresses achieved in the field of mobile robots. Beyond the transportation area, many applications may benefit from the advances in such a topic, such as workshop, defense or agriculture [START_REF] Bergerman | Robotics in Agriculture and Forestry[END_REF]. In this latter domain, the use of robots rises as a preponderant issue, since the level of production must increase while the environmental impacts have to decrease significantly, especially with respect to the use of chemical products [START_REF] Hunter | Agriculture in 2050: Recalibrating targets for sustainable intensification[END_REF]. The use of robots able to act autonomously to pull an implement or to carry out heavy loads may permit, among others, mechanical treatments while reducing the hardness of work [START_REF] Blackmore | Towards robotic agriculture[END_REF]. Nevertheless, considering autonomous motion for such an application requires a high level of accuracy, despite various conditions of evolution.

Many developments in this area are based on trajectory tracking with respect to an absolute reference [START_REF] Samson | Modeling and control of wheeled mobile robots[END_REF]. In open areas, this can be easily achieved using a GPS sensor, since it offers permanently an absolute localization for the robot navigation [START_REF] Cariou | Autonomous maneuver of a farm vehicle with a trailed implement: motion planner and lateral-longitudinal controllers[END_REF]. Nevertheless, the variety of encountered situations does not permit to use such a sensor in all cases. As a consequence, alternatives must be developed to address such a variability. Among others, edge following based on the detection of a structure can be considered. This is particularly suitable for agricultural applications, since the robot may have to follow rows of vegetation (as it is the case for precision spraying [START_REF] Gonzalez-De Soto | Autonomous systems for precise sprayingevaluation of a robotised patch sprayer[END_REF]) or when considering autonomous motion inside the farm, where the robot may have to 1 Université Clermont Auvergne, INRAe, UR TSCF, 63178 Aubière, France FirstName.Name@inrae.fr 2 Université Clermont Auvergne, CNRS, SIGMA Clermont, Institut Pascal, F-63000 Clermont-Ferrand, France benoit.thuilot@uca.fr follow walls or furniture (such as in cleaning or feeding applications [START_REF] Bergerman | Robotics in Agriculture and Forestry[END_REF]). For that purpose, numerous developments are based on artificial vision such as considered in [START_REF] Cherubini | Visual navigation of a mobile robot with laser-based collision avoidance[END_REF] or [START_REF] Astrand | A vision based row-following system for agricultural field machinery[END_REF] for agricultural applications or based on a Lidar such as achieved in [START_REF] Bayar | Localization and control of an autonomous orchard vehicle[END_REF] or [START_REF] Tourrette | Tracking footprints for agricultural applications: a low cost lidar approach[END_REF]. Based on this local reference, different control laws can then be applied (such as derived in [START_REF] D'andréa-Novel | Control of nonholonomic wheeled mobile robots by state feedback linearization[END_REF] or [START_REF] Lenain | Robust sideslip angles observer for accurate off-road path tracking control[END_REF]). Nevertheless, in such approaches, angular and lateral dynamics are not explicitly controlled, so that they are not suitable when moving autonomously in narrow spaces. In [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF], a backstepping approach has been proposed to achieve edge following with an explicit computation of the robot position, heading and control variables, allowing to consider independently longitudinal and lateral controls and separate the position and the orientation servoing of the robot. It obtains satisfactory results as long as the variations in the structure to be followed are smooth and the robot does not meet limitations (such as saturation or high dynamics). When the structure to be followed includes harsh variations, some overshoots may occur, depreciating the accuracy of the motion.

In this paper, a model predictive control [START_REF] Richalet | Elementary predictive functional control: A tutorial[END_REF] is derived, based on the backstepping approach previously developed. An optimal desired heading is computed to minimize the future lateral error along a horizon of prediction. Classical predictive control laws developed in the framework of mobile robotics (such as in [START_REF] Prez-Morales | Multi-Sensor-Based Predictive Control for Autonomous Backward Perpendicular and Diagonal Parking[END_REF] for autonomous parking or in [START_REF] Krid | A new explicit dynamic path tracking controller using generalized predictive control[END_REF] to compensate for inertia and bad grip conditions) minimize with respect to the duration. In the control approach proposed in this paper, the minimization is carried out with respect to a distance along the detected edge. This permits to consider independently the lateral and the longitudinal dynamics of the robot, allowing to achieve harsh maneuvers without a loss of accuracy. This paper is decomposed as follows. First, the model of a mobile robot with respect to an edge is presented. Based on this generic modeling, the backstepping approach previously proposed is briefly recalled. The proposed predictive approach is then detailed. Finally, the efficiency of the proposed control architecture is investigated through fullscale experiments in various conditions, using several kinds of robots.

II. MOTION MODELING

A. Assumptions and notations

This paper aims at proposing a predictive control architecture for mobile robots, dedicated to edge following. As a result, the motion is described with respect to a detected edge Γ, such as depicted in Figure 1. As classically achieved in mobile robotics, it is assumed that the robot is symmetrical, moves on a flat ground, and its wheels roll without sliding on the ground. Without loss of genericity, the point O to be controlled is chosen as the center of the driving wheels for skid-steered robots or the center of the rear axle for car-like mobile robots, see Figure 1. The distance between the point O and the closest point to O belonging to Γ (point A in Figure 1) is the lateral deviation denoted by y. The difference between the robot heading and the orientation of the tangent to Γ at point A constitutes the angular deviation and is denoted by θ . The longitudinal position along the trajectory is described by the curvilinear abscissa s. Finally, the curvature of Γ at the curvilinear abscissa s is denoted by c(s). The robot's motion is achieved thanks to two control variables:

• the robot longitudinal velocity v,

• the robot angular velocity u θ , which is directly controlled from the driving wheels velocities in the case of skid-steered mobile robots or linked to the front steering angle δ F in the case of car-like mobile robots:

u θ = v tan δ F L , (1)
with L the car-like mobile robot's wheelbase. Using these notations, the objective is to ensure the convergence and the stabilization of the lateral deviation y to some desired distance y d , whatever the robot's speed. To achieve this goal, the robot's state vector is defined as X = [s y θ ] T .

B. Motion equations of mobile robots

In order to derive the control laws, the kinematic equations of motion with respect to the edge Γ to be followed are considered (dynamical phenomena are assumed to be negligible). Referring for instance to [START_REF] Samson | Modeling and control of wheeled mobile robots[END_REF], the robot model is then:

Ẋ =         ṡ = v cos ( θ ) 1 -y c(s) ẏ = v sin( θ ) θ = u θ - c(s) v cos ( θ ) 1 -y c(s)         (2)
Lines 1 and 3 in model (2) exist under the assumption that y = 1 c(s) (control point O is not superimposed with the center of curvature of the edge to be followed), which is considered to be always satisfied. This model expresses the derivative of the state vector with respect to time. As a consequence, when control laws are derived from such a model, their tuning impose a settling time, whatever the speed of the robot is. Nevertheless, when considering edge following, it may be more relevant to impose a settling distance, in order to control explicitly the robot trajectory. To this aim, the derivative with respect to the traveled distance s is introduced: let us define by ξ = dξ ds the derivative of any variable ξ with respect to the traveled distance. Model (2) can then be rewritten as:

X =       s = 1 y = (1 -y c(s)) tan ( θ ) θ = u θ 1 -y c(s) v cos( θ ) -c(s)       . ( 3
)
This model supposes that v cos( θ ) = 0, i.e. the velocity is not null and the robot orientation is not perpendicular to the edge to be followed. At this step, one can consider that these assumptions are met. Nevertheless, it will be shown in the sequel that these conditions are not restrictive for the proposed control law design.

III. BACKSTEPPING APPROACH

A generic control framework has been proposed in [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF] for the problem of edge following for different classes of mobile robots. It consists in a backstepping control approach, illustrated in Figure 2, composed of the following steps:

• Computation of a desired angular deviation θd . The first step consists in deriving, thanks to model (3), a desired orientation for the robot allowing to reach the desired lateral deviation within a settling distance. Specifically, a first order asymptotically stable differential equation is imposed to the lateral error e(s) = y(s)y d (s) by designing θd as:

θd (e y (s)) = arctan k y e y (s) α

where α = 1y c(s) and k y is a negative gain tuning the settling distance.

• Angular deviation regulation.

The robot will converge to the desired lateral deviation if its own orientation is equal to the desired one (i.e. θ = θd ). The second step is then devoted to the regulation of the angular velocity u θ in order to ensure such a convergence. From model (2), a first order asymptotically stable differential equation is imposed to the angular error e θ (s) = tan( θ (s))tan( θ d (s)) by designing u θ as:

u θ (s) = k θ e θ (s) cos 2 θ + c(s) v cos θ α , (5)
with k θ a negative gain tuning the settling time to reach the desired orientation. Once an expression for u θ has been obtained, a skid-steered mobile robot can be directly controlled. For a car-like mobile robot, the CONFIDENTIAL. Limited circulation. For review only.

Manuscript 25 submitted to 2020 4th IEEE Conference on Control Technology and Applications (CCTA). Received January 15, 2020.

corresponding steering angle δ F can then easily be deduced from (1). As shown in [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF], this two-step control architecture imposes the convergence of the lateral deviation y to a desired set point y d , whatever the robot's speed as long as it is strictly positive. Fig. 2. Backstepping Control Architecture used in [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF] for edge following Once the lateral dynamic is regulated thanks to the previous strategy, the speed regulation can be addressed. In the case of skid-steered mobile robots, the yaw rate can be defined independently from the longitudinal speed (the robot can turn on itself for instance), so there is no constraint on velocity design. In order to enhance the capability to achieve maneuvers, the velocity is decreased when large variations of orientation have to be achieved. Specifically, the proposed velocity control law is:

v = f (e θ ) ( 6
)
where f is the function depicted in Figure 3: the robot speed is decreased when the difference between actual and desired orientation increases, so that skid-steered mobile robots can follow accurately a tricky and narrow edge. In contrast, for car-like mobile robots, the yaw rate also relies on the longitudinal velocity, so that a velocity regulation cannot be used to enhance maneuverability. Moreover, the control architecture proposed in [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF] and summarized in this section is only reactive. As a result, harsh variations of curvature imply harsh variations of the desired orientation. This leads to abrupt reductions of the speed for skid-steered mobile robots or an impossibility to follow the edge for carlike mobile robots (or at least, an important overshoot due to the actuators properties and the settling distance imposed by the gain k y in equation ( 4)). To address this drawback, this paper proposes to take part of the knowledge of the edge to be followed in order to anticipate for such variations thanks to a model predictive control approach.

IV. MODEL PREDICTIVE CONTROL LAYER

A. Limitation of feedforward approach

A natural way to anticipate for harsh variations of curvature in the framework of edge following is to consider the lateral and angular deviations of the robot ahead of its actual position. In other words, instead of using y(s) and θ (s), one can consider these variables at a curvilinear abscissa s H in front of the robot, namely y(s + s H ) and θ (s + s H ) can be considered as the inputs of control laws (4) and [START_REF] Samson | Modeling and control of wheeled mobile robots[END_REF]. In this feedforward point of view, the control laws can then be expressed as:

     u θ (s + s H ) = k θ e θ (s + s h ) cos 2 θ + c(s + s H ) v cos θ α θd (e y (s)) = arctan
k y e y (s + s H ) α (7) As it will be pointed out in the sequel, such an approach permits to anticipate for curvature variations, but is very sensitive to the choice of the horizon of prediction s H . Moreover, because of the noise attached to the edge detection and estimation, the predicted heading may be very noisy and the stability of such a feedforward approach becomes questionable. In this paper, a formal approach is favored, in order to derive a rigorous predictive control law, allowing the anticipation of harsh variations, while preserving the control stability and the accuracy of the edge following.

B. Predictive desired orientation computation

In this paper, the objective is to control accurately a mobile robot with respect to a detected edge, potentially containing harsh variations of curvature. As mentioned previously, a reactive approach does not permit to achieve such a goal satisfactorily, since it leads to:

• important reductions of the velocity in the case of skid-steered robots, in order to let the robot's heading converge to the desired one. • important overshoots in the case of car-like robots, because of settling distances and actuator limitation, potentially leading to a collision with the edge. Since the edge is detected at each iteration, one can address these drawbacks by considering predictive control. In this point of view, the objective is to derive the desired angular deviation in order to minimize the lateral error along a horizon of prediction s H . Previously, a reactive desired orientation θd has been computed with respect to the current error in order to impose its convergence within a settling distance tuned by k θ . In this paper, the objective is to find the desired heading which minimizes the predicted lateral error along s H . For that purpose, one can consider the criterion J defined by:

J = Σ n H k=0 e k y -e y (s)e -λ k.S 2 (10)
where e y (s)e -λ k.S constitutes a desired shape for the convergence of e y (s), which can be tuned by positive parameter λ . Let us regard ξ = α tan θ as a virtual control variable related to the orientation. The optimal control ξ H = α tan( θ H d ) is then the one which minimizes the criterion J. It can therefore be inferred from the condition:

∂ J ∂ ξ = 0 (11)
The formal computation of condition [START_REF] Tourrette | Tracking footprints for agricultural applications: a low cost lidar approach[END_REF] permits to find an explicit expression for the optimal desired angular deviation:

ξ H = -1 σ 2 [ e y (s) σ 1 -a 1 σ 2 -a 2 σ 3 -e y (s) σ e ] θ H d = -arctan 1 α σ 2 [ e y (s) σ 1 -a 1 σ 2 -a 2 σ 3 -e y (s) σ e ] (12)
with the following notations:

         σ 1 = Σ n H k=0 (k S) σ 2 = Σ n H k=0 (k S) 2 σ 3 = Σ n H k=0 (k S) 3 σ e = Σ n H k=0 (k S e -λ k.S ) (13)
The angular deviation θ H d , since it satisfies condition [START_REF] Tourrette | Tracking footprints for agricultural applications: a low cost lidar approach[END_REF], is the optimal angular deviation minimizing the criterion J along s H . It can be considered as a predictive term which has to be applied at present time to anticipate for edge variations. It can then be substituted to the reactive term θd previously computed in (4) in the first step of the backstepping control approach.

C. Yaw rate servoing

Once the predictive desired orientation ( 12) is applied, the robot minimizes the criterion J defined by [START_REF] Bayar | Localization and control of an autonomous orchard vehicle[END_REF] while converging to the desired lateral deviation, i.e. its lateral error is minimized along the horizon s H . The second step then remains unchanged and aims at ensuring the convergence of θ to the desired predicted one, i.e. θ H d . Similarly to [START_REF] Deremetz | A generic control framework for mobile robots edge following[END_REF], this can be achieved by imposing:

u H θ = k θ e θ cos 2 θ + c(s) v cos θ α . ( 14
)
with e θ = tan( θ )tan( θ H d ), since this indeed leads to the following differential equation:

ė θ = k θ e θ ( 15
)
and ensures the convergence θ → θ H d as soon as k θ is properly set. It can be directly applied as the yaw rate control in the case of skid-steered robots.

For car-like mobile robots, reporting ( 14) into (1) leads to:

δ F = arctan L u H θ v (16)
Control expression ( 16) thus obtained requires that velocity v is non-null. In practice, this assumption is met, since the velocity is supposed to be strictly positive in order to achieve edge following. Nevertheless, this assumption can be relaxed by considering the derivative of the yaw rate of a carlike mobile robot with respect to the curvilinear abscissa. Reporting (1) into (3) leads to:

d θ ds = α tan δ F L cos θ -c(s) (17)
Then, it can be shown that the following control expression:

δ F = arctan L cos 3 ( θ ) α k θ e θ + c(s) cos 2 ( θ ) , (18)
leads to the first order asymptotically stable differential equation:

e θ = k θ e θ , (19)
This ensures the convergence θ → θ H d within a settling distance imposed by the parameter k θ .

V. EXPERIMENTAL RESULTS

A. Experimental set-up

The proposed predictive control law has been tested on two robots in two different conditions 1 , as depicted in Figure 4: • a skid-steered mobile robot. It moves in indoor conditions and follows several kinds of structures, such as depicted in Figure 4(a). It weights 25kg, with a wheelbase of 0.6m. The Lidar is settled at 0.8m from the rear wheel and with an elevation of 0.5m; • a car-like mobile robot moving in off-road conditions, following a row of trees, see Figure 4(b). This robot has a weight of 520kg and a wheelbase of 1.2m. The Lidar is in front of the robot at a 0.65m height and is situated at 1.8m from the middle of the rear axle. Both robots are equipped with a Sick-LMS laser. The structures to be followed, quite different, are depicted by red lines in Figure 4, reflecting some harsh curvatures in both conditions. Indoor, the robot must first follow an artificial 1 Video of experiments are available at https://stratus.irstea.fr/f/acac69b39e0e4c8abe77/ In both cases, the edge is detected as a cloud of points by the Lidar. From this cloud, points are sorted in groups depending on their distance to the robot. From these groups, one can identify a second order curve using the closest points to the robot. Potentially, two curves have to be defined, one on the left and one on the right. Depending on the distance, only one is selected. For the sake of clarity, in the forthcoming experiments, only an edge on the right of the robot is detected in the indoor case, while only an edge on the left is detected in the outdoor application. The parameters detailed in the following table are used during the experiments. Selected parameters are different for the two robots, since their properties are slightly different. In particular, for reactive and feedforward approaches (where k y is used), the settling distance is much lower for the skid-steered robot, since it is able to turn on itself and may then have a reduced settling distance. One can note that in Table I, the parameter k y is only used for reactive control (i.e. first step of the backstepping control law (4)) and feedforward control [START_REF] Gonzalez-De Soto | Autonomous systems for precise sprayingevaluation of a robotised patch sprayer[END_REF]. Similarly, the horizon of prediction s H , as well as the parameter λ , are only used in the framework of the proposed predictive control [START_REF] D'andréa-Novel | Control of nonholonomic wheeled mobile robots by state feedback linearization[END_REF].

B. Skid-steered results

Let us first consider the results obtained in indoor conditions with the skid-steered mobile robot. The structure to be followed, see Figure 4(a), is composed of three phases:

• A: the following of an artificial hedge (from 0 to 4m); • B: the following of construction poles in front of the artificial hedge (from 4 to 7m); • C: the following of boxes arranged in a straight line perpendicular to the artificial hedge (after 8.5m). This structure has been followed three times and a comparison of lateral errors obtained during these tests is reported in Figure 5:

• the first test uses the reactive backstepping approach described by ( 4) and ( 5). It is reported in blue dotted line; • the second test is based on the feedforward method described by [START_REF] Gonzalez-De Soto | Autonomous systems for precise sprayingevaluation of a robotised patch sprayer[END_REF] and results are reported in red dashed line; • finally, the third run has been achieved using the proposed model predictive control ( 12) and ( 14). The results are reported in plain magenta line. The results depicted in Figure 5 show the improved accuracy obtained thanks to the proposed model predictive control law. The error converges around 0 during phase A (i.e up to 4m) and phase B (from 4 to 7m) as well. It is not exactly the case for the other approaches. The feedforward strategy is less accurate, since the relevance of the predicted error strongly depends on the quality of the edge detection: the heading of a skid-steered robot may be fast varying. In these situations, the detected edge also changes quickly and the estimated error at s + s H is then noisy. It results in some instabilities and a larger convergence distance. Eventually, the reactive approach converges to zero, but is affected by the step in the trajectory imposed by the construction poles. As it does not anticipate for this step, an error of 25cm is recorded at that point. The transition between phases B and C (from 7 to 8.5m) generates an equivalent overshoot (30cm). At this transition, the robot has to turn on itself: the edge detection is then changing quickly and the robot cannot adapt instantaneously its position. After the robot is correctly oriented, it can speed up and then converge to zero with the three approaches, even if model predictive control is still slightly more accurate. This follows from the anticipation of edge variations, allowing the robot to stay closer to the set point related to the detected edge.

This anticipation also permits to smooth the control variables and then limit the speed reduction when large variations of the desired orientation θd are required. This can be viewed in Figure 6. One can indeed see in Figure 6 that both control laws using anticipation (feedforward and model predictive control) are able to foresee the big curve imposed by the transition between phases B and C (around 8m). Even if the velocity is reduced, it does not go to zero. On the contrary, reactive control does not anticipate for the curve variation and an abrupt variation of heading is then required. As a result, the robot has to stop in order to decrease the heading error. This explains the important error of -20cm recorded using this control law at abscissa 9.2m, see Figure 5. For skid-steered mobile robots, the proposed predictive control law permits to ensure a high level of accuracy thanks to the anticipation of edge variations and enhances as well the stability of the tracking.

C. Car-like mobile robot results

Let us now consider the results related to the car-like mobile robot in the outdoor environment shown in Figure 4(b). In this configuration, the robot cannot turn on itself (i.e. turning without a longitudinal velocity). Since the target velocity is 0.8m.s -1 , the harsh curve is difficult to track due to the steering angle saturation. Reactive control (4) is indeed unable to follow the edge and tends to collide the construction barrier. This can be noticed in Figure 7, where the lateral error e y (s) for this test is reported in blue dotted line: during the tree following (up to curvilinear abscissa 13m), the robot converges to the desired distance y d . Nevertheless, after 13m the robot does not turn sufficiently to maintain the error close to zero and the error rises up to 1.5m. The robot then stops as it is about to collide the barrier.

On the contrary, the predictive control proposed in this paper succeeds in preserving a reduced overshoot when the robot is in the harsh curve imposed by the barrier: in Figure 7, the result related to the predictive control is reported in magenta plain line. Before 13m, one can see the convergence of the error to zero and next a transient overshoot of 0.4m is recorded when the curve appears. This proposed approach permits to achieve harsh maneuvers, without colliding the edge to be followed. In order to highlight the benefits of the proposed predictive control with respect to a feedforward approach, the results related to this latter control [START_REF] Gonzalez-De Soto | Autonomous systems for precise sprayingevaluation of a robotised patch sprayer[END_REF] are reported in red dashed line in the same figure. Just before the curve (at 13m), the robot starts turning in advance and consequently an overshoot inside the curve appears: the robot moves aside from -0.5m.

Moreover, during straight line parts, the feedforward anticipation of the robot orientation generates an oscillating behavior and then an inaccurate tracking. This can be checked in Figure 8, where the same control has been used to achieve tree following when the construction barriers are removed. The robot has then to follow an almost straight line, even if the irregularities of the trees generate some noise. One can see that the feedforward approach is not as stable as both reactive and predictive approaches, leading to errors up to 0.6m when the reactive and predictive controls stay accurately around a null error, with a more stable behavior for the proposed predictive approach. Table II shows the mean and standard deviation of the lateral error e y (s) during these last three tests and highlights the very satisfactory results achieved by the proposed predictive control strategy in this straight line following scenario. Its accuracy in harsh curve following had been demonstrated in Figure 7 above.

VI. CONCLUSIONS

In this paper, a generic predictive control framework is proposed for the problem of edge following, here detected thanks to a lidar. It permits to ensure a stable and high accurate tracking by anticipating for huge variations of curvature. This enables to follow an important collection of structures, even the structures difficult to track because of robot limitations (in particular, the radius of curvature of car-like robots). The genericity of the proposed algorithm gives the possibility to consider in the same way the control of different kinds of robots (skid-steered, car-like, etc.). In the proposed approach, two steps are considered to ensure the accurate convergence of the robot at the desired distance from the edge to follow. The first step is dedicated to the computation of the desired heading for the robot in order to converge to the targeted lateral deviation along a distance of prediction. This desired heading is computed thanks to a model predictive control approach, as it minimizes the lateral error along the detected edge, locally approximated as a second order curve. Once the optimal desired heading is computed, the yaw rate is controlled. A convergence with respect to time is considered (particularly interesting for a skid-steered robot, since it can turn on itself) or with respect to the traveling distance (particularly interesting for car-like mobile robots, since it permits to avoid a singularity at null speed). Such an approach permits to separate the lateral and longitudinal dynamics, allowing to control independently the robot's speed and then modulate the velocity to enhance turning capabilities. The efficiency of the proposed control law has been highlighted through full-scale experiments with two types of robots and in two different configurations (using different kinds of structures and environments). The results demonstrate the capability of anticipation and the accuracy that can be reached thanks to the model predictive approach proposed in this paper. If variations in the edge to be followed are an important perturbation, the terrain conditions are also to be accounted for edge following at higher speeds. Future works are focused on the observation of grip conditions, compatible with only the edge detection and no additive exteroceptive sensors, in order to preserve the accuracy at high speed when moving on uneven surfaces with car-like mobile robots.

Fig. 1 .

 1 Fig. 1. Modeling of a generic mobile robot in Frénet frame

Fig. 3 .

 3 Fig. 3.Dependency between the linear speed and the robot angular deviation

First

  , let us express the lateral error e y (s + ∆s), with ∆s ∈ [ 0 ; s H ]. Considering Figure 1 and robot models (2) CONFIDENTIAL. Limited circulation. For review only. Manuscript 25 submitted to 2020 4th IEEE Conference on Control Technology and Applications (CCTA). Received January 15, 2020.

and ( 3 )

 3 , one can write: e y (s + ∆s) = y(s) + α tan θ ∆sy d (s)a 1 ∆sa 2 ∆s 2 , (8) where a i , i ∈ 1, 2 are the coefficients of the local second order interpolation of the detected edge. In order to achieve the minimization, let us consider n H ∈ N points of coincidence, regularly spaced along [ 0 ; s H ] by a step S = s H n H . The error e k y at the k th coincidence point (k ∈ [ 0 ; n H ]) is obtained straightforwardly by injecting ∆s = k.S into (8):e k y = y(s) + α tan θ (k.S)y d (s)a 1 (k.S)a 2 (k.S) 2 = e y (s) + α tan θ (k.S)a 1 (k.S)a 2 (k.S) 2 (9) Using this formulation, the objective of predictive control is to find an expression for θ which minimizes the lateral error e k y for k ∈ [ 0 ; n H ].

Fig. 4 .

 4 Fig. 4. Robots used during experiments

Fig. 5 .

 5 Fig. 5. Comparison of lateral errors in indoor conditions using different control strategies

Fig. 6 .

 6 Fig. 6. Comparison of robot velocity in indoor conditions

Fig. 7 .

 7 Fig. 7. Comparison of lateral errors in outdoor conditions using different control strategies

Fig. 8 .

 8 Fig. 8. Comparison of lateral errors during trees following

  CONFIDENTIAL. Limited circulation. For review only.Manuscript 25 submitted to 2020 4th IEEE Conference on Control Technology and Applications (CCTA). Received January 15, 2020.

Manuscript 25 submitted to 2020 4th IEEE Conference on Control Technology and Applications (CCTA). Received January 15, 2020.

ACKNOWLEDGMENT

This work has been sponsored by the French National Research Agency under the grant number ANR-14-CE27-0004 attributed to Adap2E project (adap2e.irstea.fr). It has also been sponsored by the French government research program "Investissements d'Avenir" through the IDEX-ISITE initiative 16-IDEX-0001 (CAP 20-25), the IMobS3 Laboratory of Excellence (ANR-10-LABX-16-01) and the RobotEx Equipment of Excellence (ANR-10-EQPX-44). This research was also financed by the European Union through the Regional Competitiveness and Employment program -2014-2020-(ERDF -AURA region) and by the AURA region.