Dr Wei Hu

Zhen-Yu Yin

Gianvito Scaringi

Christophe Dano

Pierre-Yves Hicher

Relating fragmentation, plastic work and critical state in crushable rock clasts

Keywords: particle crushing, fragmentation, critical state, plastic work, relative breakage

come

INTRODUCTION

Clasts fragmentation (also termed grain breakage) can occur during compression and shearing of granular materials, especially under high confining stresses and dynamic loadings [START_REF] Coop | Particle breakage during shearing of a carbonate sand[END_REF][START_REF] Kaproth | Deformation band formation and strength evolution of unlithified sand: The role of grain breakage[END_REF][START_REF] Coop | Changes to particle characteristics associated with the compression of sands[END_REF][START_REF] Kimura | The Effect of Effective Normal Stress on Particle Breakage, Porosity and Permeability of Sand: Evaluation of Faults around Methane Hydrate Reservoirs[END_REF][START_REF] Kimura | Shear-induced permeability reduction and shear-zone development of sand under high vertical stress[END_REF][START_REF] Mehta | An investigation on the particle breakage of Indian River sands[END_REF]. Important breakage can also occur under low confining stress, particularly if the grains are fragile and present irregular shape (e.g., [START_REF] Nakata | A probabilistic approach to sand particle crushing in the triaxial test[END_REF]. Grain breakage is a fundamental issue for the stability of earth dams, embankments and railway ballasts, and for the bearing capacity of piles (e.g., Daouadji et al., 2011;[START_REF] Okonta | Effect of grading category on the roundness of degraded and abraded railway quartzites[END_REF][START_REF] Winter | Influences of particle characteristic and compaction degree on the shear response of clinker ash[END_REF]. Quantitative information about the shearing history can be retrieved from glacial sediments [START_REF] Altuhafi | A hypothesis on the relative roles of crushing and abrasion in the mechanical genesis of a glacial sediment[END_REF] and landslide granular deposits [START_REF] Zhang | Rock-avalanche deposits store quantitative evidence on internal shear during runout[END_REF], and can be crucial for understanding the likelihood of the latter to exhibit further instability [START_REF] Hu | Sensitivity of the initiation and runout of flowslides in loose granular deposits to the content of small particles : An ins ight from flume tests[END_REF](Hu et al. , 2018a(Hu et al. , 2018b(Hu et al. , 2018c;;[START_REF] Fan | What we have learned from the 2008 Wenchuan Earthquake and its aftermath: A decade of research and challenges Engineering Geology[END_REF]. Extensive fragmentation has been also related to frictional weakening in rock avalanches, although the issue is still debated [START_REF] Davies | A fragmentation-spreading model for long-runout rock avalanches[END_REF][START_REF] Fan | Failure mechanism and kinematics of the deadly June 24th 2017 Xinmo landslide, Maoxian, Sichuan, China[END_REF][START_REF] Scaringi | Shear-Rate-Dependent Behavior of Clayey Bimaterial Interfaces at Landslide Stress Levels[END_REF][START_REF] Scaringi | Discussion on: "Experimental study of residual strength and the index of shear strength characteristics of clay soil[END_REF][START_REF] Ren | The grain size distribution and composition of the Touzhai rock avalanche deposit in Yunnan, China[END_REF] and further, dedicated experimental investigations are required.

The amount of fragmentation can be quantified through the changes of grain size distribution (GSD). Different ways to measure its evolution have been suggested (e.g., [START_REF] Marsal | Large-scale testing of rockfills materials[END_REF][START_REF] Hardin | Crushing of soil particles[END_REF][START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF]. The effect of particle crushing on various physical and mechanical properties has been investigated (e.g., [START_REF] Mcdowell | Fractal crushing of granular materials[END_REF][START_REF] Bandini | The influence of particle breakage on the location of the critical state line of sands[END_REF][START_REF] Zhang | Particle shape factors and fractal dimension after large shear strains in carbonate sand[END_REF], and introduced into constitutive and numerical modelling (e.g., [START_REF] Ma | Modeling the fragmentation of rock grains using computed tomography and combined FDEM[END_REF][START_REF] Zhou | Discrete element modeling of crushable sands considering realistic particle shape effect[END_REF].

However, understanding how the breakage amount can be measured and how this measure can be related to mechanical properties explicitly is still challenging. The three following issues need to be investigated: [START_REF] Hu | Internal erosion controls failure and runout of loose granular deposits: Evidence from flume tests and implications for post-seismic slope healing[END_REF] what are suitable measures for the amount of grain breakage; (2) how to determine these measures based on stress-strain histories; and (3) how to relate these measures to mechanical properties of the granular assembly. Once these three points are solved, the quantity measuring the amount of grain breakage can be introduced as a variable to control the influence of particle crushing on the mechanical properties of a granular material undergoing a generic

ACCEPTED MANUSCRIPT

stress-strain path.

For the first issue, a breakage measure varying in the range 0-1 with the GSD changing seems convenient. Among the proposed measures, the modified relative breakage index B r * of [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF], and the relative uniformity B u , based on the fractality of the GSD (e.g., [START_REF] Coop | Particle breakage during shearing of a carbonate sand[END_REF] and on the coefficient of uniformity C u , seem good choices. Various methods to quantify the GSD change under different loading conditions were developed and validated through tests under monotonic loading: i.e. the crushing surface approach [START_REF] Kikumoto | Particle crushing and deformation behaviour[END_REF], and the energy approach [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF]. So far, these two approaches have not been fully examined under different loading conditions, including under cyclic loading. Therefore, in this paper we propose a unified approach to quantify the change of GSD (in terms of B u and B r * ) for different loading conditions.

Changes of the GSD can be correlated to the evolution of mechanical characteristics. Several properties (stress-dilatancy, peak shear strength, plastic modulus, etc.) depend on the distance of the current stress state (p', e) from the corresponding critical state (p', e c ) of the granular material.

Hence, understanding how fragmentation-induced changes of GSD affect the location of the critical state line (CSL) becomes crucial for relating particle crushing to the mechanical behaviour. CSL-GSD relationships have been demonstrated by experiments on crushable sands [START_REF] Bandini | The influence of particle breakage on the location of the critical state line of sands[END_REF] and by discrete element simulations on crushable granular materials [START_REF] Bolton | Micro-and macro-mechanical behaviour of DEM crushable materials[END_REF].

In the following sections, we first analyse the measures of particle crushing amount. We carry out laboratory tests on crushable limestone grain assemblies under various stress-strain paths, with both monotonic and cyclic loadings to propose a unified determination method for the measure of particle crushing. Secondly, we use the test results to propose a relation between these measures and the CSL location. Finally, we apply the relations obtained on limestone fragments to other materials to validate the method.

MEASURES OF GRAIN BREAKAGE

Different measures have been proposed, such as: (1) B 15 =D 15i /D 15f , where D 15 is the grain diameter for a percentage of passing equal to 15% in weight, and subscript "i" and "f" represent the initial and final values, e.g. before and after testing [START_REF] Lee | Compressibility and crushing of granular soil in anisotropic triaxial compression[END_REF]Fig. 1a); (2) R/100, that is the increase in % passing the sieve having the largest increase, measured from the GSDs before and after testing [START_REF] Marsal | Large-scale testing of rockfills materials[END_REF]Fig. 1b); (3) the relative breakage B r =B t /B p , where B p is the area between the line d = 0.074 mm and the portion of the initial GSD with d > 0.074 mm, and B t is the area between the initial GSD and the GSD after loading with d > 0.074 mm [START_REF] Hardin | Crushing of soil particles[END_REF]Fig. 1c); (4) the particle breakage factor B 10 =1-D 10f /D 10i [START_REF] Lade | Significance of particle crushing in granular materials[END_REF]Fig. 1d); (5) the change of the coefficient of uniformity C u =D 60 /D 10 , related to the amount of grain breakage [START_REF] Biarez | Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires[END_REF]Fig. 1e); ( 6) the particle breakage factor B f =R/100, where R is the percentage of particles smaller, after testing, than the smallest particle size in the initial GSD [START_REF] Nakata | A probabilistic approach to sand particle crushing in the triaxial test[END_REF]Fig. 1f) [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF]Fig. 1g); (8) the grading state index I G =B t /B p , with the initial grading curve determined by lower limit of the biggest particle size [START_REF] Kikumoto | Particle crushing and deformation behaviour[END_REF]Fig. 1h). As the grading proceeds toward a fractal distribution, only Einav's measure B r * evolves from 0 to 1 following the evolution of the GSD, therefore it is selected in this study to describe the breakage amount.

Another measure can be defined, that can be named "relative uniformity" B u and employs the coefficient of uniformity C u =D 60 /D 10 in the following way:

u ui u uf ui CC B CC    (1)
where C ui is the initial coefficient of uniformity, and C uf is the ultimate coefficient of uniformity proper of the final, fractal grading. Hence, B u also evolves from 0 to 1. Note that B u includes also the change of D 10 which is a factor influencing the mechanical properties as shown by [START_REF] Lade | Significance of particle crushing in granular materials[END_REF]. Since the definition of B u is based on C u , which is a common coefficient in engineering practice and can be easily measured, the use of B u is convenient to describe the grain breakage amount, provided that it can be related to mechanical variables influencing grain breakage.

Therefore, in addition to B * r , B u is also examined in this study in connection with the amount of grain breakage.

A UNIFIED APPROACH

As grain breakage is produced by mechanical loading, the breakage amount can be related to the input energy. Based on this idea, [START_REF] Lade | Significance of particle crushing in granular materials[END_REF] proposed a hyperbolic relation between the particle breakage factor B 10 and the total energy input calculated by stresses and total strains. [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF] proposed to link the breakage energy, using stresses and total strains, to the relative breakage B r * for compression tests. Different from the energy methods, [START_REF] Simonini | Analys is of Behavior of Sand Surrounding Pile Tips[END_REF] and [START_REF] Kikumoto | Particle crushing and deformation behaviour[END_REF] defined a crushing surface to determine the breakage amount based on the distance from the stress state to this crushing surface. However, the accumulation of grain breakage during cyclic loading cannot be described by these different approaches.

Another measure in relationship with the amount of crushing is the plastic work of [START_REF] Daouadji | An elastoplastic model for granular materials taking into account grain breakage[END_REF]. To avoid the non-uniqueness of the summation of the plastic work along a loading path where stresses or strains are reversed, we adopted the following different definition of the plastic work:

p p w d     (2)
where  and d p are stress and incremental plastic strain tensors, respectively, and the use of

Macaulay's brackets implies that   2 F F F  .
In this way, under cyclic loading, the plastic work can be cumulated [START_REF] Hu | Influence of grain breakage on critical state[END_REF]. Thus, linking the plastic work to the measure of breakage amount makes it possible to take the additional grain breakage after the first loading into account.

MATERIALS AND METHODS

In order to establish the relation between the measures of breakage ( B u and B r * ) and the plastic work, we performed drained triaxial tests under both monotonic and cyclic conditions.

For each test, the GSD was measured after testing and the plastic work was calculated along the loading path.

We first employed crushable limestone fragments from a quarry in Prefontaines, France. The

ACCEPTED MANUSCRIPT

shape of the grains is sub-angular according to the D2488-10 ASTM standard. [START_REF] Head | Manual of soil laboratory testing[END_REF].

The relations between breakage measures, plastic work and CSL location obtained on the limestone grains are also applied to granite fragments, Cambria sand (e.g., [START_REF] Yamamuro | Drained sand behavior in axisymmetric tests at high pressures[END_REF] and Sydney sand [START_REF] Russell | A bounding surface plasticity model for sands exhibiting particle c rushing[END_REF]), as will be shown in Section 5.3.

The granite fragments consist of sub-angular grains of quartz with feldspars, with sizes between 0.3 mm and 2 mm. The material has an initial uniformity coefficient C u = 2.0, a median grain size D 50 = 0.85 mm and a specific gravity G s = 2.65. Minimum and maximum void ratios are 0.83 and 1.23. The initial e ref0 = 1.28 and  c = 0.088 were obtained adopting the correlation by [START_REF] Biarez | Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires[END_REF]. The relative uniformity, modified relative breakage index, CSL location and applied plastic work were obtained from drained triaxial tests carried out by [START_REF] Kim | Etude expérimentale du comportement mécanique des matériaux granulaires sous fortes contraintes[END_REF] with confining stress of 1, 5, 10, 30 and 60 MPa.

The Cambria sand is composed of sub-rounded grains of quartz (with small quantity of gypsum) with sizes between 0.83 mm and 2 mm, a median grain size D 50 = 1.43 mm and an initial uniformity coefficient C u = 1.3. Minimum and maximum void ratios are 0.503 and 0.792.

The specific gravity is G s = 2.69. The initial e ref0 = 0.73 and  c = 0.088 were measured from drained and undrained triaxial tests under low confining pressures (100, 200, 500 kPa) with no change of grading [START_REF] Yamamuro | Drained sand behavior in axisymmetric tests at high pressures[END_REF], thus Biarez and Hicher's correlation is not needed. The relative uniformity, modified relative breakage index, CSL location and applied plastic work were obtained using drained and undrained triaxial tests under high confining

ACCEPTED MANUSCRIPT

pressures up to 52 MPa [START_REF] Lade | Relative density effects on drained sand behaviour at high pressures[END_REF][START_REF] Bopp | Relative density effects on undrained sand behavior at high pressures[END_REF][START_REF] Yamamuro | Drained sand behavior in axisymmetric tests at high pressures[END_REF].

The Sydney sand is mostly made of quartz. It was sampled from sand dunes at Kurnell, Sydney. It has a median grain size D 50 = 0.31 mm, initial uniformity coefficient C u = 1.83, specific gravity G s = 2.65. Minimum and maximum void ratios are 0.60 and 0.92. The initial e ref0 = 0.97 and  c = 0.07 were obtained from drained and undrained triaxial tests under low confining pressures (50, 150 and 200 kPa) with no change of grading [START_REF] Russell | A bounding surface plasticity model for sands exhibiting particle c rushing[END_REF]. The relative uniformity, modified relative breakage index, CSL location and applied plastic work were obtained using drained and undrained triaxial tests under high confining pressures up to 7.8

MPa [START_REF] Russell | A bounding surface plasticity model for sands exhibiting particle c rushing[END_REF].

Photographs of the limestone material before testing and after a drained triaxial test under constant mean effective stress p' = 400 kPa are shown in Fig. 2. By comparing Fig. 2a and Fig. 2b, one can see that some grains were crushed into smaller pieces, and some grains were damaged with visible cracks, while some others appear still intact. The breakage amount can be estimated by the GSD measured after testing. Figure 3a presents the results of two drained triaxial tests with confining stresses p' 0 = 50 kPa and 100 kPa. The elastic properties can be measured on the stress-strain curves. We assume a non-linear elasticity and we adopt the formulation proposed by [START_REF] Richart | Vibration of Soils and Foundations[END_REF]:

  0.5 2 0 2.97 1 at at e p G G p ep        (3)   0.5 2 0 2.97 1 at at e p K K p ep        (4)
where p at = 101.325 kPa is the atmosphere pressure.

The shear modulus was measured on the stress-strain curves up to a deviatoric strain of 0.1% in Fig. 3a. We found G 0 = 80. The bulk modulus was measured from the isotropic compression tests (Fig. 3b). We found K 0 = 105, which corresponds to a Poisson's ratio of 0.2. Taking into consideration the decrease of the shear modulus with the strain level found for various sands (e.g., [START_REF] Clayton | Stiffness at small strain: research and practice[END_REF], we can assume that the elastic modulus at very small strains is

ACCEPTED MANUSCRIPT

approximately five times bigger than the elastic modulus at 0.1%. Therefore, G 0 = 400 and K 0 = 525 corresponding to the very small strain level were used to estimate the elastic strain component.

According to [START_REF] Mcdowell | Fractal crushing of granular materials[END_REF], the GSD for both initial and tested samples can be approximately expressed as follows

  3 ( ) / M F D D d    (5)
where F(D) is the percentage of passing corresponding to the grain size D. In agreement with [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF], a fractal GSD was assumed with  = 2.6 according to [START_REF] Coop | Particle breakage during shearing of a carbonate sand[END_REF].

Therefore, the relative breakage index B r * for each test can be obtained as proposed by [START_REF] Einav | Breakage mechanics-Part I: Theory[END_REF]. For the measure of B u , the current C u and the ultimate C uf are needed besides C ui . C u can be obtained from the grading curve at the end of each test. C uf corresponding to fractal GSD can be determined using Eq. ( 5) with  = 2.6 as follows: for F(D) = 60% and 10%, 60% = (D 60 /d M ) 3- and 10% = (D 10 /d M ) 3-2.6 , then C uf = D 60 /D 10 = 6 1/(3-) = 88.18 is obtained. Therefore, based on Eq. ( 1), the value of B u can be obtained once C u is measured after the test.

The corrections of the volumetric change due to the membrane penetration is important for accuracy of analysis. We used the method proposed by Baldi and Nova (1984) to correct the test results. Baldi and Nova (1984) investigated the membrane penetration in triaxial tests and found that m V  depended strongly on D 50 , ' 3  , and the diameter of the specimen (D s ) as well as the membrane characteristics. The following equation was used to account for the membrane penetration effect:

1 ' 3 50 3 50 0 2 m s m m D D V V D E t                     (6)
where m E is the Young's modulus of the membrane material and m t is the thickness of the membrane.

ACCEPTED MANUSCRIPT

RESULTS

TESTS UNDER MONOTONIC LOADING

Tests at different confining stresses

Drained triaxial tests under various constant confining stresses (constant-' 3 tests) varying from 22 kPa to 800 kPa were carried out (Fig. 4). One can notice that, for elevated confining stresses, the material appears very ductile and the stress ratio q/p' continues to increase even at large strains (Fig 4a). This is due to a large amount of particle breakage developing during loading, as shown in Fig. 4c. Figure 4b presents the volume changes during the drained triaxial tests. They show a dilative behaviour under low confining stresses of 22 kPa and 50 kPa, and a contractive behaviour under higher confining stresses of 100 kPa, 400 kPa and 800 kPa. For samples under high confining stresses, the large contraction is also a consequence of a significant grain breakage. The GSDs after testing were measured, as shown in Fig. 4c. Significant crushing occurs even under low confining stresses. The relative uniformity B u and the modified relative breakage index B r * were then obtained for each test, and the values summarised in this figure.

The results show that a higher confining stress induces a higher amount of grain breakage, thus raising the values of B u and B r * .

Tests with different loading paths

A drained triaxial test under constant mean effective stress (constant-p' test) p' = 400 kPa up to  a = 25 % was carried out, as shown in Fig. 5. The results are compared to the constant-' 3 test results at the same consolidation stress. During the constant-p' test, the lateral stress is reduced in order to keep p' constant, which also results in a smaller amplitude of the deviatoric stress. As a consequence, less grain breakage occurred compared to the constant-' 3 test (Fig. 5c). This leads to a stiffer stress-strain relationship and a smaller volume change for the constant-p' test. The evolution of the stress-strain relationship with the grading change is similar to the one obtained for the previous constant-' 3 tests. Note that the variability of the grading measurements for a given stress history is mainly due to the amount of fine particles (e.g. < 2 mm). In Fig. 5c, the difference in the gradation curves between the constant p' and the constant

ACCEPTED MANUSCRIPT

' 3 tests is statistically significant, especially for particles bigger than 2 mm. The values of B u and B r * for the constant-p' test are smaller than those for the constant-' 3 test.

Tests with different consolidation histories

A conventional drained triaxial test at ' 3 = 100 kPa up to  a = 25 % on a sample isotropically consolidated up to 800 kPa and unloaded to 100 kPa was carried out, and compared to the constant-' 3 test with ' 3 = 100 kPa without an over-consolidation history (Fig. 6).

Although the whole difference of grading between the two tests is small (Fig. 6c), the difference of distribution for particles larger than 2 mm was adequately measured. The results show the influence of stress history on the grain breakage: the test with the over-consolidation process results in slightly more grain breakage than the test without the over-consolidation history (B u and B r * for the former are higher than for the latter). This is essentially due to the isotropic loading phase up to 800 kPa where plastic strains developed (see Fig. 3b). The shearing phase did not bring any significant additional grain breakage and, therefore, the stress-strain relationship and the volume change for these two tests are in agreement with what is usually observed by comparing normally and over-consolidated material behaviours: a stiffer and less contractive behaviour for the over-consolidated sample.

CYCLIC LOADING TESTS

In order to seek a unified approach for determining the breakage amount under different loading conditions, three types of drained cyclic triaxial tests on samples isotropically consolidated at 400 kPa were carried out:

 Series 1: two tests with constant confining stress, CYC-CCS1 in which q max and q min increase from one cycle to another, CYC-CCS2 with constant q max and q min , in order to investigate the influence of the axial loading amplitude on grain breakage;

 Series 2: three tests under constant p' with different number of cycles N cycle = 0, 5, 20

(CYC-CP0, CYC-CP5 and CYC-CP20) to investigate the influence of the number of cycles on grain breakage;

ACCEPTED MANUSCRIPT

 Series 3: a test with constant deviatoric stress with varying p' (CYC-CDS) to investigate the influence of the stress ratio on grain breakage.

Tests at constant confining stress

Figure 7 shows stress-strain patterns and GSDs for two tests at constant confining stress (Series 1). The test CYC-CCS1 was strain-controlled with axial strains up to 1%, -1%, 2%, -2%, 3%, -3%, 4%, and ending by -4%. q max reached 370, 570, 770 and 890 kPa, and q min reaches to -250, -310, -310, -310 kPa. The test CYC-CCS2 was stress-controlled with the axial stress varying between q max = 865 kPa and q min = -310 kPa for 4 cycles. The stress-strain curves in Fig. 7a and the evolution of the volumetric strain in Fig. 7b can be used to calculate the plastic work during the tests. The amount of grain breakage in Test CYC-CCS1 is lower than the one in Test CYC-CCS2 (Fig. 7c), because the values of q max during the first three cycles for Test CYC-CCS1 are smaller than those applied in Test CYC-CCS2 (Fig. 7a).

Tests with constant mean effective stress

Figure 8 shows the GSDs for tests at constant p' with different number of cycles (Series 2).

For the test CYC-CP0, the sample was only monotonically sheared up to q = 730 kPa with  a = 10 % and q/p' = 1.83. The other two tests CYC-CP5 and CYC-CP20 were carried out with q max = 730 kPa and q min = -420 kPa (q/p' between 1.83 to -1.05) up to 5 and 20 cycles, respectively. The stress-strain curves in Fig. 8a and the evolution of the volumetric strains in Fig. 8b can be used to calculate the plastic work. The amount of grain breakage observed after testing increases with the applied number of cycles (Fig. 8c). Figure 8d shows the evolution of B u and B r * with the number of cycles N cycle , which demonstrates that grain breakage occurs at fast pace at the beginning of the loading, and then slows down. Based on this figure, the evolution of these two variables can be expressed by

    * * max 0 max 0 ** 0 0 B ; r r cycle u u cycle u u r r cycle cycle B B N B B N B B B a N b N         (

Tests with constant deviatoric stress

The GSD for the test at constant deviatoric stress CYC-CDS (Series 3) is shown in Fig. 9. In this test, the sample was sheared up to q = 180 kPa with  a = 0.1 % and p' = 460 kPa after an isotropic compression up to 400 kPa; then it was stress-controlled by changing the axial stress and the confining stress simultaneously, keeping constant the value of q and varying p' from 460 kPa to 140 kPa (q/p' varies from 0.39 to 1.29, less than q/p'=1.46 of the constant-' 3 test at ' 3 =400 kPa) with three cycles. Such low values of q and p' max were applied to ensure that little grain breakage would be caused by the applied stresses during the first loading. Therefore, the test can be used to investigate the influence of the stress ratio (changing between 0.39 and 1.29) on grain breakage. Note that the deviatoric stress varies from 190 to 170 kPa with an average of 180 kPa due to the lack of very accurate test control; this average value is considered for the interpretation. The stress-strain curve in Fig. 9a and the volume change in Fig. 9b show that only small strains developed during the test. The increase of the volumetric strain during the stages of decreasing p' is mainly due to elastic unloading. The amount of grain breakage observed after testing was very small (Fig. 9c), which demonstrates that the cyclic change of stress ratios with moderate amplitudes had little influence on grain breakage.

INTERPRETATION OF THE RESULTS AND DISCUSSION

RELATION BETWEEN GRAIN BREAKAGE MEASURES AND PLASTIC WORK

An energy input is needed for grains to break [START_REF] Mcdowell | Fractal crushing of granular materials[END_REF]. We decided to adopt the plastic work as a measure of this energy and to link it to the breakage measures B u and B r * representing the evolution of GSD. Under axisymmetric loading condition, the plastic work w p defined in Eq. ( 2) can be expressed as follows

d d p p p v d w p q      (8)
where d v p and d d p are the volumetric and deviatoric plastic strain increments. The plastic strain increment can be calculated by subtracting the elastic strain increment to the total strain increment (d p = d-d e ) with d v e = dp'/K and d d e = dq/(3G). Under axisymmetric condition we have  d = 2( a - r )/3 and q = ( a - r ). The Macaulay function in Eq. ( 8) implies that the negative

ACCEPTED MANUSCRIPT

plastic work due to shear induced dilation (d v p <0) does not influence grain breakage.

The plastic work was calculated through Eqs. ( 3), (4), and ( 8) from the stress-strain patterns.

Then, B u and B r * were plotted against the plastic work w p (Fig. 10). They can be expressed univocally, regardless of the loading pattern as:

* ; p p u r p p w w B B a w b w     (9)
with the fitting parameters a = 6000 and b = 1000. Eq. ( 9) implies that the two breakage indices increase asymptotically to the value of 1 when the plastic work increases. There is no further breakage once the ultimate grading is reached, which agrees with [START_REF] Coop | Particle breakage during shearing of a carbonate sand[END_REF]. Therefore, the plastic work appears to be a suitable variable to determine B u and B r * , and thus the evolution of the GSD, independently of the (monotonic or cyclic) stress path.

INFLUENCE OF GRAIN BREAKAGE ON THE EVOLUTION OF THE CSL

In the literature, three main methods for describing the effect of particle crushing on the CSL have been developed: (1) movement of CSL with the amount of plastic work (e.g., [START_REF] Daouadji | An elastoplastic model for granular materials taking into account grain breakage[END_REF][START_REF] Biarez | Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires[END_REF]; (2) direct description of CSL by one or several segments for different stress levels (e.g., [START_REF] Russell | A bounding surface plasticity model for sands exhibiting particle c rushing[END_REF]; (3) movement of CSL with the change of the grading state index I G [START_REF] Kikumoto | Particle crushing and deformation behaviour[END_REF]. The concept of the CSL relies on the assumption that, at the critical state, the material keeps its volume while undergoing deformation.

If the CSL moves with the grading changing because of grain crushing, this concept becomes dynamic and its evolution needs to be specified. To investigate this, we carried out drained triaxial tests under a confining stress of 100 kPa and different stress paths (Fig. 11a), and tests with constant ' 3 and constant p' (Fig. 11c) (see Section 4.1.3). The hypothesis of a movement of the CSL can be justified on the basis of the e-logp' patterns of the two tests shown by Fig. 11a with different amounts of crushing (Fig. 6c). If we assume that the CSL evolves with grain breakage, we need to determine its initial position. For our tested material, as well as for any easily crushable materials, it is difficult to measure the initial location of the CSL because significant breakage can occur even during triaxial tests at low confinement. We assumed that

ACCEPTED MANUSCRIPT

tests performed at 22 kPa could be used to indicate the location of the initial CSL because crushing, if not absent, was very limited and could not change the GSD significantly (Fig. 4c).

Furthermore, it is still unclear whether the slope of CSL is dependent or not on the GSD [START_REF] Daouadji | An elastoplastic model for granular materials taking into account grain breakage[END_REF][START_REF] Biarez | Influence de la granulométrie et de son évolution par ruptures de grains sur le comportement mécanique de matériaux granulaires[END_REF][START_REF] Bandini | The influence of particle breakage on the location of the critical state line of sands[END_REF]. Here for simplicity we hypothesized a constant slope of the CSL, so that the CSL can be expressed by: ln

c ref c at p ee p       ( 10
)
where e ref is the critical void ratio for a reference stress (for convenience equal to p at ), e c stands for the critical void ratio relative to the current p', and  c represents the slope of CSL in the e-lnp' plane (Fig. 11b).

Test results under different confining stresses (Fig. 11b) show the shift of the CSL with the confining stress increasing due to the increasing amount of grain breakage (Fig. 4c). Fig. 11c shows that the final value of the void ratio for the test under constant ' 3 = 400 kPa is lower than that for the test under constant p' = 400 kPa, due to the higher amount of grain breakage in the former test (Fig. 5c). Note that for tests under high confining stresses the critical state is not easily measurable (Fig. 4). [START_REF] Bandini | The influence of particle breakage on the location of the critical state line of sands[END_REF] proposed to perform a second (i.e.

unloading and shearing) loading phase to measure the critical state line under an appropriate, much lower confining stress, under which only negligible further crushing occurs. In this study, the final states at a strain level of 25% for all tests were chosen (arbitrarily) to represent the "dynamic" critical state. This specific treatment does not affect the trend of the shift in the CSL.

If we plot e ref for each test against the breakage measures B u and B r * (Fig. 12b), we can relate the evolutions of the CSL location and the GSD in the following way:

          (12)
In Eq. ( 12), the values of ( 1+c) and ( 1+d and Fig. 12b). It should be noted that the CSL location cannot be directly measured in cyclic tests.

If we assume that the results obtained from monotonic tests hold for more general loading conditions, than the CSL location for cyclic tests can be estimated from Fig. 12c and Eq. ( 12).

This seems a reasonable assumption, since the grain breakage measures appear to be a univocal function of w p for various loading conditions including monotonic and some cyclic loadings (Eq. 8). Therefore, it appears suitable to use the breakage measures B u and B r * to determine the CSL location, which is of practical importance for constitutive modelling. However, it should be pointed out that Eq. ( 8) should be still validated also through larger numbers of loading cycles under high confinement stresses

APPLICATION TO OTHER GRANULAR MATERIALS

The above relations between breakage measures, plastic work and CSL location obtained on the limestone grains were also used to interpret triaxial tests carried out on granite fragments, Cambria sand and Sydney sand. For the three granular materials, we plotted the relative uniformity and the modified relative breakage index versus the applied plastic work in Fig. 13a and Fig. 13b using Eq. ( 8). Furthermore, using Eq. ( 11), the reference critical void ratio was plotted with the relative uniformity and the modified relative breakage index in Fig. 13c and Fig. 13d, respectively. Comparisons between experimental results and calculations show that the proposed relationships based on limestone grains are also suitable for other crushable granular materials.

CONCLUSIONS

In this work, we employed the modified relative breakage index B r * and the relative uniformity B u as convenient measures for the engineering practice. Through triaxial tests on limestone fragments under monotonic loading, we demonstrated that the evolutions of these measures appear to be hyperbolic functions of the plastic work w p (which highlights their asymptotic nature). This relation was found to hold also under some cyclic loading conditions, suggesting that it should remain valid regardless of the loading nature. However, further tests remain necessary to confirm the relation also after larger numbers of loading cycles under high confinement stresses.

We find the location of the critical state line changing as breakage proceeds. In absence of further constraints, we assumed for simplicity that a parallel shift of the critical state line occurs as breakage proceeds. This is a strong assumption that must be verified carefully through further tests. If deemed acceptable, the location of the critical state line can be found by evaluating the plastic work at any stage of loading. We validated our approach by applying the relations we found for the limestone to other granular materials (granite fragments, Cambria sand and Sydney sand), obtaining satisfactory predictions.

FIGURES

  7) with B * r0 = 0.065, B * rmax = 0.215, a = 4; and B u0 = 0.0054, B umax = 0.037, b = 11 obtained by curve fitting.

  = 1.43, e ref = 0.5, and c = 0.006 and d = 0.035 by best fitting. The proposed formulation for the evolution of the CSL location with the grain breakage measures takes a different form with respect to the one proposed by Kikumoto et al. (2010) with the grading state index I G . Substituting Eq. (9) into Eq. (11), e ref can be also related directly to w p as follows:

  ) can be taken approximately equal to 1 because the values of c and d are much smaller than 1. Thus, only a•c and b•d control the evolution rate of the CSL location with the plastic work. It was found that the fitting curves based on B u and B r * in Fig. 12c are very close with a•c = 39 (based on Fig. 10a and Fig. 12a) and b•d = 35 (based on Fig. 10b

Fig. 1 .

 1 Fig. 1. Illustration of different measures of grain breakage.

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 23456789 Fig. 2. Photographs of a tested sample: (a) before testing, and (b) after testing under constant p' = 400 kPa.

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. (a) Relative uniformity and (b) modified relative breakage index versus plastic work.

Fig. 13 .

 13 Fig. 13. (a) Relative uniformity and (b) modified relative breakage index versus plastic work; (c) e ref against the relative uniformity, and (d) against the modified relative breakage index.

ACKNOWLEDGEMENTS

This research was supported by the National Basic Research Program of China: the funds for creative research groups of China (41521002) and the basic research funds (41790433). The work was also supported by the Sichuan funding for young researchers (2016JQ0021).


Maxime Lhuillier
email: maxime.lhuillier@free.frmaxime.lhuillier.free.fr

Improving thin structures in surface reconstruction from sparse point cloud

Keywords: Thin structures, Sparse features, 3D Delaunay triangulation, Visibility, Environment modeling

Methods were proposed to estimate a surface from a sparse cloud of points reconstructed from images. These methods are interesting in several contexts including large scale scenes, limited computational resources and initialization of dense stereo. However they are deficient in presence of thin structures such as posts, which are often present in both urban and natural scenes: these scene components can be partly or even completely removed. Here we reduce this problem by introducing a pre-processing, assuming that (1) some of the points form polygonal chains approximating curves and occluding contours of the scene and (2) the direction of the thin structures is roughly known (e.g. vertical). In the experiments, our pre-processing improves the results of two different surface reconstruction methods applied on videos taken by helmet-held 360 cameras.

Introduction

The automatic surface reconstruction from images is still an active research topic. One of the difficulties is caused by thin structures: they are partly or even completely removed by the methods that do not include specific processing. There are two kinds of thin structures: planar ones like road signs and rectilinear ones like posts. This paper only focuses on the rectilinear case, which is frequent in both urban and natural scenes. It also focuses on an input point cloud reconstructed from sparse features in images. This is useful not only for applications that do not need the high level of details of computationally expensive dense stereo, but also to initialize dense stereo in other cases. Since we have two difficulties (thin structures and sparse input) and would like a tractable problem, our method uses further information available in the images: we reconstruct curves as polygonal chains in 3D from image gradient. We also assume that the thin structures have a known direction with large tolerance. The paper restricts to the vertical direction which is the dominant direction of the reconstructed curves in our terrestrial image sequences (this is partly due to the aperture problem), but it is obvious to consider several directions by alternating.

Previous surface reconstruction methods

We summarize methods that use constraints similar to ours, whether they deal with thin structures or not, without or with a sparse point cloud.

Using lines or curves or image gradient

Methods use lines or curves and organize them as graphs, which are used as scaffolds for surface reconstruction. In [START_REF] Usumezbas | The surfacing of multiview 3d drawings via lofting and occlusion reasoning[END_REF], curves are reconstructed from image gradient and the surface is obtained by curve interpolation (lofting) and occlusion reasoning. In [START_REF] Wang | Image-based building regularization using structural linear features[END_REF], an initial surface mesh is regularized by back-projecting linear structures that are semi-automatically selected from the images. Other methods estimate dense depth maps including pixels at occluding contours, then merge them in a voxel grid. In [START_REF] Shan | Occluding contours for multi-view stereo[END_REF], quasi-dense depth maps of internet images are densified by encouraging depth discontinuities at image contours. In [START_REF] Yucer | Depth from gradients in dense light fields for object reconstruction[END_REF], the depths of video sequences are first computed in high-gradient regions including silhouettes, then they are propagated to low-gradient regions. Using a point cloud reconstructed from images of a scene with planar structures, [START_REF] Holzmann | Regularized 3d modeling from noisy building reconstructions[END_REF] segments into inside and outside a tetrahedralization of the points by using graph-cut and two regularizations: horizontal slicing and a smoothness term based on image lines.

A very recent work [START_REF] Li | Reconstructing thin structures of manifold surfaces by integrating spatial curves[END_REF] deals with thin rectilinear structures. First it reconstructs curves from image gradient, then integrates both dense scattered points and curve vertices in a 3D Delaunay triangulation, last segments into inside and outside by using graph-cut.

Reconstructing trees and roots

Methods estimate thin structures like trees [START_REF] Tan | Image-based tree modeling[END_REF][START_REF] Livny | Automatic reconstruction of tree skeletal structures from point clouds[END_REF] and plant roots [START_REF] Zheng | Detailed reconstruction of 3d plant root shape[END_REF]. In [START_REF] Livny | Automatic reconstruction of tree skeletal structures from point clouds[END_REF], several trees are automatically reconstructed from laser scan data assuming that the major branch structures are sampled. First every tree is modeled by a minimal spanning tree whose root is on high point density of the ground plane, vertices are input points and edge weights are Euclidean distances. Then the spanning trees are globally optimized using constraints on length, thickness and smoothness derived from biological properties. A previous method [START_REF] Tan | Image-based tree modeling[END_REF] reconstructs realistic-looking trees from pre-segmented images into branches/leaves/ background and a more difficult point cloud estimated by quasi-dense structurefrom-motion. The computation of minimum spanning trees is done for subsets of branches using weights combining two distances (3D Euclidean and image gradient-based). Little user interaction is needed to select branch subsets and to connect subsets. The method in [START_REF] Zheng | Detailed reconstruction of 3d plant root shape[END_REF] is shortly summarized in two steps: first robustly estimate a visual hull of the plant roots (They are mostly black with white background thanks to a transparent gel container.) in a voxel grid, then repair the connectivity of the inside voxels by computing a minimal spanning tree. In these methods and [START_REF] Li | Reconstructing thin structures of manifold surfaces by integrating spatial curves[END_REF], the connectedness constraint is essential to recover missing parts of the thin structures.

Other surface reconstruction methods

Other methods are not in the classifications above. In [START_REF] Chauve | Robust piecewise-planar 3d reconstruction and completion from large-scale unstructured point data[END_REF], the input point cloud is first approximated by planar segments, then these segments are prolonged and assembled into a well-behaved surface using visibility constraints. The prolongations are done using prior on urban scenes including the prevalence of the vertical structures. In [START_REF] Jancosek | Multi-view reconstruction preserving weakly-supported surfaces[END_REF], a graph-cut method segments a tetrahedralization of the input points into inside and outside by using visibility information. Compared to [START_REF] Vu | High accuracy and visibility-consistent dense multi-view stereo[END_REF], it favors the inside where the visibility gradient is high and improves the thin structure reconstructions, but it ignores image gradient and tends to produce more surfaces that do not correspond to real surfaces. In [START_REF] Oswald | Generalized connectivity constraints for spatio-temporal 3d reconstruction[END_REF], thin structures such as stick and rope are reconstructed in a voxel grid by enforcing the inside connectedness. First an approximate visual hull including the final surface is estimated, then a voxel labeling is computed by convex optimization subject to linear constraints on the labeling derivatives such that the inside and the visual hull have similar topology. Thin planar structures [START_REF] Ummenhofer | Point-based 3d reconstruction of thin object[END_REF] are not the paper topic.

Our contribution

There are several differences between our method and the previous ones that reconstruct thin (rectilinear) structures. First we reconstruct 3D curves from image gradient (including those that approximate occluding contours) and use them to detect the thin vertical structures. Their points also complete the sparse input point cloud reconstructed by structure-from-motion (SfM). Second we partition the space by using a 3D Delaunay triangulation T of the input points with tetrahedra labeled freespace or matter. Voxel grid is used in most previous works for thin structures but is inadequate for both large scale scenes and sparse input point cloud. Third we detect the thin structures in T as vertices of vertical curves that have small matter neighborhoods in the horizontal directions. These vertices are regrouped in large enough sets to avoid false positive detections. There is only one thin structure per set. Fourth we use the connectedness constraint to complete the thin structures in the vertical direction: force to matter vertical series of tetrahedra connecting vertices of a same set. Last we adapt two previous surface reconstruction methods such that their regularizations do not cancel our matter completion. Now we detail differences with the closest and recent work [START_REF] Li | Reconstructing thin structures of manifold surfaces by integrating spatial curves[END_REF]. It adds to T a lot of Steiner vertices near 3D curves, which are not in the original input and do not have visibility information (i.e. rays). This is not adequate for incremental surface reconstruction methods like [START_REF] Litvinov | Incremental solid modeling from sparse structure-frommotion data with improved visual artifact removal[END_REF][START_REF] Piazza | Real-time cpu-based large-scale threedimensional mesh reconstruction[END_REF] and is not obvious to adapt to nongraph-cut methods. Our thin structure processing does not use such points and does not have these drawbacks. The price to pay is (1) an approximate knowledge of the direction of thin structures and (2) involved detection-completion steps of the matter tetrahedra of the thin structures. Furthermore we start with terrestrial video sequences with sparse point clouds. We believe that this experimental context is more challenging (or "wild") than that of [START_REF] Li | Reconstructing thin structures of manifold surfaces by integrating spatial curves[END_REF], which starts from denser point clouds and still images whose view-points are selected by UAV.

Overview

The three first steps in Secs. 3.1, 3.2 and 3.3 estimate the curves, the vertical direction and initialize the 3D Delaunay triangulation from the output of a sparse SfM. The two next steps in Secs. 3.4 and 3.5 are the main contributions: they detect the thin vertical structures and improve their matter labeling. Last Secs. 3.6 and 3.7 explain how to update two surface reconstruction methods to take into account our matter improvement: a graph-cut method in [START_REF] Vu | High accuracy and visibility-consistent dense multi-view stereo[END_REF] and a manifold method in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF].

Step 1: estimate points and polygonal chains

We use standard methods: detect both Harris points and Canny curves, match them in consecutive images of the sequence using correlation and the epipolar constraint [START_REF] Lhuillier | Match propagation for image-based modeling and rendering[END_REF]. A curve is approximated by a polygonal chain of points in an image such that the distance between two consecutive points is roughly constant (4 pixels in our implementation). These points are matched with points of Canny curves in other images. Then the chain points are reconstructed independently and we obtain a polygonal chain in 3D. Last all 3D points are filtered by thresholding on angles formed by their rays. (Reject a point if its angles are smaller than 10 degrees.) Better methods exist to reconstruct lines and curves [START_REF] Kahl | Multiview reconstruction of space curves[END_REF][START_REF] Mai | 3d curves reconstruction from multiple images[END_REF][START_REF] Hofer | Line-based 3d reconstruction of wiry objects[END_REF][START_REF] Fabbri | 3d curve sketch: flexible curve-based stereo reconstruction and calibration[END_REF], but this is not the paper focus. The resulting polygonal chains approximate curves in 3D and have both redundancy and inaccuracy for several reasons: failures of detection and matching, matching of 2D points that do not correspond to a single 3D point due to occluding contours, inaccurate camera calibration.

Step 2: estimate the vertical direction

In the terrestrial imaging context, we observe that the main direction of the edges of the polygonal chains is vertical. (There are also a lot of horizontal curves in usual scenes, but most of them cannot be reconstructed due to the aperture problem.) Thus we compute the density of the edge directions in a histogram and assume that the direction with the largest density is vertical.

Step 3: initialize the 3D Delaunay triangulation

Step 3 is a by-product of the surface reconstruction methods based on 3D Delaunay triangulation and visibility constraints. First a 3D Delaunay triangulation T is build from the points generated in Sec. 3.1. Then its tetrahedra are labeled freespace or matter. In more details, each Delaunay vertex a has several rays and a ray is a line segment ac between a and a viewpoint c whose image includes a detection of a. The matter M is the set of tetrahedra that are not intersected by ray(s) and T \ M is the freespace. Note that we not only have isolated points but also edges of polygonal chains. Following [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF][START_REF] Sugiura | 3d surface reconstruction from point-and-line cloud[END_REF], a tetrahedron is freespace if it is intersected by a stereo triangle abc where ab is such an edge and c is a common viewpoint of a and b (more details in the supplementary material).

Step 4: detect vertices of the thin structures

For each thin structure, we would like to find (as large as possible) sets of vertices of T that are in the structure. Such a set is computed as a connected component of a graph G with vertex set V and edge set E. Since structure is vertical, all edges in E are almost vertical. The set E not only includes edges of the polygonal chains but also edges of the Delaunay triangulation. (A same structure can have many polygonal chains.) Since structure is matter, every vertex in V is in at least one matter tetrahedron in M . Since structure is thin along the vertical direction, such a vertex has a small connected component of matter tetrahedra in a horizontal neighborhood (in a slice between two horizontal planes). Sec. [START_REF] Hofer | Line-based 3d reconstruction of wiry objects[END_REF].2 details the definition of G.

Evidence of thin structure increases with connected component size. Thus we reject the smallest connected components of G to avoid false positive detections of thin structures. (If they have less than 6 vertices in our implementation.)

Step 5: improve the matter of the thin structures

Thin vertical structure should be approximated in T by connected set of matter tetrahedra with large size in the vertical direction and small size in the horizontal directions. However, this connectedness does not hold in practice due to several reasons: noise, lack of points and bad points that generate false positive freespace tetrahedra. We would like to improve this connectedness by using the selected connected components of G (Sec. 3.4).

Since all vertices of a connected component C of G are in a same structure, we force to matter (i.e. add to M ) series of tetrahedra in T connecting vertices of C. Since structure is vertical, the forced tetrahedra cover line segments that are almost vertical. Since structure is thin, we only force to matter tetrahedra with small size in the horizontal directions. Thus there are two steps: first estimate the width (measuring the size in the horizontal directions) of the structure including C, then complete M by considering tetrahedra between pairs of vertices in C if the width of these tetrahedra are not too large compared to that of the structure. This process is detailed in Sec. 4.3.

Last we set a Boolean b ∆ for every tetrahedron ∆: b ∆ is true if and only if ∆ is forced to matter by Step 5. (Whether we have ∆ ∈ M or ∆ / ∈ M in Step 3.)

Step 6a: graph-cut surface reconstruction

The graph-cut method [START_REF] Vu | High accuracy and visibility-consistent dense multi-view stereo[END_REF] labels the tetrahedra of T with inside or outside.

(The surface is the triangle set between outside and inside tetrahedra.) However, it removes thin structures by labeling outside several tetrahedra that are in M . There are two reasons: it has a regularization and it ignores our matter improvement. Thus we modify it such that the set of the inside tetrahedra includes the set of the tetrahedra with true Boolean b ∆ (Sec. 3.5). We remind that it computes a minimum s -t cut of a graph. This graph is the adjacency graph of T augmented by two vertices {s, t} and edges {vs, vt} for every vertex v of the adjacency graph. Cut is partition of the vertex set in two subsets: a subset including s and a subset including t. The details on the definition of the augmented graph and its edge weights are in [START_REF] Vu | High accuracy and visibility-consistent dense multi-view stereo[END_REF]. We only need to know that if vt has a large weight, v and t are in the same subset, which implies that the tetrahedron associated to v is labeled inside.

Thus we modify the augmented graph (before computing the minimum cut) by adding a large value to every edge vt such that the tetrahedron ∆ associated to v has a true Boolean b ∆ .

Step 6b: manifold surface reconstruction

The manifold method [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] also labels the tetrahedra of T with inside or outside. It uses our matter improvement (i.e. the set M computed in Sec. 3.5), but it removes thin structures due to a regularization, which labels outside several tetrahedra that are in M .

We summarize it to correct the regularization. Let O be the set of outside tetrahedra and ∂O be its boundary, i.e. the set of triangles between outside and inside tetrahedra. First O grows in T \ M to meet visibility constraints, then O evolves in T by the regularization. Both steps maintain the manifoldness of ∂O. The regularization removes peaks, i.e. rings of triangles in ∂O around a vertex a which form small solid angles. This is done in two cases: remove from O all tetrahedra including a, or add to O all tetrahedra including a. However parts of thin vertical structure can be removed in the second case: some of their matter tetrahedra can become outside the surface.

We solve this problem as follows: the regularization is prohibited in the second case if at least one of the tetrahedra added to O has a true Boolean b ∆ .

Detection and improvement of the structures

Here we detail the contributions of the paper summarized in Secs. 3.4 and 3.5.

Definitions

We remind that M is the set of the matter tetrahedra in T , which should include the thin vertical structures. Let v ∈ R 3 be the vertical direction. (We have ||v|| = 1.)

The small matter slice M (∆ 0 ) started from a tetrahedron ∆ 0 ∈ M is computed as follows. Let P 1 and P 2 be the two horizontal planes that enclose ∆ 0 . (The vertex of ∆ 0 with the highest or lowest altitude is in P 1 or P 2 , both altitude and horizontal directions are defined using v.) We initialize M (∆ 0 ) = {∆ 0 } and complete it by a graph traversal of the adjacency graph of M by only considering tetrahedra that intersect the area between P 1 and P 2 . In more details, the following process is repeated: take a tetrahedron ∆ ∈ M (∆ 0 ) and add to M (∆ 0 ) every neighbor tetrahedron of ∆ in M (∆ has no more than four neighbors) that intersects the area between P 1 and P 2 . We only need matter slices whose number of tetrahedra is less than a small threshold M max , i.e. which can be included in a thin vertical structure. (We use M max = 20.) Thus we stop the graph traversal if M (∆ 0 ) becomes larger than M max , then we reset M (∆ 0 ) = ∅ (no slice found).

Let a and b be in

R 3 . The horizontal distance w(a, b) is w(a, b) = ||a -b|| 2 -(v (a -b)) 2 . (1)
Let v(A) be the vertex set of a tetrahedron set A ⊆ T . The maximum width of A is measured from a base point b ∈ R 3 :

w(A, b) = max{w(a, b), a ∈ v(A)}. (2)
We also define the path p(∆, a) between a tetrahedron ∆ ∈ T and a vertex a ∈ v(T ): this is the set of tetrahedra in T whose interiors are intersected by the line segment joining a and the barycenter of ∆.

Vertex detection

We only need to detail the definition of the graph G = (V, E) introduced in Sec. 3.4. Let α be an angle threshold. Let E be the set of the edges ab in the polygonal chains that are almost vertical, i.e. if

|v b -a ||b -a|| | > cos α. (3)
We use a large tolerance α = 20 • is our implementation. Let V be the vertex set of E . We have a ∈ V if and only if a ∈ V and there is a tetrahedron ∆ 0 ∈ M including a such that M (∆ 0 ) = ∅. We have ab ∈ E if and only if a ∈ V and b ∈ V and Eq. 3 is meet and the following condition is meet: ab ∈ E or ab is an edge of the Delaunay triangulation or there is c ∈ V such that ac ∈ E and cb ∈ E . The last condition using c ∈ V is useful to connect two parts of a thin structure separated by a vertex that is not in a matter tetrahedron.

Matter improvement

First we estimate a median width w C for every connected component C of G by using the horizontal distance between a vertex b ∈ C and other vertices of the small matter slices M (∆ 0 ) such that b ∈ ∆ 0 . We collect the small matter slices attached to b ∈ C in

M b = ∪ b∈∆0∈M M (∆ 0 ) (4)
and compute

w C = median{w(a, b), b ∈ C, a ∈ v(M b ), a = b}. ( 5
)
Second the vertices a i ∈ C are ordered along the vertical direction v: they meet v a i ≤ v a i+1 for all i. We obtain one composite polygonal chain a 1 a 2 , a 2 a 3 , • • • for C, which should be included in matter tetrahedra of the thin structure (it is not expected to be exactly on the surface of the thin structure).

Third we collect in a list L for several pairs (i, j) a path between a i and a j that will be forced to matter. Let b = (a i + a j )/2. We consider several paths between a i and a j and take one that minimizes the maximum width with respect to b. These paths are p(∆, a j ) where a i ∈ ∆ ∈ M and p(∆, a i ) where a j ∈ ∆ ∈ M . For example, p(∆, a j ) minimizes the maximum width for a given (i, j). If w(p(∆, a j ), b) < cw C where c is a threshold, we decide that p(∆, a j ) is thin enough to be forced to matter and we add (∆, a j ) to the list L. This is done for all C. In the experiments, we use c = 2 and only take j ∈ {i + 1, i + 2}.

Last the matter is improved: we add p(∆, a j ) to M for every (∆, a j ) ∈ L. Fig. 1 summarizes the construction and the use of a composite polygonal chain. Second the vertices of these chains form a connected component of G and we draw the resulting composite polygonal chain (eight bold edges). Third we draw all line segments linking the barycenter of matter triangles and vertices such that j = i + 1. Fourth all triangles intersected by these segments are forced to matter.

Experiment

Steps 1-5 (our pre-processing) are described in Sec. 5.1. Secs. 5.2 and 5.3 focus on step 6a (graph-cut methods) and step 6b (manifold methods), respectively. Last Sec. 5.4 shows a lot of results for all surface reconstructions methods. For fair comparisons, the experimented methods start from the same 3D Delaunay tetrahedralization T and end by the same Laplacian smoothing of the surface.

First we take two 2496 × 2496 videos at 30Hz by walking during 473s in a town using a helmet-held 360 camera (Garmin Virb 360). Then we apply generic SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] followed by detection and closure of loops inspired by [START_REF] Yousif | An overview to visual odometry and visual slam: application to mobile robotics[END_REF][START_REF] Younes | Keyframe-based monocular slam: design, survey, and future directions[END_REF]. Last we experiment on the 1334 keyframes selected by SfM. Fig. 2 shows our setup, a keyframe, a top view of the SfM result, a view of the 3D model by one of our method, and the set of vertical edges that we use.

Our pre-processing

Step 1 generates 106k polygonal chains having 249k vertices. Step 3 builds the 3D Delaunay triangulation T . It has 1.2M vertices reconstructed from image Step 5 is the matter improvement step. It forces 10k tetrahedra to be matter. (6.5k of them were in T \ M before.) Fig. 3 shows our pre-processing: it has a noticeable improvement on a low-textured thin structure with missing matter tetrahedra and a negligible effect on a well-textured thick structure with a lot of matter tetrahedra. These results are those expected.

The detection and improvement steps take less than 0.6s on a standard laptop. The feature reconstruction (both Canny and Harris points) is the most time consuming step: 40 minutes.

Graph-cut surface reconstruction methods

Here we compare methods "base graph-cut" and "our graph-cut". The former is in [START_REF] Vu | High accuracy and visibility-consistent dense multi-view stereo[END_REF] and the latter is an update of the former (Sec. 3.6) preceded by our pre-processing (Sec. 3). We remind that base graph-cut estimates a surface minimizing an energy E vis + λ qual E qual with a visibility term E vis (using α vis = 1) and a quality term E qual . The greater the value of λ qual , the better the surface quality in a regularization sense. Fig. 4 shows a light post reconstructed by our graph-cut using λ qual = 2 and base graph-cut using λ qual = 0.5 and λ qual = 2. The surface of the complete scene has 1.888M, 2.183M and 1.885M Fig. 3. Our pre-processing for a post (top) and a trunk (bottom). From left to right: one input image, reconstructed vertical edges (using α = 20 • ), matter tetrahedra before our pre-processing, composite polygonal chains, matter tetrahedra after our pre-processing, final surface. triangles, respectively. First we see that choosing a small λ qual is a first solution to improve the post. However this solution is not satisfactory for two reasons: the post is incomplete (about 30% is lacking) and the noise increases in the background house. See the gray level variations in the left part. The smaller the value of λ qual , the greater the number of surface vertices. (They include inaccurate points.) Second, we see that our graph cut provides the most complete post at the price of very small increases of triangle number and computation time. In the next experiments using base or our graph-cut, λ qual = 2.

Manifold surface reconstruction methods

Now we compare methods "base manifold" and "our manifold". The former is in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] and the latter is an update of the former (Sec. 3.7) preceded by our preprocessing (Sec. 3). Fig. 5 shows the same light post as Fig. 4 reconstructed by base manifold and our manifold. The surface of the complete scene has 1.734M and 1.739M triangles, respectively. Almost all parts of the post are removed by base manifold. In contrast to this, the post is almost complete by our manifold (at the price of very small increases of triangle number and computation time). There are two spurious handles that connect the post and the house. Thanks to other examples, we observe that the probability of such an artifact increases when the distance between a post and another scene component (house or wall or another post) decreases. Fig. 5 also shows base manifold with the regularization update in Sec. 3.7, and our manifold without the regularization update (i.e. base manifold with our pre-processing). In both cases, the post reconstruction is less complete than that of our manifold. This shows that both our pre-processing and regularization update are useful.

Other experiments

In other examples we compare base graph-cut and our graph-cut, base manifold and our manifold. Fig. 6 shows results for electric posts similar to that in the experiments above with different backgrounds. Our methods provide the most complete post reconstructions. Examining the background house in the top of the figure and assuming its planarity, we compare noises: the graph-cut noise is greater than the manifold noise although the total numbers of surface triangles are similar. Examining the post in the bottom of the figure, we note a shrinking artifact due to the Laplacian smoothing and a fattening artifact due to low texture in the background. (There is a of points to carve the freespace near the post bottom, thus the foreground post and background are not separated.) This holds for all methods. Figs. 7 and8 show thin structures in another dataset: a 2.5km long sequence using four helmet-held Gopro cameras biking in a city [START_REF] Nguyen | Self-calibration of omnidirectional multi-camera including synchronization and rolling shutter[END_REF]. Although our results on thin structures are not perfect, they are clearly better than those of the two base methods which omit important parts. Furthermore we remind that they are obtained in a difficult context: sparse point clouds from videos sequences taken in the wild.

We also examine the differences between base methods and our methods in overall reconstructed scenes, and see few blunders of our preprocessing due to excesses of the matter completion. There are several reasons. A first one is the conservative selection of vertices in Sec. 3.4: the small connected components of G are discarded for robustness (for both detection and w C computation). Second our matter completion forces only 10k tetrahedra to matter. This number is small compared to the 3.1M tetrahedra that are matter by default, i.e. that are not intersected by rays or stereo triangles in Sec. 3.3. Third the tetrahedra forced to matter are selected in Sec. 4.3 such that they have smallest widths as possible in the horizontal directions. This limits the size of a potential blunder.

Conclusion

This paper presents a new pre-processing to improve the thin rectilinear structures in a surface reconstructed from an image sequence using a sparse cloud of features and their visibility. Starting from the traditional 3D Delaunay triangulation, we first detect vertical series of vertices in a same thin structure by reconstructing curves. Then we complete the matter label in tetrahedra between such vertices by using the connectedness constraint. Last we easily adapt two previous surface reconstruction methods such that their regularizations do not remove the thin structures that we detect and complete.

Our method has limitations and can be improved by several ways. The thin structures cannot be detected if they don't have gradient edges in the images. Since the detection step relies on polygonal chains approximating curves and occluding contours in the scene, this step can be improved thanks to previous work on curve reconstruction from images. Furthermore, the completion step should not only enforce connectedness but also geometric regularity like constant or smooth width along the thin structures. Our contributions can be applied to the incremental surface reconstruction and several directions of thin structures (although we only experiment batch surface reconstruction methods with the vertical direction). Last our pre-processing should also be generalized to other thin structures such as road signs, which are not rectilinear but planar.

Fig. 1 .

 1 Fig.1. Construction and use of a composite polygonal chain in a simplified case: in 2D using j = i + 1 without the width test. (The vertical direction is horizontal in this figure.) First there are three reconstructed polygonal chains (bold edges) in the Delaunay triangulation. The gray triangles are matter and the remainder is freespace. Second the vertices of these chains form a connected component of G and we draw the resulting composite polygonal chain (eight bold edges). Third we draw all line segments linking the barycenter of matter triangles and vertices such that j = i + 1. Fourth all triangles intersected by these segments are forced to matter.

Fig. 2 .

 2 Fig. 2. Left: the helmet-held Garmin Virb 360 camera and one keyframe. Middle: a top view of the SfM point cloud and the camera trajectory. Right: A view of the final surface reconstructed by step 6b and the vertical edges (the set E ) that we use.

features and 7 .

 7 4M tetrahedra, 3.1M of them are matter. Step 4 is the detection step of thin structures. It only retains 25k edges of the polygonal chains that are almost vertical (the set E ). Then it finds 13k vertices in the vertex set of E that have a small matter slice in their neighborhood (the set V ). Using edges in both Delaunay triangulation and E , the graph G has 7k connected components. Last Step 4 only retains 230 of them by ignoring the smallest ones (very end of Sec. 3.4).

Fig. 4 .

 4 Fig. 4. Comparing three graph-cut methods. From left to right: base graph-cut using λ qual = 0.5, base graph-cut using λ qual = 2, our graph-cut using λ qual = 2, a local view of an input image.

Fig. 5 .

 5 Fig. 5. Comparing four manifold methods. From left to right: base manifold, base manifold with updated regularization, our manifold with naive regularization, our manifold, a local view of an input image.

Fig. 6 .

 6 Fig. 6. Comparisons between several methods. From left to right: base graph-cut, our graph-cut, local view of input image, base manifold, our manifold.

Fig. 7 .

 7 Fig. 7. Comparisons between several methods. From left to right: base graph-cut, our graph-cut, our manifold (textured), base manifold, our manifold.

Fig. 8 .

 8 Fig. 8. Comparisons between several methods. Top: base graph-cut (left) and our graph-cut (right). Middle: base manifold (left) and our manifold (right). Bottom: local view of an input image.