Fabien Treyssède
email: fabien.treyssede@lcpc.fr

Free linear vibrations of cables under thermal stress

This paper aims at investigating the effect of thermal stress on free linear vibrations of cables. An analytical model extending Irvine's theory to thermoelasticity is proposed. In addition to the sag-extensibility parameter, this model yields a second dimensionless independent parameter related to the temperature change. Some general two-dimensional contour plots are presented giving the relative change in natural frequencies due to temperature with respect to the only two independent parameters of the problem. Depending on their values, this change can be quite significant. For in-plane symmetric modes, it is shown that a positive (resp. negative) temperature change does not necessarily decrease (resp. increase) natural frequencies and that thermoelastic effects can lead to cross-over due to the modification of the initial sag.

Introduction

Cables are widely used in modern structures (bridges, post-tensioned concrete, suspended roofs, skilifts, elevators). The study of their vibrations has a long history, since the eighteenth century (see Refs. [START_REF] Irvine | Cable Structures[END_REF][START_REF] Triantafyllou | Linear dynamics of cables and chains[END_REF]). Irvine's work [START_REF] Irvine | Cable Structures[END_REF] is probably one of the most important contributions, its main advantage being that only one dimensionless parameter (namely, the sag-extensibility parameter) is needed to determine cable natural frequencies. Since then, many studies have aimed at taking into account more complicated mechanics by including bending stiffness [START_REF] Mehrabi | Unified finite difference formulation for free vibration of cables[END_REF][START_REF] Ni | Dynamic analysis of large-diameter sagged cables taking into account flexural rigidity[END_REF] or non-linear dynamics effects [START_REF] Rega | Nonlinear vibrations of suspended cables. part i: Modeling and analysis; part ii: Deterministic phenomena[END_REF][START_REF] Srinil | Three-dimensional non-linear coupling and dynamic tension in large-diameter free vibrations of arbitrary sagged cables[END_REF][START_REF] Zhu | Elastodynamic analysis of low tension cables using a new curved beam element[END_REF][START_REF] Srinil | Nonlinear longitudinal/transversal modal interactions in highly extensible suspended cables[END_REF][START_REF] Rega | Experimental and numerical studies of inclined cables: Free and parametrically-forced vibrations[END_REF] for instance. However to the author's knowledge, the investigation of thermal effects on cable dynamics Fig. 1. Cable profile for the initial static equilibrium state (thick line), thermally stressed state (thin line) and dynamic state (dashed line). The figure gives an example for a positive thermal change, corresponding to a decrease in the applied tension at the end and an increase of sag. The dynamic profile here corresponds to a second symmetric in-plane mode.

has surprisingly not received a similar attention in the literature. The goal of this paper it to extend Irvine's model to thermally stressed cables.

One of the motivations is the potential need of adequate models for vibration based methods in structural health monitoring (SHM). These methods are potentially attractive for tension estimation [START_REF] Ren | Empirical formulas to estimate cable tension by cable fundamental frequency[END_REF][START_REF] Kim | Estimation of cable tension force using the frequency-based system identification method[END_REF] or damage detection [START_REF] Bouaanani | Numerical investigation of the modal sensitivity of suspended cables with localized damage[END_REF][START_REF] Lepidi | Static and dynamic response of elastic suspended cables with damage[END_REF] in cables, but are likely to suffer a lack of robustness because of environmental temperature change (affecting the initial cable tension and in turn its modal parameters). It is well-known that differentiating changes due to the environment from changes due to damage is still a challenging task [START_REF] Rohrmann | Structural causes of temperature affected modal data of civil structures obtained by long time monitoring[END_REF][START_REF] Peeters | Vibration-based damage detection in civil engineering: excitation sources and temperature effects[END_REF][START_REF] Balmès | Merging sensor data from multiple temperature scenarios for vibration monitoring of civil structures[END_REF][START_REF] Deraemaecker | Vibration-based structural health monitoring using output-only measurements under changing environment[END_REF]. In addition to SHM, one could also note that temperature changes may also affect the robustness of vibration control strategies [START_REF] Na | Dynamics and active bending vibration control of turbomachinery rotating blades featuring temperature-dependent material properties[END_REF][START_REF] Parisse | Dynamic response of a system driven by thermal actuation[END_REF].

This study focuses on moderate thermal stress typically due to climatic variations. In Sec. 2, the thermoelastic equilibrium equations are derived both for cable statics and linear dynamics. A second independent parameter, related to temperature change, is naturally introduced. Section 3 gives general contour plots of the relative change in natural frequencies with respect to the only two independent parameters of the problem. Section 4 finally concludes this paper.

Theory

Figure 1 depicts the cable profile for its three equilibrium states: the initial static state (without thermal change), the thermally stressed static state and the superimposed dynamic state. The present study is restricted to small linear superimposed vibrations.

Statics without thermal stress

We consider a suspended cable anchored on supports at the same level, l being its span. Due to its own weight, the cable has a sag, denoted d. H is the horizontal component of the external applied tension. τ is the tension inside the cable (the notation T is usually chosen in the literature but, in this paper, it will be kept for temperature). mg is the self-weight of the cable per unit length. The static equilibrium of the cable is given by [START_REF] Irvine | Cable Structures[END_REF]:

d ds τ dx ds = 0, d ds τ dz ds = -mg (1)
yielding:

τ dx ds = cst = H, H d 2 z dx 2 = -mg 1 + dz dx 2 (2)
where we have used the equality ds/dx = 1 + (dz/dx) 2 , obtained from the geometric constraint (dx/ds) 2 + (dz/ds) 2 = 1.

For small sag (typically, d/l < 1/8), one can assume (dz/dx) 2 << 1 and the static equilibrium equation simplifies into [START_REF] Irvine | Cable Structures[END_REF]:

H d 2 z dx 2 = -mg (3)
The boundary conditions are:

z(0) = z(l) = 0 (4)
The solution for the cable profile is then:

z(x) = mgl 2 2H x l 1 - x l (5)
and the sag is:

d = z(l/2) = mgl 2 /8H (6)
One must have H > mgl (in order for d/l to be less than 1/8). The cable length is L = l 0 1 + (dz/dx) 2 dx, which can be approximated by a first order Taylor expansion as:

L = l 1 + m 2 g 2 l 2 24H 2 (7)

Thermally stressed state

Under thermal stress, an increment h (positive or negative depending on the sign of temperature change) is added to the initial applied tension H. The static equilibrium now reads:

d ds (τ + τ ′ ) dx ds + du ds = 0, d ds (τ + τ ′ ) dz ds + dw ds = -mg (8)
where u and w are the additional longitudinal and vertical cable displacement owing to temperature effects. τ ′ is the additional tension in the cable. Neglecting u and taking into account Eq. ( 2), the first equation in [START_REF] Srinil | Nonlinear longitudinal/transversal modal interactions in highly extensible suspended cables[END_REF] becomes:

τ ′ dx ds = h (9)
and using the same simplifications as before, Eqs. ( 8) yield the following equilibrium equation:

(H + h) d 2 dx 2 (z + w) = -mg (10)
with boundary conditions:

w(0) = w(l) = 0 (11)
Taking into account the equilibrium state (3) into [START_REF] Ren | Empirical formulas to estimate cable tension by cable fundamental frequency[END_REF], integrating and using boundary conditions [START_REF] Kim | Estimation of cable tension force using the frequency-based system identification method[END_REF] yields:

w(x) = - h H + h z(x) (12)
At this stage, h is still unknown and must be evaluated. Denoting ∆T as the temperature change (assumed uniform throughout this paper), the generalized Hooke's law of linear thermoelasticity stipulates that [START_REF] Boley | Theory of Thermal Stresses[END_REF]:

τ + τ ′ = EA(ǫ + ǫ ′ -α∆T ) (13)
where τ = EAǫ from the initial equilibrium state (ǫ being the initial strain). E, A and α respectively denote the Young modulus, the cross-section area and the thermal expansion coefficient. ǫ ′ is the additional thermal strain:

ǫ ′ = ds ′2 -ds 2 2ds 2 ( 14
)
where ds is the original length of the element before temperature change and ds ′ is its new length:

ds 2 = dx 2 + dz 2 , ds ′2 = (dx + du) 2 + (dz + dw) 2 (15)
Using Eqs. ( 14), ( 15) and ( 9), and ignoring the high order term (du/dx) 2 /2, Eq. ( 13) yields after multiplication by (ds/dx)

Integrating from 0 to l and by parts the second term of the right hand side gives:

hL e EA = u(l) -u(0) + mg H l 0 wdx + 1 2 l 0 dw dx 2 dx -α∆T L t ( 17
)
where L e = l 0 (ds/dx) 3 dx and L t = l 0 (ds/dx) 2 dx. From ds/dx = 1 + (dz/dx) 2 , we have (Taylor expansion):

L e = l 1 + m 2 g 2 l 2 8H 2 , L t = l 1 + m 2 g 2 l 2 12H 2 (18)
It might be useful to consider flexible horizontal supports at the ends, simply modelled as linear springs of stiffness k 0 and k l :

h = k 0 u(0), h = -k l u(l) (19)
Using Eqs. ( 5), ( 12) and ( 19) into (17) finally gives the following equation for the additional tension:

h * 3 + (2 + θ + λ 2 /24)h * 2 + (1 + 2θ + λ 2 /12)h * + θ = 0 ( 20
)
where h * = h/H, λ 2 = (mgl/H) 2 l/(HL e /E eq A eq ), θ = α∆T L t /(HL e /E eq A eq ), and E eq A eq = EA/{1 + (1/k 0 + 1/k l )EA/L e }. λ 2 is the Irvine sag-extensibility parameter. θ is a dimensionless parameter related to the temperature change.

Note that θ = 0 yields h * = 0. It is emphasized that Eqs. ( 12) and ( 20) are equivalent to the result given without proof in Exercise 2.5 of Irvine's book [START_REF] Irvine | Cable Structures[END_REF] (the present form being more compact).

Free linear vibrations under thermal stress

The dynamic equilibrium of the cable is given by the following equations:

d ds (τ + τ ′ + τ ) dx ds + du ds + dũ ds = m ∂ 2 ũ ∂t 2 , (21)
d ds (τ + τ ′ + τ ) dṽ ds = m ∂ 2 ṽ ∂t 2 , (22)
d ds (τ + τ ′ + τ ) dz ds + dw ds + d w ds = m ∂ 2 w ∂t 2 -mg (23)
where ũ, ṽ, w are the longitudinal, out-of plane and vertical dynamic components. τ is the additional dynamic tension inside the cable. We apply the same simplifications as before. The longitudinal component displacement is neglected as in the previous subsection. The equations are linearised by assuming that superimposed dynamic perturbations are small. After substituting the equations of static equilibrium, the equations of motion with harmonic time dependence are finally reduced to:

(H + h) d 2 ṽ dx 2 + mω 2 ṽ = 0, ( 24
) (H + h) d 2 w dx 2 + mω 2 w = mg H + h h (25)
h is the external dynamic tension. The boundary conditions for the dynamics are given by: ũ(0) = h/k 0 , ũ(l) = -h/k l , ṽ(0) = ṽ(l) = 0, w(0) = w(l) = 0 (26)

Following the same procedure as in the previous subsection, the application of Hooke's law yields after linearisation the following equation for h:

h EA ds dx 3 = dũ dx + dz dx + dw dx d w dx (27)
and after integration: hL e

E eq A eq = mg H + h l 0 wdx (28)
Out-of-plane motion The eigensolutions of Eq. ( 24) together with the boundary conditions (26) are straightforward:

Ω ′ n = nπ, ṽn (x) = C n sin nπ l x (29)
where Ω = ωl(H/m) -1/2 and Ω ′ = Ω(1 + h * ) -1/2 respectively denote the initial and modified dimensionless angular frequencies. n is a positive integer.

In-plane motion For the solutions of Eq. ( 25) with boundary conditions (26), it is convenient to distinguish symmetric modes from antisymmetric ones. From Eq. ( 28), antisymmetric modes do not induce additional dynamic tension ( h = 0), so that their natural frequencies and modeshapes are:

Ω ′ n = 2nπ, wn (x) = A n sin 2nπ l x (30)
For the calculation of symmetric modes, one needs to find a peculiar solution of Eq. ( 25), which can simply be the constant hg/ω 2 (H + h). After application of boundary conditions, it can be shown that the general solution for symmetric modes is:

wn (x) = mgl 2 h (H + h) 2 Ω ′2 n 1 -tan Ω ′ n 2 sin Ω ′ n x l -cos Ω ′ n x l (31)
Using Eq (31) into (28), we obtained the following transcendental equation for the natural frequencies Ω ′ n of symmetric modes:

tan Ω ′ 2 = Ω ′ 2 - 4 λ ′2 Ω ′ 2 3 ( 32
)
with λ ′2 = λ 2 /(1 + h * ) 3 . One remarks that this equation is the same as the well-known transcendental equation of Irvine, but with differently defined dimensionless parameters Ω ′ and λ ′2 (for h * = 0, Irvine's equation is exactly recovered).

Results

The quantity of interest is the relative change in natural frequency due to temperature for a given cable (fixed λ 2 ). Frequency values without thermal stress are well-known [START_REF] Irvine | Cable Structures[END_REF]. For clarity, their evolution with respect to the sagextensibility parameter is recalled in Fig. 2 (given Eq. (32), this figure also corresponds to the evolution of Ω ′ versus λ ′2 ). Provided that one is interested in relatively low temperature change (typically due to climatic variations), the influence of temperature on material properties is neglected in the following results (without loss of generality for the proposed theory).

The only independent parameters of the problem are λ 2 and θ, so that quite general results are obtained through two-dimensional contour plots of ∆ω/ω = f (λ 2 , θ), where ∆ω/ω = ω(λ 2 , θ)/ω(λ 2 , 0) -1 is the relative change in natural frequency under the influence of temperature. It is intuitively clear that for a fixed temperature change ∆T , the ratio H/EA must be low enough in order to have a non-negligible value of θ -and hence a non-negligible effect of temperature on frequencies. In other words, cables that are strongly pretensioned will be less sensitive to temperature change. In this paper, Eqs. ( 20) and ( 32) are solved with a Newton-Raphson algorithm for λ 2 varying from 1 to 200 and θ varying from -1 to +1.

For out-of-plane and antisymmetric in-plane modes, the relative change in natural frequency is analytically obtained in terms of h * as:

∆ω/ω = √ 1 + h * -1 (33)
Figure 3 gives the contour plots for √ 1 + h * -1 (x100) in terms of θ and log λ 2 . For every value of λ 2 , the eigenfrequencies of such modes decrease as the temperature change increases. One can observe a slightly higher sensitivity to cooling than heating. However for fixed values of θ, this sensitivity to temperature becomes lower and lower for cables having a larger and larger value of λ 2 . Note that comparing temperature sensitivity for different values of λ 2 at fixed θ implies that the ratio between the working stress H/A and Young's modulus E should remain almost constant. For a given cable material, comparisons for fixed θ are hence indeed made for an almost constant safety factor, which is meaningful.

For symmetric in-plane modes, Eq. ( 32) must be solved. The effect of temperature on frequencies manifests itself through a modification of the axial tension (as for out-of-plane and antisymmetric modes) as well as through a sag modification, yielding a modified parameter λ ′2 as shown by Eq. (32). Figure 4 plots the relative change of the first symmetric natural frequency. The thermal behaviour of this mode is strongly dependent on λ 2 and quite different from Fig. 3. Roughly, the frequency relative change is rather small above log λ 2 = 4.5 (λ 2 ≃ 90). Between 2.5 and 4.5 (λ 2 between 12 and 90), it is more pronounced and the frequency is counterintuitively increasing with temperature. As shown later, such increase can lead to cross-over.

Below log λ 2 = 2.5 (λ 2 ≃ 12), the behaviour changes again. For θ < 0 (cooling), the first frequency is increasing as the temperature is decreasing. However for θ > 0 (heating), the modal behaviour is more complex: the frequency might increase or decrease depending on the value of λ 2 as well as of θ. Note that in certain region, the frequency is not monotically varying with respect to temperature change, as it is the case for log λ 2 = 2 (λ 2 ≃ 7) for instance: the frequency tends to increase for any negative or positive temperature change.

The fact that the frequency can increase with temperature is explained as follows. A temperature change modifies the axial tension as well as the cable sag. When the temperature increases, the axial tension always decreases while the sag-to-span always increases (as shown by Eq. ( 12) and (33) together with the results of Fig. 3). A decrease of tension tends to decrease all natural frequencies, while the increase of sag tends to increase the frequencies of symmetric in-plane eigenmodes: tension and sag actions have hence a counteracting effect. This can be formally proven from Fig. 2: one has Ω ′ n = f (λ ′2 ), where f is a monotically increasing function, which can be rewritten as:

Ω n = f λ 2 (1 + h * ) 3 √ 1 + h * (34)
Hence when h * is decreasing (increasing temperature), the term √ 1 + h * decreases while the function f (•) increases: this means that the decrease of tension can be compensated by the increase of the sag-extensibility parameter. Note that for taut strings (λ 2 = 0), frequencies are given by Ω n = nπ √ 1 + h * and hence always decrease for positive thermal changes (their relative change is the same as in Fig. 3).

Figs. 5 and 6 also give the contour plots for the relative change in frequency of the second and third in-plane symmetric modes. As for the first mode, their thermal behaviour is still significantly different from anti-symmetric modes. However, differences become weaker as the mode order increases and the behaviour tends to become identical to that of antisymmetric and out-of-plane modes (compare Fig. 6 with 3).

As a last example, one considers a cable subjected to a temperature change ranging from -40 to +40 K. The cable has the following characteristics: l=200 m, A=7.069e-2 m 2 , E=2.0e+11 Pa, ρ=7800 kg.m -3 , α=1.2e-5 K -1 , g=9.81 m.s -2 . The initial applied tension is H=0.938e+7 N so that λ 2 =20.01 and θ ranges from -0.723 to +0.723. Horizontal supports are supposed to have infinite stiffness. Figure 7 exhibits the evolution of the first in-plane eigenfrequencies versus temperature change. Excepted the first one, all natural frequencies are decreased with heating. The frequency relative change for the first and second symmetric modes respectively varies from -7.0 to +11.6% and +15.2 to -9.5%, while it varies from +16.1 to -11.3% for antisymmetric modes (these values can be also directly found from Figs. 345). Also sketched in Fig. 7, the evolution of the first modeshape appears to be significant. Indeed, one cross-over between the first symmetric and anti-symmetric modes occurs near +40 K. This due to the fact that for ∆T varying from -40 to +40 K, the modified Irvine parameter λ ′2 ranges from 8.17 to 40.97, which is greater than the well-known cross-over value of 4π 2 . In the present example, it should be noticed that cable end supports were supposed to have infinite stiffness and to remain fixed with respect to temperature, which maximizes thermal effects. In real situations, the behaviour of end supports is more complicated and usually depends on the thermomechanical behaviour of the cable structure as a whole.

In order to check the validity of the above results, a comparison with finite element solutions is finally performed (see Fig. 7 also). Finite element solutions have been obtained from a planar Euler-Bernoulli beam model developed by the author in Ref. [START_REF] Treyssède | Prebending effects upon the vibrational modes of thermally prestressed planar beams[END_REF] for studying the effect upon vibrations of thermal stress including initial bend. The beam has been discretized into one hundred elements. Simple supports have been used as well as a negligible bending stiffness (I=5.8908e-5 m 4 , yielding a high bending stiffness parameter [START_REF] Mehrabi | Unified finite difference formulation for free vibration of cables[END_REF][START_REF] Ni | Dynamic analysis of large-diameter sagged cables taking into account flexural rigidity[END_REF] ξ = Hl 2 /EI=178). As can be observed in Fig. 7, results obtained with the proposed cable model are in quite good agreement with numerical results.

Conclusion

In this paper, the analytical model of Irvine has been extended to thermoelasticity. The effect of thermal stress typically due to climatic change on cable free vibrations has been investigated. Some general two-dimensional contour plots have been presented for the relative change in natural frequencies due to temperature with respect to the only two independent parameters of the problem. Depending on the values of λ 2 and θ, thermal stress is likely to have a quite significant effect on natural frequencies as well as modeshapes. For in-plane symmetric modes, it has been shown that a positive (resp. negative) temperature change does not necessarily decrease (resp. increase) natural frequencies and that thermoelastic effects can produce sag modifications that can lead to cross-over. As a consequence, the thermoelastic behaviour of cables are likely to affect the robustness of vibration based methods in SHM. Future studies should deal with the influence of bending stiffness, which cannot generally be neglected for higher order modes or large-diameter cables.

2 )Fig. 2 .

 22 Fig. 2. Irvine's natural frequencies vs. sag-extensibility parameter for the first three symmetric modes (solid lines) and the first two antisymmetric modes (dashed lines). Out-of-plane modes are not shown.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Relative change in natural frequency (x100) for out-of-plane and in-plane antisymmetric modes.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Relative change in natural frequency (x100) for the 2nd in-plane symmetric mode.

Fig. 7 .

 7 Fig. 7. Natural frequencies vs. temperature change for the first three symmetric modes (solid lines) and the first two antisymmetric modes (dashed lines). Out-ofplane modes not shown for clarity. Insets: first symmetric modeshapes for -40, 0 and +40 K. x-mark: finite element solution


Sylvie Chambon
email: chambon@ifsttar.fr

Jean-Marc Moliard
email: moliard@ifsttar.fr

Automatic Road Pavement Assessment with Image Processing: Review and Comparison

In the field of noninvasive sensing techniques for civil infrastructures monitoring, this paper addresses the problem of automatic crack detection by automatic analysis of optical images of the surface of the French national roads. The first contribution is a state of the art of the image processing tools applied to civil engineering. Second, we describe the proposed method to detect fine defects in pavement surface. This approach is based on a multi-scale extraction and a Markovian segmentation. Third, an evaluation and comparison protocol which has been designed for evaluating this difficult task -the road pavement crack detection -is introduced. Finally, the proposed method is validated, analysed and compared to a detection approach based on morphological tools.

Introduction

The evaluation of road quality is an important task in many countries, like, for example, in France, where the national roads are inspected each three years in order to estimate the needed reparations and constructions. To estimate the quality, these aspects can be taken into account: the adherence, the micro-texture, the macro-texture and the surface degradations. Before 1980, all these inspections were accomplished manually, but, it can be automated with noninvasive techniques, like, image processing. To be more comfortable, less dangerous for employees and users of the road but also more efficient and less expensive, many systems have been proposed, based on ground penetrating radar [START_REF] Forest | Non destructive crack depth measurements with ground penetrating radar[END_REF] or laser system [START_REF] Payne | Automating road surface analysis[END_REF]. However, for noninvasive evaluation of surface degradations, the recent research results seem more promising with optical image processing approaches for these reasons [START_REF] Schmidt | Automated pavement cracking assessment equipment -state of the art[END_REF]:

(1) The acquisition systems based on optical devices are easier to design and to use than other kinds of systems (they are less sensitive to mouvement or vibrations).

(2) They also allow a dense acquisition (each millimiter), i.e. the acquisition can be realized for the entire road surface, whereas for the other systems, like laser, the measurements are available every 4 millimeters at normal speed (90 km/h) 1 . 1 It has been determined from the most recent systems.

(3) As the acquisition system is more dense, the measurement of the defects is more precise than with other systems.

(4) Even if the images are not always well contrasted, they are more contrasted than the images/signals that can be given by other devices, i.e. the ratio between noise and signal is greater with optical sensor than with other kind of sensors.

Nowadays, many acquisition systems are available [4,[START_REF] Schmidt | Automated pavement cracking assessment equipment -state of the art[END_REF], see Table 1 (interested readers can find details about the evaluation of such systems in [START_REF] Mraz | Experimental evaluation of a pavement imaging system : Florida department of transportation's multipurpose survey vehicle[END_REF][START_REF] Sokolic | Criteria to evaluate the quality of pavement camera systems in automated evaluation vehicules[END_REF]). Moreover, to the best of our knowledge, many semiautomatic detection of road defects can be found in the literature but only one is commercialized (by INO 2 ). From all the approaches proposed, it is difficult to know which ones is the most adapted to the task and what are the actual methods that are favoured, this is why, the first goal of this paper is to present a state-of-the-art of existing methods in noninvasive control based on image processing. This task of crack detection is difficult for the special case of road crack detection because it needs to detect a signal weakly represented (1.5% of the whole image) and weakly contrasted (the road possesses a texture that hides the crack). Recent methods have shown their limits: the detection contains a lot of false detections (induced by the particular texture of the road), the detection is not enough precise (we have a region of detection and not the skeleton plus the width of the crack). The main default of the existing methods is that they do not take into account the specific geometry of the crack: it is a thin and linear object. In consequence, the second aim of this work is to introduce a new method that take into account some geometric properties of the cracks.

Even if this problem is hard, and very important in the field of civil engineering, as far as we are concerned, there is no protocol for evaluating and comparing existing methods and it is difficult to known what kind of methods has to be chosen for this task. So, we think that, with the multiple methods proposed in the literature, it is important to evaluate and to compare the various methods in order to validate previous work and to identify the approaches that can be employed and/or the methods that need improvements. So, the third aspect discussed is the introduction of such a protocol.

In consequence, the objectives of this publication are, first, to give a state of the art of existing methods in noninvasive control based on image processing for estimating the quality of the road surface, second, to present our method and, third, to introduce a protocol of evaluation and comparison that allows to highlight the advantages and drawbacks of each method.

Automatic Road Crack Detection

In the literature, many papers have introduced approaches to detect thin objects in textured images, like in medical imagery, for the detection of blood vessels [START_REF] Chaudhuri | Detection of blood vessels in retinal images using two-dimensional matched filters[END_REF], satellite imagery, for road network detection [START_REF] Geman | An Active Testing Model for Tracking Roads in Satellite Images[END_REF]. Since 1990, algorithms have been proposed for semi-automatic detection of road cracks (interested readers can see [START_REF] Scheffy | Asphalt concrete fatigue crack monitoring and analysis using digital image analysis techniques[END_REF] for details about road imaging system and their limits in 1999). For the detection of cracks, three components have to be taken into account:

(1) Acquisition (see. In this paper, only the last step is studied but the choices for the two first steps are important for the success of the image treatment. Moreover, most of the references are given in the field of road quality assessment, but, some of them come from different applications, like cracks and defects in concrete (for bridges or pipelines), on ceramics or on metallic surfaces (for industrial applications). For road cracks, most of the time, these hypotheses can be exploited:

(1) Photometric hypotheses

• (H p 1 )
The crack pixels are darker than the road pixels. • (H p 2 ) The gray-level distributions of road crack and road surface are distinguishable.

(2) Geometric hypotheses

• (H g 1 ) A crack is a fine continuous object.

• (H g 2 ) A crack is a set of connected segments that have different orientations. • (H g 3 ) A crack does not have a constant width on the whole length.

(3) Photometric and geometric hypotheses

• (H pg 1 )
The points inside a crack can be considered as points of interest, from a photometric and/or geometric point of view.

These different hypotheses can be complementary, like (H p 1 ) and (H p 2 ) or (H g 1 ) and (H g 3 ), but some of them are opposite, like (H g 1 ) and (H g 2 ). The hypothesis (H pg 1 ) combines two kinds of constraint because, the definition of a point of interest (POI), that is a significant point in a scene, can be expressed both with photometric constraints (the distribution of gray levels near POI has some particularities) and geometric constraints (the point of interest can be a corner, an edge, or any kind of geometric structure).

For image processing, we enumerate semiautomatic and automatic detection approaches, and, five families can be considered, see Table 3: (1) based on histogram analysis (hyptotheses H p 1 , H p 2 and H p 2 ): the most ancient and the most popular ones. These methods use a thresholding based on histogram analysis [START_REF] Acosta | Low-Cost Video Image Processing System for Evaluating Pavement Surface Distress[END_REF][START_REF] Mahler | Pavement Distress Analysis Using Image Processing Techniques[END_REF][START_REF] Teomete | Digital image processing for pavement distress analyses[END_REF], with Gaussian hypotheses [START_REF] Koutsopoulos | Primitive-Based Classification of Pavement Cracking Images[END_REF] and/or adaptive or local thresholding [START_REF] Darwin | Automated crack indentification for cement paste[END_REF][START_REF] Elbehiery | Surface Defects Detection for Ceramic Tiles Using Image Processing and Morphological Techniques[END_REF]. These approaches are simple and not time consuming, but they also give many false detections. In fact, these methods assume that the two gray level distributions (the road pavement distribution and the crack distribution) can be separated based on a global level statistics (histogram) 3 . In Figure 1, we can see that most of the time, this hypothesis is not valid.

(2) based on mathematical morphological tools [START_REF] Bhagvati | Gaussian normalization of morphological size distribution for increasing sensitivity to texture variations and its application to pavement distress classification[END_REF][START_REF] Coster | Image analysis and mathematical morphology for civil engineering materials[END_REF][START_REF] Elbehiery | Surface Defects Detection for Ceramic Tiles Using Image Processing and Morphological Techniques[END_REF][START_REF] Ito | Accurate extraction and measurement of fine cracks from concrete block surface image[END_REF][START_REF] Iyer | A robust approach for automatic detection and segmentation of cracks in underground pipeline images[END_REF][START_REF] Oliveira | Automatic road crack segmentation using entropy and image dynamic thresholding[END_REF][START_REF] Tanaka | A Crack Detection Method in Road Surface Images Using Morphology[END_REF] (hypotheses H p 1 and H g 1 ).

An initial thresholding is needed and the results contain less false detections than methods based on histogram analysis. However, the major drawback of this kind of techniques is the strong dependence to the parameter choices.

(3) based on a learning phase in order to alleviate the problems of the two first groups of methods [START_REF] Oliveira | Supervised strategies for cracks detection in images of road pavement flexible surfaces[END_REF][START_REF] Oliveira | Identifying and retrieving distress images from road pavement surveys[END_REF] (hypotheses H p 1 and H p 2 ). Most of them are based on neural networks [START_REF] Chou | Pavement distress classification using neural networks[END_REF][START_REF] Hsu | Airport pavement distress image classification using moment invariant neural network[END_REF][START_REF] Kaseko | A neural network-based methodology for pavement crack detection and classification[END_REF]. The drawback is the learning step that can not allow a fast and fully automatic analysis.

(4) based on filtering, the most recent ones (hypotheses H p 1 , H g 1 and H g 3 ). Using edge detections by fixed scale filtering is not adapted to the task of the detection of road cracks because the width of the crack is not constant and this is why many methods are based on wavelet detections [START_REF] Abdel-Qader | Analysis of Edge-Detection Techniques for Crack Identification in Bridges[END_REF][START_REF] Boukouvalas | ASSIST: automatic system for surface inspection and sorting of tiles[END_REF][START_REF] Wang | Waveletbased pavement distress image edge detection with " À trous" algorithm[END_REF][START_REF] Zhou | Waveletbased pavement distress classification[END_REF] with adaptive filtering [START_REF] Chambon | Introduction of a wavelet transform based on 2d matched filter in a markov random field for fine structure extraction: Application on road crack detection[END_REF][START_REF] Chaudhuri | Detection of blood vessels in retinal images using two-dimensional matched filters[END_REF][START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF] (these approaches will be detailed in the § 4), contourlets [START_REF] Ma | Pavement distress detection based on nonsubsampled contourlet transform[END_REF], Gabor filters [START_REF] Abas | Classification of painting cracks for content-based analysis[END_REF], Finite Impulse Response filter (FIR) [START_REF] Delagnes | A Markov random field for rectilinear structure extraction in pavement distress image analysis[END_REF] and methods using models based on partial differential equations (PDE) [START_REF] Augereau | A Differential Approach for Fissures Detection on Road Surface Images[END_REF][START_REF] Nishimura | A boundary integral equation for an inverse problem related to crack detection[END_REF]. Some techniques also use auto-correlation filtering [START_REF] Lee | New Crack-Imaging Procedure Using Spatial Autocorrelation Func-tion[END_REF][START_REF] Shirataki | A study of road crack detection by image processing[END_REF] (similarity measures are estimated between the targets that simulate cracks and targets of the original image). An other kind of algorithms is based on texture analysis [START_REF] Petrou | Automatic surface crack detection on textured materials[END_REF][START_REF] Song | Wigner based crack detection in textured images[END_REF] (the crack is considered as a noise inside a texture).

(5) based on an analysis of a model [START_REF] Nguyen | Automatic detection and classification of defect on raod pavement using anisotropy measure[END_REF][START_REF] Chambon | Detection of points of interest for geodesic contours: Application on road images for crack detection[END_REF] (hypotheses H p 1 , H g 1 -g 3 and H pg 1 ). These approaches are based on a local analysis versus a global analysis in order to take into account the local constraints and the global constraints of a crack. By multi-scale analysis of texture combined with an algorithm of minimal path [START_REF] Nguyen | Automatic detection and classification of defect on raod pavement using anisotropy measure[END_REF] or by local detection of point of interest combined with geodesic contours [START_REF] Chambon | Detection of points of interest for geodesic contours: Application on road images for crack detection[END_REF].

In conclusion, we can notice that: i.e. the performance of the method is dependent on the road texture.

• Old methods based on histogram studies, even those that are local, do not design correctly the problem, i.e. they do not take into account geometric characteristics of the cracks and photometric characteristics of the road pavement.

• Learning methods are efficient but the learning step is expensive (the time and the investments from the users that are not expert in image processing).

For all these reasons, even if learning methods have been used in our previous work, this paper focuses on the presentation of two methods that try to alleviate the limits of the old ones: we try to obtain a denser detection with a low rate of false detections.

Proposed and Compared Methods

Some preliminary works about methods adapted to this task (the detection of road cracks), have included experiences on a learning method. A neuron-based method has been tested [START_REF] Meignen | One Application of Neural Networks for Detection of Defects Using Video Data Bases: Identification of Road Distresses[END_REF], on the real images of size 768 × 512 presented in § 4.2. Results are interesting but learning methods are not easy to use for non-specialist in image processing, and, they cost a lot of time to the users. The main goal is to propose a system that facilitates the work of users and not a system that induces a lost of time by including a learning phase, and a maintenance each year in order to maintain the performances of the system 4 . In consequence, we have now focused our work on methods that allow automatic processing and, in particular, we present two approaches:

(1) The first, Morph, belongs to the families ( 1) and ( 2) because it combines thresholding and refinement by morphological analysis.

(2) The second, GaMM, of families ( 4) and ( 5), is based of the advantage of multi-scale analysis and local modelling of the crack.

Morph has been proposed before GaMM and is quite near the method presented in [START_REF] Tanaka | A Crack Detection Method in Road Surface Images Using Morphology[END_REF]. The contributions of this section are about GaMM, we propose a new model for the sites and the potentials used in the Markovian model. The advantages of this new method will be illustrated with qualitative and quantitative results in § 5.

Morphological method (Morph)

The chosen approach is based on hypotheses H p 1 , H g 1 and H g 3 and it follows these steps:

(1) Pre-processing of the images: to reduce the texture and increase the contrast between the road pavement and the crack;

(2) Binarization by thresholding (the threshold is different in the various variant and a local threshold can be used);

(3) Refinement by closing;

(4) Segmentation with shape analysis;

(5) Extraction of the crack characteristics.

For step 1, three variants are introduced by combining these local tools: an erosion in gray levels, a conditional median filtering, a histogram equalization, a mean filtering (these pre-processings are detailed in § 5.1). The step 4 is realised in two passes: first a connected components labeling is realized and, second, the size and the shape of each component is determined in order to remove components where shape is not similar to a crack: the shape of a crack has to look like a thin object, the width, w, and the height, h, are used for this task. More precisely, from an expert point of view, a crack is not significant if h < 50cm but, we can suppose that we manage to detect only a small part of the crack and this constraint becomes h < 7.5mm. Moreover, the mean width w min < 3mm and the maximal width w max < 6.5mm. All these thresholds are empirically set. In Figure 2, we illustrate the kind of results obtained at each step for the 3 variants. The final proposed method, named Morph, merges the 3 results (with a weighted sum and the weights are chosen with a learning phase) and refines the result by computing the closing in gray levels of the fusion result.

Adaptive Filtering and Markovian Modelling (GaMM)

More recently, our work focused on the field of wavelet decomposition. As it is difficult to chose the mother wavelet (useful for generating the wavelet family for multi-scale analysis) well adapted to the detection of road cracks, the adaptive filter theory seems convenient and, in particular, it allows to build a mother wavelet adapted to our task. First, we present the first step of the algorithm based on adaptive filtering (hypotheses H p 1 and H g 3 ) and the second on Markovian segmentation that can take into account the particular geometry of the crack (H g 2 and H g 3 ).

Algorithm

The goal of this algorithm, presented in Figure 3, is to obtain, step 1, a binarization (black pixels for background and white pixels for the crack) and a refinement of this detection by using a Markovian segmentation, step 2. Using adaptive filtering is important in order to allow the detection of non-constant width of the crack (which is realistic) (hypothesis H g 3 ). The number of scales for the adaptive filtering has to be chosen and depends on the resolution of the image. By supposing a resolution of 1 mm per pixel, by choosing 5 scales, a crack with a width from 2 mm to 1 cm can be detected. Moreover, the number of directions (for the filtering) also has to be chosen and, it seems natural to take these four directions: [0, π 4 , π 2 , 3π 4 ] that correspond to the four usual directions used for crack classification. The adaptive filtering is applied in each scale, each directions and then all the results are merged on each scale (mean of the coefficients). The results of this filtering is used to initialize the Markovian modelling used for the segmentation step.

Input

Adaptive filtering

Some details are provided in order to realize the step 1a and 1b in Figure 3. The ψ ∈ L 2 (IR 2 )5 function is a wavelet if:

IR 2 |Ψ(x)| 2 x 2 dx < ∞, with x = (i, j), ( 1
)
where Ψ is the Fourier transform of ψ. The equation (1) induces that IR 2 ψ(x)dx = 0. The wavelet family is defined for each scale s and for each position u, by :

ψ(x) = 1 2 ψ(R θ ((x -u)/s)), (2)
where R θ is a rotation of angle θ.

One of the main difficulties to apply a wavelet decomposition is the choice of the mother wavelet ψ. The adaptive filter h of s is defined by:

h = φ -1 bb s. (3)
The crack signal depends on the definition of the crack. In this paper, like in most of the papers of this domain, crack pixels correspond to black pixels surrounded by background pixels (road pixels). This is why, in [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF], a crack is a piecewise constant function f , defined for each position x ∈ IR by:

f (x) = -a If x ∈ [-T 2 , T 2 ] 0 elsewhere, (4)
where the factor a and the threshold T have to be determined. It does not correspond to a realistic representation of the crack. Because of sub-sampling, lights, orientation of the camera, the signal is more like a Gaussian function with zero mean:

f (x) = -a e -1 2 ( x σ ) 2 , ( 5
)
where a is the size of the crack and depends on σ, the deviation of the Gaussian law, i.e. a =

1 σ √ (2π)
. Consequently, the term σ allows to fix the width of the crack (like threshold T in equation ( 4)). Finally, for the step 1, h is estimated for each size of signal (determined by σ) for the 5 scales, as explaines in the beginning of § 3.2.1, and φ bb is interpolated in order to have the same size. Then the filter is transformed by rotation in order to cover the 4 orientations.

Segmentation

The goal of this part is to extract shapes, i.e. cracks, using the detection maps estimated at the first stage of the algorithm (step 2a of the algorithm of the Figure 3). For the first step of segmentation (initialization), the sites are of size 3 × 3, consequently, a regular grid is considered in the image. In [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF], four configurations are possible and represented in Figure 4 (the part inside the rectangle with low gray levels). The initialization of the sites is based on the configuration that maximizes the coefficients obtained with the adaptive filtering. More formally, if we denoted

γ 2,0 , γ 2, π 4 , γ 2, π 2 and γ 2, 3π 4
, the four configurations, the best configuration γ best is:

γ best = argmax α∈[0,..., 3π 4 ] m 2,α , (6)
where m 2,α is the mean of the coefficients on the considered configuration γ 2,α . These four configurations do not represent all the possibilities and are not realistic configurations. In fact, all these four configurations are centered, whereas, it is possible to have some non-centered configurations. Consequently, we use the set of sixteen configurations illustrated in Figure 4 (all the presented sites). By modifying the number of configurations, we need to adapt the initialization of sites and equation ( 6) becomes:

γ best = argmax i∈[0;4], α∈[0,..., 3π 4
]

m i,α , (7)
where m i,α is the mean of the coefficients on the considered configuration γ i,α . The image is considered as a finite set of sites denoted S = {s 1 , . . . , s N }.

For each site, the neighborhood is defined by:

V s = s ′ |s / ∈ V s ′ & s ′ ∈ V s ⇒ s ∈ V s ′ . A clique c
is defined as a subset of sites in S whose every pair of distinct sites are neighbors. These random fields are considered:

(1) The observation field Y = {y s } with s ∈ S. Here, y s is the mean of the coefficients on the site.

γ 3, 3ψ 4 γ 0, 3π 4 γ 1, 3π 4 γ 2, 3π 4
γ 4, 3π 4 γ 1,0 γ 3,0 γ 2,0 γ 0, π 4 γ 2, π 2 γ 3, π 2 γ 1, π 2 γ 3, π 4 γ 1, π 4 γ 2, π 4 γ 4, π 4
Figure 4: The sixteen configurations in order to improve the modeling of sites -The four initial configurations proposed in [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF] are in the bold rectangle, the sites are represented by the clearer gray levels, and for the proposed configurations the sites are represented by the darker gray levels.

(2) The descriptor field L = {l s } with s ∈ S. If there is a crack l s = 1 elsewhere l s = 0.

At each iteration, a global cost, or a sum of potentials, that depends on the values of the sites and the links between neighborhoods, is updated. This global cost takes into account the coefficients of the sites (computed from the coefficients estimated during the first part of the algorithm: adaptive filtering) and the configurations of each site and its neighbor sites (the 8 neighbors). More formally, the global cost is the sum of all the potential functions of the sites. This potential function contains two terms:

u s = α 1 u 1 (s) + (1 -α 1 ) ∑ s ′ ∈V s u 2 (s, s ′ ), (8)
The first term, u 1 , corresponds to the data term, and it evaluates how a site is similar to a crack from a photometric point of view (hypotheses H p 1 and H p 2 ). This term is based on the results given by the adaptive filtering. The second term, u 2 , represents the constraints induced by the neighbors of the site. More precisely, it estimates the consistency between a site and each neighbor site and it takes into account the geometric hypotheses H g 2 and H g 3 . The choice of the value α 1 depends on the importance of each part of the equation ( 8) and it will be discussed in § 5.1.1.

The function u 1 is given by:

u 1 (y s , l s = 1) = e ξ 1 (k-y s ) 2 If y s ≥ k 1 elsewhere and u 1 (y s , l s = 0) = e ξ 2 (y s -k) 2 If y s < k 1 elsewhere, (9)
The parameters ξ 1 , ξ 2 and k have to be fixed 6 . For the definition of u 2 , we have to determine the number of cliques. In [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF], 4 cliques are possible and the 8-connexity is considered. The potential function proposed in the precedent work only considers the difference of orientations between two neighborhoods and not the position between the two sites of the clique, see Table 4. Some cases are not penalized with the

γ 2,0 γ 2, π 4 γ 2, π 2 γ 2, 3π 4
γ 2,0 β 1 β 2 β 3 β 2 γ 2, π 4 β 2 β 1 β 2 β 3 γ 2, π 2 β 3 β 2 β 1 β 2 γ 2, 3π 4 β 2 β 3 β 2 β 1
β 1 = -2, β 2 = -1 and β 3 = 2.
old configuration. For example, these two unfavorable cases are not penalized:

• two sites with the same orientation but with no connection between them ; • two sites with the same orientation but their position makes them parallel. This is why, with the sixteen configurations that are presented in Figure 4, the potential has to take into account the differences of orientations between two sites (there are 16 × 16 possibilities) and the position of the two sites (there is 8 possibilities because we consider the 8 neighbors). Consequently, the new potential function u 2 follows these two important rules:

(R 1 ) The lower the difference of orientations between two sites, the lower the potential.

(R 2 ) The lower the distance between two sites, the lower the potential (in this case, the distance implies the minimal distance between the extremities of the two segments). 6 The choice of k is related to the maximal number of pixels that belong to a crack (it depends on the resolution of images and hypothesis about the size and configuration of cracks). We have chosen k in order to consider at most 5% of the image as a crack. Moreover, our experiments have brought us to take ξ 1 =ξ 2 =100.

s 0 s 1 s 2 u 2 s 3 s 0 s 3 s 1 s 2 s 3 s 0 s 1 s 2 θ 0 = θ 1 = θ 2 = π 4 , θ 3 = π 2 θ e1 = 0, θ e2 = π 4 , θ e3 = π 2
For s 1 :

u 2 = 0 -1 + 0 For s 2 : u 2 = π 4 -1 + 0 For s 3 : u 2 = π 2 -0 + 2 θ 0 = θ 1 = θ 2 = θ 3 = π 2 θ e1 = 0, θ e2 = θ e3 = π 2
For s 1 :

u 2 = 0 -1 + 0 For s 2 : u 2 = π -3 + 0 For s 3 : u 2 = 0 -0 + 2
Figure 5: Examples of the function u 2 -These two examples of sites with their respective neighbors show the behavior of the potential u 2 with the two considered aspects: orientation and distance. In example (a), with the help of the orientation term, the configuration s 3 is penalized and s 2 is less penalized than s 3 . In example (b), with the help of the two terms about the distance, the site s 3 is penalized, compared to s 1 . On the contrary, the particular case of s 2 is favorable and it compensates the penalty given by the orientations.

More formally, if:

• d denotes the Euclidean distance between the two closest extremities of the sites, with d ∈ [0, d max ]7 ;

• θ 1 and θ 2 are the orientations of respectively s = {p i } i=1..N s , and s ′ = {s ′ j } j=1..N p ′ where, p i , respectively p ′ i , is the pixel i of the N s , respectively N ′ s , pixels that composes the site s, respectively s ′ ;

• θ e is the angle between the two sites; the u 2 function is defined by:

u 2 (s ′ , s) = α 2 |2θ e -θ 1 -θ 2 | 2π + (1 -α 2 ) J(NbC) min i, j (d(s i ,s ′ j )) d max -NbC 3 .
(10) where NbC indicates the number of connected pixels between the two sites s and s ′ and J(x) equals 1 if

x = 0 and 0 elsewhere. The first term is induced by the rule about the orientations, (R 1 ). This term equals zero when the sites have the same orientation and this orientation is the same as the orientation between the sites, i.e. θ e = θ 2 = θ 1 . This first term penalizes the configurations where the sites do not have the same orientation but also the particular case where they are parallel, see example (a) in Figure 5. The second term and the third term express the rule (R 2 ) about the distances. Two aspects have to be distinguished: the number of connected pixels, when the sites are connected, and, on the contrary, i.e. when the sites are not connected, the distance between the sites. It allows to give low influence at disconnected sites and also to increase the cost of sites that are parallel but connected, see example (b) in Figure 5. To study the influence of all these terms, the equation has been normalized and the different terms have been weighted (using α 2 , the choices for α 2 will be discussed in § 5.1.1).

Evaluation protocol

For the evaluation of image processing methods, nowadays, there is no protocol used by the community of road pavement analysis. However, all the approaches for characterizing the quality of the road advise to take into account the severity of the defect. For this, we need to know precisely the size, the width and the location of the cracks. This is why, now, it seems important to precisely evaluate the performances of all these automatic existing methods. For designing this protocol, we have to determine which kind of images have to be tested and which kind of criteria of efficiency can be used. In the community, to the best of our knowledge, there is no reference images and most of the time, the evaluation of the performances is qualitative.

For the evaluation of automatic crack detection, to the best of our knowledge, no evaluation and comparison protocol has been proposed in the community. However, in all the countries, for estimating the quality of the road surface, it is important to know exactly the size and the width of defects, i.e. to detect precisely the defect. This is why, it seems important to characterize quantitatively the performances of the methods. For building this kind of protocol, it is necessary, first, to choose the tested images, second, to choose how to build reference segmentations, and, third, to determine the criteria used for the quantitative analysis. For the reference segmentations, two approaches can be used to estimate them: (1) To compute synthetic images with synthetic defects -In this case, the exact position of the defects is known and these reference segmentations can be considered as ground truth.

(2) To propose reliable segmentations of real images -It supposes that we are able to provide a segmentation that is reliable enough to be employed as a reference. For evaluation, these segmentation can be called "pseudo-ground truth". The two solutions are studied, and, this is why we have to explain how the manual segmentations (that are our references) are done. Before, we briefly describe the acquisition system. In some cases, the sensors do not have the same settings and the global illuminations are different, so, it can generate some "false cracks". This aspect has been easily taken into account in a pre-processing step by eliminating the junction area in the region of interest.

Acquisition

The acquisition system used for the dataset of our experiments is described in Figure 6. It contains 4 video cameras with 3 sensors in gray levels in the backside of the car and 1 color one in front of the car. The first camera is needed to determine the environment conditions (weather, location, traffic) whereas the three other ones are used for the crack detection.

The resolution of this one is smaller than the 3 others and moreover, the optical axis is not perpendicular to the road surfaces, on the contrary of the 3 others. The 3 cameras have been physically synchronized directly during the acquisition. To be independent of the illumination problems, nine stroboscopic lights have been added. The position of the lights is perpendicular to the road plane and distant of 1 meter from the surface. The light power has been chosen in order to not deteriorate the visualisation of the road pavement and of the defects.

Reference images

The most difficult is to propose images with a reference segmentation. On the first hand, we introduce synthetic images with a simulated crack (the size of these images are 256 × 256) 8 As shown in Figure 7, the result is not realistic enough. It does not seem realistic because the contrast is too important between the road and the crack. Moreover, the interuptions of the crack, the changes in the direction, the presence of many paths, etc, in the default, are no simulated. In order to be more realistic, it seems that we have to design and to implement a complex heuristic to simulate the crack and it is too much effort for finally having only a synthetic default. This is why, on the second hand, we have simulated different defects on real images that previously contain no defect (the size of these images are 768 × 512 and 1920 × 480). The result is more realistic but the shape and the photometric aspect of the cracks (which are randomly chosen) does not seem realistic enough. This is why, it appears important to propose a set of real images (size 768 × 512 and 1920 × 480) with manual segmentations that are reliable enough to be considered as reference segmentations. To summarize, the two first kinds of images allow to propose an exact evaluation and to illustrate theoretically the behavior of the method whereas the last kind of images allows to validate the work on real images with a pseudo-ground truth.

Reference segmentations

For real images, we briefly explain how the manual segmentations are validated. Four experts have man- 8 The road is a random texture, i.e. each intensity is randomly chosen by supposing a uniform distribution of intensities in [0; 255]. Then, the user gives the position of the the length and the orientation (vertical, horizontal or oblique) of the crack. The crack points are built by randomly selecting the next point in the neighborhood and the intensity in [0; 100].

Synthetic image Ground truth

Real image + simulated defect Ground truth

Real image manually segmented Pseudo ground truth ually segmented the images with the same tools 9 and in the same conditions. Then, the four segmentations are merged, following these rules:

(1) A pixel marked as a crack by more than two experts is considered as a crack pixel;

(2) Every pixel marked as a crack and next to a pixel kept by step (1) or ( 2) is also considered as a crack.

The second rule is iterative and stops when no pixel is added. Then, the result is dilated with a squared structuring element of size 3 × 3. To evaluate the reliability of the reference segmentations, we estimate, first, the percentage of covering between each operator, and, second, the mean distance, D, between each pixel (detected by only one expert and not kept in the 9 We use a "home-made" software that proposes an interface that helps the person to segment the default. The principle is that the user has to select points on the crack. These points have to be close enough (from 5 to 20 pixels of distance). Then, the path between two close points is automatically detected by using a simple heuristic: the path that minimizes the mean intensity is selected. The interface is complete enough to allow the displacement of the points, the removing of some points, the removing of some cracks. The user can also select the width of the path (crack). Some filters are also proposed to improve the contrast between the crack and the road in order to help the user. reference image) and the reference segmentation.

Table 5 shows some results for 5 of the 42 images manually segmented. We have distinguished 5 families: the first one contains images acquired in static whereas the four other ones are acquired in dynamic. Moreover, we have 4 different kinds of road pavement acquired in dynamic. The 10 images have been taken in order to show results on each of these families. We can notice that the first 4 images are the most reliable because the mean error is less than 2 pixels. The precision of these results is satisfactory. On the contrary, the last 6 images show the important variabilities between operators and how it is difficult to extract a segmentation for these images, and, in particular, in the image 936, where the error is due to a bad interpretation of one of the four operators who finds a defect that does not exist.

IMAGES F (%) 2 (%) 3 (%) 4 (%) S (%) D (PIX) The comparison of the 4 manual segmentations for estimating the final reference segmentations -For each image, are presented: the percentage of pixels in the whole image that are preserved as crack pixels in the final reference segmentation (F), the percentage of recovering between 2, 3 and 4 manual segmentations and the sum of these 3 percentages (S). For all the crack pixels that are not preserved in the final reference segmentation, the mean distance to this segmentation is given (D).

By analyzing the results for the criterion D, presented in Table 5, we can classify the 42 tested images in 3 categories, i.e. images with:

(1) A reliable segmentation: The criterion D < T r . It means that all the operators have built segmentations that are quite near to each other.

(2) A segmentation that is moderately reliable : The criterion T r ≤ D < T a . It means that some parts of the crack are not easy to segment and there are local errors.

(3) An ambiguous segmentation: The criterion D ≥ T r . It clearly shows that the images are difficult to segment and in most of the cases, it means that some parts are detected as a crack whereas they are not and reversely.

The threshold have been empirically chosen and T r = 2, T a = 4. In Figure 8, first, we present the mean distance, D i , i ∈ {1, 2, 3, 4}, between the final reference segmentation, S r , and each manual segmentation, S i (obtained for the four operators) and, second, the criterion D for each real images of our protocol. The first graph illustrates how it is important to combine the four manual segmentations instead of using only one manual segmentation. Indeed, we can notice that each operator, alternately, gives an interpretation that is different from the three others. The second graph explains how the thresholds are chosen for determining the detections that are "accepted" for the evaluation, see § 4.4 for explanations about accepted detections.

All these remarks are illustrated in Figure 9. Overall, the four segmentations are near to each other and if the segmentations are combined it permits to detect the width of the crack. However, these examples also present the difficulties of this task: areas where the cracks are less visible, the texture elements that have the same size and/or the same gray levels as the crack pixels. Thus, in some cases, one operator extends the crack or gives a different shape. In some extreme cases, the operator can even confound a crack and an other object of the scene (a piece of wood, for example). In another way, these last examples highlight the interest in combining different segmentations in order to obtain reference segmentations as reliable as possible.

Criteria of efficiency

In this section, we introduce how the reference segmentation and the estimated segmentation are compared. In Figure 10, we present common evaluation criteria that are used for segmentation evaluation:

(1) the percentage of correct detections (true positives) (TP) ;

(2) the percentage of false positives (FP) ;

(3) the percentage of false negatives (FN) and

(4) the similarity coefficient (DICE).

This last criterion seems to be the most significant because it evaluates the ratio between the FP and the FN, and it well resumes the results of all the criteria. Moreover, it directly expresses what is important to evaluate: how the evaluated method can reduce errors of detection whilst increasing the density of good detection.

For real images, the detections that are "accepted" have been added in order to tolerate a small error on the localization of crack pixels. This criterion is needed because perfect detection seems, for the moment, difficult to reach, see the results in Table 5, that illustrate this aspect. In consequence, these accepted pixels have been included in the estimation of the similarity coefficient or DICE. The threshold for accepted Figure 8: The variations between each manual segmentations that are used for building pseudo-ground truth -The first graph represents for each operator (one curve for one operator) and each image (x-coordinates), cf. Figure 17 for the corresponding images, the mean distance, D i , i ∈ {1, 2, 3, 4}, whereas, the second graph, presents D, cf. s § 4.3 and Table 5. This graph allows us to distinguish the different categories of images (red axes): reliable (D < T r ), moderately reliable (T r ≤ D < T a ), ambiguous (T a ≤ D). In the two graphs, the purple axes represent the five different samples of images (each sample corresponds to a kind of road pavement). The first ones were acquired with a static system whereas the four others were acquired with a dynamic system. pixels equals 0 for synthetic images whereas it depends on the mean distances, see D in table 5, for the real images.

Experimental results

In this section, two aspects are studied:

(1) the evaluation of the method based on an adaptive filtering and a Markovian modelling in order to characterize its behavior, to estimate the best parameters and to determine the best variant ;

(2) the comparison to the Morph method.

Adaptive filtering and Markovian Modelling

We want to determine, first, how to fix the different parameters, second, the pre-processing steps that are necessary, and, finally, which variant is the most efficient. In consequence, these points have been studied:

• Parameter values -The weights α 1 , equation ( 8), and α 2 , equation [START_REF] Chambon | Detection of points of interest for geodesic contours: Application on road images for crack detection[END_REF], are tested from 0 to 1 with a step of 0.1.

• Pre-processings -pre-processings have been experimented to reduce noises induced by texture, to increase the contrast of the defect and, to reduce the light halo in some images:

(1) Threshold -This pre-processing has been proposed in order to reduce the light halo in the six last images presented in Figure 7 and in all the im-

Influence of parameters

Among the results, two conclusions can be done:

(1) For each variant and each pre-processing, the weights between the term for adaptive filtering and the term for the Markovian modelling should be the same, see equation ( 8), i.e. α 1 = 0.5. However, when more weight is given to adaptive filtering, the quality of the results is lower than when more weight is given to the Markovian segmentation. It means that in this kind of application, the geometric information is more reliable than the photometric information. It seems coherent with the difficulties of the acquisition.

(2) For the Markovian modelling, we have noticed that the results are the best when the weights are the same between the orientation term and the distance term, see equation [START_REF] Chambon | Detection of points of interest for geodesic contours: Application on road images for crack detection[END_REF], i.e. α 2 = 0.5. However, better results are obtained when the weight of the orientation is greater than the distance one instead of the reverse. It implies that the orientation characteristics are more reliable than the distance ones and this remark is coherent with the fact that cracks present strong spatial constraints. Moreover, it is also linked with the difficulties induced by the acquisition (the lighting system makes the photometric information less reliable).

Pre-processing

These tests have been done with real images, because, the synthetic images do not need pre-processings. The results are given by: Init Gaus InMM GaMM Restoration Restoration Threshold Erosion However, for the four first images (acquired with lighting conditions more comfortable than the lighting conditions of the next 6 ones), the pre-processing is not significant for increasing the quality of the results. Moreover, with the new Markovian modelling, the pre-processing step does not significantly increase the quality of the results.

Variants

The results are presented for two different cases:

(1) with synthetic images and

(2) with real images.

For the first category, the ground truth is available whereas for the second category, a pseudo-ground truth is used and the detections which are accepted are taken into account in the evaluation, i.e. a threshold is applied to the distance between the segmentation estimated by the evaluated method and the pseudoground truth segmentation. The thresholds applied on the distance for the accepted detections are determined with the results given in Table 5, column D. In Figure 11, the evolution of the similarity coefficients, or DICE, for the 11 synthetic images, 11.(a) and 10 of the real images is presented, 11.(b). With synthetic images, the method GaMM is clearly the best for most of the images. However, for one image (the fifth), the results are worse than the results of the method Gaus but they are still correct (DICE=0.72).

On the contrary, for the most difficult images (the 3 first ones that contain a real road background), the method GaMM obtains acceptable results (DICE > 0.5) whereas the other methods are not efficient at all. Illustrations are given in Figures 12 and13: they show how the method GaMM can reduce false detections.

Image

Ground truth Init

Gaus InMM GaMM

Figure 12: The segmentation results on some synthetic images presented in Figure 7 -These are the results obtained with the four variants and these illustrations show how the method GaMM gives the clearest result. We can also notice the good results of the method InMM.

Results and comparison with Morph

Finally, we have compared the results of GaMM on each of the complementary dataset (32 images) with Morph. The mean DICE is 0.6 with GaMM whereas it is 0.49 with Morph, see Figure 14. It shows how GaMM can outperform Morph. However, if we compare image per image, the results show that in 50% of the cases GaMM is the best, see illustrations of these results in Figures 15 and16. More precisely, GaMM seems more efficient with ambiguous images, whereas Morph is the best with reliable images. Finally, we can also precise the execution time for the two methods: about 1 minute for GaMM and 5 seconds for Morph with a processor intel R core TM 2 duo of 2 GHz). These execution times give only some indications because the implementation, in particular for GaMM, has not been optimized.

Conclusions

In conclusion, this paper gives a review about image processing methods for the crack detection of road pavement. It can help the researchers who want to choose and to adapt an auscultation method to the constraints of the transport structure that is studied (it depends on the quality of the surface, the needs Init Gaus InMM GaMM 13: The segmentation results on some real images -These are the results obtained with the real images presented in Figure 7. The method InMM obtains the clearest (i.e. with less false detection) but we can also noticed the good quality of the detection map with the method GaMM.

of the auscultation). Moreover, a new method for the detection of road cracks has been introduced and we have presented a new evaluation and comparison protocol for automatic detection of road cracks. As far as we are concerned, we proposed real images with ground truth for the first time in the community. The new method, GaMM, has been validated by the proposed protocol and compared to a previous one, Morph. This evaluation shows the complementarity of the two methods: the Morph method obtains more true positives than the GaMM method whereas this ones reduces the percentage of false positives.

Our first improvements of this work will focus on the evaluation and comparison protocol. We want to increase our data set by taking into account the different qualities of road surface or road texture (because for the moment, each proposed method seems very dependent on the quality of the road texture). In a second step, our future work will include new experiments about the acquisition system. Indeed, the acquisitions and the results obtained with the acquisition system presented have shown its limits, for example, in Figure 7, some parts of the crack are not "visible". It comes from the fact that to highlight the crack, it depends on the orientation of the lights and of the sen-Figure 14: The comparison of the similarity coefficients between GaMM and Morph -The dotted lines illustrate the five sets of tested images. For the first set that corresponds to real images with no illumination problems, the results are mixed whereas for the four other sets, GaMM is the best. The mean of this criterion is 0.6 (variance = 0.0257) for GaMM whereas it is 0.49 (variance = 0.0750) for Morph. However, this method has one step of characterization of the cracks (not introduced in GaMM) and this step can remove cracks that do not respect the characteristics of a cracks (in length, size and shape). This step contributes to reduce errors, but, in some cases, it decreases the performances of the detection, compared to GaMM. sors. Using one single sensor and one light always in the same position/orientation, we can sometimes miss some defaults in the acquisition. So, it seems important to study other kinds of system to improve the quality of the automatic treatments.

Third, we want to improve the GaMM method, to begin, by adding the extraction of the crack characteristics, like in Morph.

( 2 ) 3 )

 23 Storage and (Image processing.

Figure 1 :

 1 Figure 1: Examples of histograms (second line) of images (first line) with cracks -These images present only one mode in their histograms and it is impossible to separate the gray level distribution of cracks from the gray level distribution of the road pavement.

1 )

 1 For each scale do 1a) For each direction do Estimate Adaptive Filter (AF) 1b) Merge AF in all the directions 2) For each scale do 2a) Initialization of the sites (Markov) 2b) While not (stop condition) do Updating of the sites 3) Fusion of the results on each scale

Figure 3 :

 3 Figure 3: The studied algorithm for the method based on adaptive filtering and segmentation by Markovian modelling -The step 1 leads to a binary image using adaptive filtering, while step 2 refines this result with a Markovian modelling.

Figure 2 :

 2 Figure2: The different steps of the method morph -The conditional filtering is applied when the gray level is higher than 40 (to prevent the removal of the crack). The last step is proposed in order to reduce false detections and to complete the detection.

Figure 6 :

 6 Figure6: The acquisition system used for the evaluation -In (a), The acquisition system is illustrated. In (c) and (d), an example of the final images is given. In (d), we can see the road that is visible in (c). The processing is done 1 meter by 1 meter, i.e. independently on each image presented in (d). The surface contains two reparations of vertical cracks. In some cases, the sensors do not have the same settings and the global illuminations are different, so, it can generate some "false cracks". This aspect has been easily taken into account in a pre-processing step by eliminating the junction area in the region of interest.

Figure 7 :

 7 Figure 7: Examples of the images tested.

Figure 11 :

 11 Figure 11: The variations of the similarity coefficients, see Figure 10, for the 4 variants -The first graph shows the results for synthetic images (the 3 first ones are obtained from real images with simulated defect) and the second graph presents the results with real images. The good performances of the methods InMM and GaMM can be noticed.

Figure 17 :

 17 Figure 17: The reference segmentations (or pseudo-ground truth) used in our evaluation protocol for real images -In each original image, the reference segmentation is highlighted in red. Visually, these manual segmentations seem correct and useful for being the reference segmentations in our comparison protocol.

Table 2

 2

for details);

Table 1 :

 1

	NAME	COUNTRY L	CLS LC	IPT
	ADDA, Automated Distress Data Acquisition (1991)	USA			-
	ACM, Automated Crack Monitor (1991)	USA			-
	SIRANO (1991)	France			-
	HARRIS, Highways Agency Road Research Information System (1999) a	UK	-		-
	ADA, Automated Distress Analyser (2007) b	USA	-		-
	AIGLE RN (2008) c	France		-
	AMAC (2004) d	France	-	-	-	-
	Profilograph and laser (2007) e	Danemark	-	-	-
	REAL, Road Excellent Automatic Logging f (1992)	Japan	-	-	-	-
	RoadCrack (2008) g	Australia		-	-	-
	ARAN, Automatic Road Analyser (1993) and then FuGro ADVantage (2003) h	Canada	-	-		-
	PAVUE j	Sweden		-		-

The quality evaluation systems of the road pavement based on optical devices -L, CLS, LC and IPT mean respectively, that the described system involves a Laser, a Controlled Lighting System, a Linear Camera and/or some Image Processing Tools. The first set corresponds to systems without image processing tools. While the second group of systems are able to provide non-dense measurement, the third group of systems can give a dense evaluation of the quality of the road.

Table 2

 2

	: The imaging systems -For each kind of configuration, the given details are: the orientation of the sensor (O.) that
	can be perpendicular or not to the road plane, the number of sensors, the presence (-) or not of a laser, the presence (-) or not
	of a overlapping images (R.) between the different acquisition. For the sensors, we distinguish 2D sensors (CCD, Charged
	Couple Device) from 1D sensor (Linear camera coupled with laser).
	HISTOGRAM	LEARNING	MORPHOLOGY	FILTERING	MODEL-BASED
	(H p 1 , H p 2 ) [3, 14, 16, 19, 21, 23, 30, 33, 36, 38, 39, 40, 41, 44, 55, 67, 72]	(H p 1 , H p 2 ) [8, 17, 18, 20, 32, 37, 45, 50, 51, 57, 58]	(H p 1 , H g 1 ) [6, 15, 25, 26, 34, 35, 47, 56, 64, 66]	(H p 1 , H g 1 -g 3 ) [1, 2, 5, 7, 13, 22, 24, 29, 42, 43, 49, 52, 54, 61, 63, 65, 69, 71, 73]	(H p 1 , H g 1 -g 3 , H pg 1 ) [9, 10, 11, 12, 48, 68, 70]

Table 3 :

 3 The classification of crack detection (or thin objects in textured images) into five different families. For each family, the hypotheses that are employed are specified.

Table 4 :

 4 The function u 2 used in [65] -This table presents the values u 2 (s ′ , s) for the sites in low gray levels in Figure 4. In our experiments, like the authors, we have chosen

Table 5 :

 5

http://www.ino.ca/fr-CA/Realisations/Description/ project-p/systeme-laser-mesure-fissures.html

This separation would be possible based on local statistics around the crack.

This maintenance is necessary because conditions and systems of acquisition can changed every year, and the road pavement also evolves.

L 2 is the square integrable space.

As the sites are of size 3 × 3,d max = 5 √ 2)

a

MODERATELY RELIABLE SEGMENTATIONS

AMBIGUOUS SEGMENTATIONS Figure 9: The levels of difficulties of the tested images -On these images, we present the four examples of manual segmentations. The codes are: red (light and dark), blue, green, for each of the four operators. The parts in yellow correspond to the parts of the cracks detected by more than one operator. There are two examples per category of segmentations. To better visualize, only a part of the images is shown. In the ambiguous images, we clearly see the mistake of the operator in dark red.

ages of the four last categories in Figure 8. Each pixel lower than a given threshold is replaced by the local average of the gray levels.

(2) Smoothing -A mean filter of size 3 × 3 is applied to reduce the granularity of the texture. (3) Erosion -An erosion (in gray level) with a square structuring element of size 3 × 3 is also applied to reduce the granularity of the texture. In order to preserve the crack signal, each pixel under a given threshold is not filtered 10 .

• Algorithm variants -Four variants are compared:

(1) Init -This is the initial method proposed in [START_REF] Subirats | Automation of pavement surface crack detection with a matched filtering to define the mother wavelet function used[END_REF].

(2) Gaus -This variant supposes that the distribution of the gray levels inside a crack follows a Gaussian function, see § 3. • Comparison -We have compared this method with the method based on morphological tools and that is quite similar to [START_REF] Tanaka | A Crack Detection Method in Road Surface Images Using Morphology[END_REF], Morph.

Synthetic images

Ground truth

Morph GaMM

Real images Pseudo ground truth