Fabrizio Bianchi

Tien-Cuong Dinh

EVERY COMPLEX H ÉNON MAP IS EXPONENTIALLY MIXING OF ALL ORDERS AND SATISFIES THE CLT

Keywords: 2010 Mathematics Subject Classification. 37F80 (primary), 32U05, 32H50, 37A25, 60F05 (secondary) Complex Hénon maps, Exponential mixing of all orders, Central Limit Theorem

We show that the measure of maximal entropy of every complex Hénon map is exponentially mixing of all orders for Hölder observables. As a consequence, the Central Limit Theorem holds for all Hölder observables.

Notation. The pairing •, • is used for the integral of a function with respect to a measure or more generally the value of a current at a test form. By (p, p)-currents we mean currents of bi-digree (p, p). Given k ≥ 1, we denote by ω FS the Fubini-Study form on P k = P k (C). The mass of a positive closed (p, p)-current R on P k is equal to R, ω k-p FS and is denoted by R . The notations and stand for inequalities up to a multiplicative constant. If R and S are two real currents of the same bi-degree, we write |R| ≤ S when S ± R ≥ 0. Observe that this forces S to be positive.

Introduction

Hénon maps are among the most studied dynamical systems that exhibit interesting chaotic behaviour. They were introduced by Michel Hénon in the real setting as a simplified model of the Poincaré section for the Lorenz model, see, e.g., [START_REF] Benedicks | The dynamics of the Hénon map[END_REF][START_REF] Hénon | A two-dimensional mapping with a strange attractor[END_REF]. Hénon maps are also actively studied in the complex setting, where complex analysis offers additional powerful tools. The reader can find in the work of Bedford, Fornaess, Lyubich, Sibony, Smillie, and the second author fundamental dynamical properties of these systems, see [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . III. Ergodicity, exponents and entropy of the equilibrium measure[END_REF][START_REF] Dinh | Rigidity of Julia sets for Hénon type maps[END_REF][START_REF] Fornaess | Dynamics in several complex variables[END_REF][START_REF] Fornaess | Complex Hénon mappings in C 2 and Fatou-Bieberbach domains[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF] and the references therein. It is shown in [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF] that the measure of maximal entropy of such systems is Bernoulli. In particular, it is mixing of all orders. On the other hand, the control of the speed of mixing for general dynamical systems is a challenging problem, and usually one can obtain it only under strong hyperbolicity assumptions on the system. The main goal of this work is prove that the measure of maximal entropy of any complex Hénon map is exponentially mixing of all orders with respect to Hölder observables. As a consequence, we also solve a longstanding question proving the Central Limit Theorem for all Hölder observables with respect to the maximal entropy measures of complex Hénon maps.

Let us first recall the following general definition. Definition 1.1. Let (X, f ) be a dynamical system and ν an f -invariant measure. Let (E, • E ) be a normed space of real functions on X with • L p (ν)

• E for all 1 ≤ p < ∞. We say that ν is exponentially mixing of order κ ∈ N * for observables in E if there exist constants C κ > 0 and 0 < θ κ < 1 such that, for all g 0 , . . . , g κ in E and integers 0 =:

n 0 ≤ n 1 ≤ • • • ≤ n κ , we have ν, g 0 (g 1 • f n 1 ) . . . (g κ • f nκ ) - κ j=0 ν, g j ≤ C κ • κ j=0 g j E • θ min 0≤j≤κ-1 (n j+1 -n j ) κ .
We say that ν is exponentially mixing of all orders for observables in E if it is exponentially mixing of order κ for every κ ∈ N.

A recent major result by Dolgopyat, Kanigowski, and Rodriguez-Hertz [START_REF] Dolgopyat | Exponential mixing implies Bernoulli[END_REF] ensures that, under suitable assumptions on the system, the exponential mixing of order 1 implies that the system is Bernoulli. In particular, it implies the mixing of all orders. (with no control on the rate of decay of correlation). It is a main open question whether the exponential mixing of order 1 implies the exponential mixing of all orders, see for instance [START_REF] Dolgopyat | Exponential mixing implies Bernoulli[END_REF]Question 1.5].

Let now f be a complex Hénon map on C 2 . It is a polynomial diffeomorphism of C 2 . We can associate to f its unique measure of maximal entropy µ [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . III. Ergodicity, exponents and entropy of the equilibrium measure[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF], see Section 2 for details. It was established by the second author in [START_REF] Dinh | Decay of correlations for Hénon maps[END_REF] that such measure if exponential mixing of order 1 for Hölder observables, see also Vigny [START_REF] Vigny | Exponential decay of correlations for generic birational maps of P k[END_REF] and Wu [START_REF] Wu | Exponential mixing property for Hénon-Sibony maps of C k[END_REF]. Similar results were obtained by Liverani [START_REF] Liverani | Decay of correlations[END_REF] in the case of uniformly hyperbolic diffeomorphisms and Dolgopyat [START_REF] Dolgopyat | On Decay of Correlations in Anosov Flows[END_REF] for Anosov flows.

Theorem 1.2. Let f be a complex Hénon map and µ its measure of maximal entropy. Then, for every κ ∈ N * , µ is exponential mixing of order κ as in Definition 1.1 for

C γ observables (0 < γ ≤ 2), with θ κ = d -(γ/2) κ+1 /2 .
For endomorphisms of P k (C), the exponential mixing for all orders for the measure of maximal entropy and Hölder observables was established in [START_REF] Dinh | Exponential estimates for plurisubharmonic functions[END_REF]. We recently proved such property for a large class of invariant measures with strictly positive Lyapunov exponents [START_REF] Bianchi | Equilibrium states of endomorphisms of P k II: spectral stability and limit theorems[END_REF]. This was done by constructing a suitable (semi-)norm on functions that turns the so-called Ruelle-Perron-Frobenius operator (suitably normalized) into a contraction. As far as we know, the present paper gives the first instance where the exponential mixing of all orders is established for general dynamical systems with both positive and negative Lyapunov exponents.

The exponential mixing of all orders is one of the strongest properties in dynamics. It was recently shown to imply a number of statistical properties, see for instance [START_REF] Björklund | Central Limit Theorems for group actions which are exponentially mixing of all orders[END_REF][START_REF] Dolgopyat | Multiple Borel Cantelli Lemma in dynamics and MultiLog law for recurrence[END_REF]. As an example, a consequence of Theorem 1.2 is the following result. Take u ∈ L 1 (µ). As µ is ergodic, Birkhoff's ergodic theorem states that

n -1 S n (u) := n -1 u(x) + u • f (x) + • • • + u • f n-1 (x) → µ, u
for µ -a.e. x ∈ X.

We say that u satisfies the Central Limit Theorem (CLT) with variance σ 2 ≥ 0 with respect to µ if n -1/2 (S n (u) -n µ, u ) → N (0, σ 2 ) in law, where N (0, σ 2 ) denotes the (possibly degenerate, for σ = 0) Gaussian distribution with mean 0 and variance σ 2 , i.e., for any interval I ⊂ R we have

lim n→∞ ν S n (u) -n µ, u √ n ∈ I =    1 when I is of the form I = (-δ, δ) if σ 2 = 0, 1 √ 2πσ 2 I e -t 2 /(2σ 2 ) dt if σ 2 > 0.
By [START_REF] Björklund | Central Limit Theorems for group actions which are exponentially mixing of all orders[END_REF], the following is then a consequence of Theorem 1.2. We refer to [START_REF] Bianchi | Equilibrium states of endomorphisms of P k II: spectral stability and limit theorems[END_REF][START_REF] Denker | On the transfer operator for rational functions on the Riemann sphere[END_REF][START_REF] Dinh | Decay of correlations and the central limit theorem for meromorphic maps[END_REF][START_REF] Dupont | Bernoulli coding map and almost sure invariance principle for endomorphisms of P k[END_REF][START_REF] Przytycki | Statistical properties of topological Collet-Eckmann maps[END_REF][START_REF] Szostakiewicz | Stochastics and thermodynamics for equilibrium measures of holomorphic endomorphisms on complex projective spaces[END_REF][START_REF] Szostakiewicz | Fine inducing and equilibrium measures for rational functions of the Riemann sphere[END_REF] for other cases where the CLT for Hölder observables was established in holomorphic dynamics. As is the case for Theorem 1.2, this is the first time that this is done for systems with both positive and negative Lyapunov exponents.

Corollary 1.3. Let f be a complex Hénon map and µ its measure of maximal entropy. Then all Hölder observables u satisfy the Central Limit Theorem with respect to µ with

σ 2 = n∈Z µ, ũ(ũ • f n ) = lim n→∞ 1 n X (ũ + ũ • f + . . . + ũ • f n-1 ) 2 dµ,
where ũ := u -µ, u .

Theorem 1.2 and Corollary 1.3 in particular apply to any real Hénon map of maximal entropy [START_REF] Bedford | Real polynomial diffeomorphisms with maximal entropy: Tangencies[END_REF], i.e., complex Hénon maps with real coefficients and whose measure of maximal entropy is supported by R 2 . They hold also in the larger settings of regular automorphisms of C k in any dimension [START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF], see Definition 2.1 and Remark 3.3, and invertible horizontal-like maps in any dimension [START_REF] Dinh | Dynamics of horizontal-like maps in higher dimension[END_REF][START_REF] Dinh | Geometry of currents, intersection theory and dynamics of horizontal-like maps[END_REF], see Remark 3.4. We postpone the case of automorphisms of compact Kähler manifolds to the forthcoming paper [START_REF] Bianchi | Exponential mixing of all orders and CLT for automorphisms of compact Kähler manifolds[END_REF], see Remark 3.5.

Our method to prove Theorem 1.2 relies on pluripotential theory and on the theory of positive closed currents. The idea is as follows. Using the classical theory of interpolation [START_REF] Triebel | Interpolation theory, Function Spaces, Differential Operators[END_REF], we can reduce the problem to the case γ = 2. For simplicity, assume that g j C 2 ≤ 1 for all j. The measure of maximal entropy µ of a Hénon map f of C 2 of algebraic degree d ≥ 2 is the intersection µ = T + ∧ T -of the two Green currents T + and T -of f [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF]. If we identify C 2 to an affine chart of P 2 in the standard way, these currents are the unique positive closed (1, 1)-currents of mass 1 on P 2 , without mass at infinity, satisfying f * T + = dT + and f * T -= dT -.

Consider the automorphism F of C 4 given by F := (f, f -1 ). Such automorphism also admits Green currents

T + = T + ⊗ T -and T -= T -⊗ T + . These currents satisfy (F n ) * T + = d 2 T + and (F n ) * T -= d 2 T -.
Under mild assumptions on their support, other positive closed (2, 2)-currents S of mass 1 of P 4 satisfy the estimate

(1.1) | d -2n (F n ) * (S) -T -, Φ | ≤ c S,Φ d -n
when Φ is a sufficiently smooth test form. Here, c S,Φ is a constant depending on S and Φ.

We show that proving the exponential mixing for κ + 1 observables g 0 , . . . , g κ with g j C 2 ≤ 1 can be reduced to proving the convergence (we assume that n 1 is even for simplicity)

(1.2) | d -n 1 (F n 1 /2 ) * [∆] -T -, Θ | d -min 0≤j≤κ-1 (n j+1 -n j )/2 ,
where [∆] denotes the current of integration on the diagonal ∆ of

C 2 × C 2 , (z, w) denote the coordinates on C 2 × C 2 and Θ := g 0 (w)g 1 (z)(g 2 • f n 2 -n 1 (z)) . . . (g κ • f nκ-n 1 (z))T + .
A crucial point here is that the estimate should not only be uniform in the g j 's, but also in the n j 's. Note also that the current [∆] is singular and the dependence of the constant c S,Φ in (1.1) from S makes it difficult to employ regularization techniques to deduce the convergence (1.2) from (1.1).

The key point here is to notice that, when dd c Φ ≥ 0 (on a suitable open set), one can also get the following variation of (1.1):

(1.3) d -2n (F n ) * (S) -T -, Φ ≤ c Φ d -n .
With respect to (1.1), only the bound from above is present, but the constant c Φ is now independent of S. This permits to regularize ∆ and work as if this current were smooth. Note also that, although Θ is not smooth, we can handle it using a similar regularization.

Working by induction, we show that it is possible to replace both Θ and -Θ in (1.2) with currents Θ ± satisfying dd c Θ ± ≥ 0. This permits to deduce the estimate (1.2) from two upper bounds given by (1.3) for Θ ± , completing the proof. Definition 2.1. We say that F is a regular automorphism of C k if I ± = ∅ and I + ∩ I -= ∅.

Given F a regular automorphism of C k , it is clear that F -1 is also regular. We denote by d + (F ) and d -(F ) the algebraic degrees of F and F -1 respectively. Observe that d

± (F ) ≥ 2, d + (F ) = d -(F -1 ) and d -(F ) = d + (F -1
). Later, we will drop the letter F and just write d ± instead of d ± (F ) for simplicity. We will recall here some basic properties of F and refer the reader to [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Dinh | Rigidity of Julia sets for Hénon type maps[END_REF][START_REF] Fornaess | Dynamics in several complex variables[END_REF][START_REF] Fornaess | Complex Hénon mappings in C 2 and Fatou-Bieberbach domains[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF] for details. Proposition 2.2. Let F be a regular automorphism of C k as above.

(i) There exists an integer

1 ≤ p ≤ k -1 such that dim I + = k -p -1, dim I -= p -1, and d + (F ) p = d -(F ) k-p . (ii)
The analytic sets I ± are irreducible and we have The set I + (resp. I -) is attracting for F -1 (resp. F ). Let W ± be the basin of attraction of

F (H ∞ \ I + ) = F (I -) = I -and F -1 (H ∞ \ I -) = F -1 (I + ) = I + .
I ± . Set W ± := W ± ∩ C k . Then the sets K + := C k \ W -and K -:= C k \ W +
are the sets of points (in C k ) with bounded orbit for F and F -1 , respectively. We have K + = K + ∪ I + and K -= K -∪ I -where the closures are taken in P k . We also define

K := K + ∩ K -which is a compact subset of C k .
In the terminology of [START_REF] Dinh | Rigidity of Julia sets for Hénon type maps[END_REF], the set K + (resp. K -) is p-rigid (resp. (k -p)-rigid ): it supports a unique positive closed (p, p)-current (resp. (k -p, k -p)-current) of mass 1, that we denote by T + (resp. T -). The currents T ± have no mass on H ∞ and satisfy the invariance relations

F * (T + ) = d p + T + and F * (T -) = d k-p
-T - as currents on C k or P k . We call them the main Green currents of F . They can be obtained as intersections of positive closed (1, 1)-currents with local Hölder continuous potentials in C k . Therefore, the measure T + ∧ T -is well-defined and supported by the compact set K. This is the unique invariant probability measure of maximal entropy [START_REF] De Thélin | Sur les automorphismes réguliers de C k[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF], see also [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . III. Ergodicity, exponents and entropy of the equilibrium measure[END_REF][START_REF] Dujardin | Hénon-like mappings in C 2[END_REF] for the case of dimension k = 2.

Using the above description of the dynamics of F , we can fix neighbourhoods

U 1 , U 2 of K + and V 1 , V 2 of K -such that F -1 (U i ) U i , U 1 U 2 P k \ I -, F (V i ) V i , V 1 V 2 P k \ I + , and U 2 ∩ V 2 C k .
Let Ω be a real (p + 1, p + 1)-current with compact support in U 1 . Assume that there exists a positive closed (p + 1, p + 1)-current Ω with compact support in

U 1 such that |Ω| ≤ Ω . Define the norm Ω * ,U 1 of Ω as Ω * ,U 1 := inf{ Ω : |Ω| ≤ Ω }.
Observe that when Ω is a d-exact current we can write Ω = Ω -(Ω -Ω), which is the difference of two positive closed current in the same cohomology class in H p+1,p+1 (P k , R). Therefore, the norm • * ,U 1 is equivalent to the norm given by inf Ω ± , where Ω ± are positive closed currents with compact support in U 1 such that Ω = Ω + -Ω -. Note that Ω + and Ω -have the same mass as they belong to the same cohomology class.

The following property was obtained by the second author, see [12, Proposition 2.1].

Proposition 2.4. Let R be a positive closed (k-p, k-p)-current of mass 1 with compact support in V 1 and smooth on C k . Let Φ be a real-valued (p, p)-form of class C 2 with compact support in U 1 ∩ C k . Assume that dd c Φ ≥ 0 on V 2 . Then there exists a constant c > 0 independent of R and Φ such that

d -(k-p)n - (F n ) * (R) -T -, Φ ≤ c d -n -dd c Φ * ,U 1 for all n ≥ 0.
Note that in what follows, since T -is an intersection of positive closed (1, 1)-currents with local continuous potentials [START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF], the intersections R ∧ T -and T + ∧ T -are well-defined and the former depends continuously on R. In particular, the pairing in the next statement is meaningful and depends continuously on R.

Corollary 2.5. Let R be a positive closed (k -p, k -p)-current of mass 1 supported in V 1 .
Let φ be a C 2 function with compact support on C k such that dd c φ ≥ 0 in a neighbourhood of K + ∩ V 2 . Then there exists a constant c > 0 independent of R and φ such that

(2.1) d -(k-p)n - (F n ) * (R) -T -, φT + ≤ c d -n -dd c φ ∧ T + * ,U 1
for all n ≥ 0. Proof. As P k is homogeneous, we will use the group PGL(k + 1, C) of automorphisms of P k and suitable convolutions in order regularize the currents R and φT + and deduce the result from Proposition 2.4. Choose local coordinates centered at the identity id ∈ PGL(k + 1, C) so that a small neighbourhood of id in PGL(k + 1, C) is identified to the unit ball B of C k 2 +2k . Here, a point of coordinates represents an automorphism of P k that we denote by τ . Thus, τ 0 = id.

Consider a smooth non-negative function ρ with compact support on B and of integral 1 with respect to the Lebesgue measure and, for 0 < r ≤ 1, define ρ r ( ) := r -2k 2 -4k ρ(r -1 ), which is supported by {| | ≤ r} . This function allows us to define an approximation of the Dirac mass at 0 ∈ B when r → 0. We define Ψ := φT + and consider the following regularized currents

R r := ρ r ( )(τ ) * (R) and Ψ r := ρ r ( )(τ ) * (Ψ) = ρ r ( )(φ • τ )(τ ) * (T + ),
where the integrals are with respect to the Lebesgue measure on ∈ B.

When r is small enough and goes to 0, the current R r is smooth, positive, closed, with compact support in V 1 , and converges to R. Since the RHS of (2.1) depends continuously on R, we can replace R by R r and assume that R is smooth. When goes to 0, φ • τ converges uniformly to φ and (τ ) * (T + ) converges to T + . Using that R is smooth and T -is a product of (1, 1)-currents with continuous potentials, we deduce that the LHS of (2.1) is equal to

lim r→0 d -(k-p)n - (F n ) * (R) -T -, Ψ r .
Since T + is supported by K + and we have dd c φ ≥ 0 on a neighbourhood of K + ∩ V 2 , we deduce that dd c Ψ ≥ 0 on V 2 . By reducing slightly V 2 , we still have dd c Ψ r ≥ 0 on V 2 for r small enough. We will use the last limit and Proposition 2.4 for Ψ r instead of Φ and U 2 instead of U 1 . Observe that for small enough, since U 1 U 2 , we have (τ ) * (dd c Ψ) * ,U 2 ≤ dd c Ψ * ,U 1 . We deduce that the LHS of (2.1) is smaller than or equal to

lim r→0 c d -n -dd c Ψ r * ,U 2 ≤ c d -n -dd c Ψ * ,U 1 = c d -n -dd c φ ∧ T + * ,U 1 .
This completes the proof of the corollary.

In order to use the above corollary, we will need the following lemmas.

Lemma 2.6. Let κ ≥ 1 be an integer and g 0 , . . . , g κ compactly supported functions on C k with g j C 2 ≤ 1. Then there is a constant c κ > 0 independent of the g j 's such that for all 0 , . . . , κ ≥ 0 we have

dd c (g 0 • F 0 ) . . . (g κ • F κ ) ∧ T + * ,U 1 ≤ c κ .
Proof. Set gj := g j • f j for simplicity. We have As we have dd c g0 . . . gκ = 0 near H ∞ , its intersection with T + can be computed on C k . We deduce from the above inequalities and

d k-p - = d p + that (2.2) dd c (g 0 • F 0 ) . . . (g κ • F κ ) ∧ T + κ j=0 F j ) * (ω FS ∧ T + = κ j=0 (F j ) * (ω FS ) ∧ d -p j + (F j ) * T + = κ j=0 d -(k-p) j - (F j ) * ω FS ∧ T + .
We will use that the (p + 1, p + 1)-current ω FS ∧ T + is positive, closed, of mass 1, and its support is contained in

K + ⊂ U 1 . We have (F j ) * ω FS ∧ T + = (F j ) * ω FS ∧ T + , ω k-p-1 FS = ω FS ∧ T + , (F -j ) * (ω k-p-1 FS ) ,
where the last form is positive closed and smooth outside I -. The last pairing only depends on the cohomology classes of ω FS , T + , and (F -j ) * (ω k-p-1

FS

). Hence, it is equal to the mass of (F -j ) * (ω k-p-1

FS

), which is equal to d (k-p-1) j -, see [START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF]. It follows that each term in the last sum in (2.2) is bounded by 1, which implies that the sum is bounded by κ + 1. The lemma follows.

Lemma 2.7. Let D D be two bounded domains in C k . Let g be a function with compact support in D and such that g C 2 ≤ 1. Then there are a constant A > 0 independent of g and functions g ± with compact supports in D and g ± C 2 ≤ 1 such that g = A(g + -g -),

i∂g + ∧ ∂g + ≤ dd c g + on D and i∂g -∧ ∂g -≤ dd c g -on D.

Proof. Let ρ be a smooth non-negative function, compactly supported on D and equal to 1 in a neighbourhood of D. Observe that ρg = g. We denote by z the coordinates of C k . Since g C 2 ≤ 1 and g has compact support in D, there exists a constant

A 1 > 0 independent of g such that |dd c g| ≤ A 1 dd c ( z 2 ). Set g + := A -1 ρ(g + 2A 1 z 2 ) and g -:= 2A -1 A 1 ρ z 2 for some constant A > 0. It is not difficult to check that we have dd c g ± ≥ A -1 A 1 dd c ( z 2 ) on D, i∂g ± ∧ ∂g ± ∞ = O(A -2 ) on D, g = A(g + -g -), and g ± C 2 = O(A -1
). Taking A large enough gives the lemma.

Exponential mixing of all orders for Hénon maps and further remarks

Throughout this section (except for Remarks 3.2, 3.3, and 3.4), f denotes a Hénon map on C 2 of algebraic degree

d = d + = d -≥ 2. Define F := (f, f -1 ). It is not difficult to check that F is a regular automorphism of C 4 = C 2 × C 2 .
We will use the notations and the results of Section 2 with k = 4 and p = 2. We denote in this section by T ± the Green (1, 1)-currents of f , and reserve the notation T ± for the main Green currents of F . Observe that [START_REF] Federer | Geometric measure theory[END_REF]Section 4.1.8] for the tensor (or cartesian) product of currents. We denote by K ± the sets of points of bounded orbit for f ±1 . The wedge product µ := T + ∧ T -is well defined, and is the measure of maximal entropy of f [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . III. Ergodicity, exponents and entropy of the equilibrium measure[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF]. Its support is contained in the compact set K = K + ∩ K -. We have K + = K + × K -and K -= K -× K + . Note also that the diagonal ∆ of C 2 × C 2 satisfies ∆ ∩ I + = ∅ and ∆ ∩ I -= ∅ in P 4 , see also [START_REF] Dinh | Decay of correlations for Hénon maps[END_REF].

T + = T + ⊗ T - and T -= T -⊗ T + , see
We now prove Theorem 1.2. By a standard interpolation [START_REF] Triebel | Interpolation theory, Function Spaces, Differential Operators[END_REF] (see for instance [12, pp. 262-263] and [START_REF] Dolgopyat | On Decay of Correlations in Anosov Flows[END_REF]Corollary 1] for similar occurrences) it is enough to prove the statement for γ = 2, i.e., in the case where all the functions g j are of class C 2 . The statement is clear for κ = 0, i.e., for one test function. By induction, we can assume that the statement holds for up to κ test functions and prove it for κ + 1 ≥ 1 test functions, i.e., show that µ, g 0 (g

1 • f n 1 ) . . . (g κ • f nκ ) - κ j=0 µ, g j κ j=0 g j C 2 • d -min 0≤j≤κ-1 (n j+1 -n j )/2 .
Recall that n 0 = 0. The induction assumption implies that we are allowed to modify each g j by adding a constant. Moreover, using the invariance of ν, the desired estimate does not change if we replace n j by n j -1 for 1 ≤ j ≤ κ and g 0 by g 0 • f -1 . Therefore, we can for convenience assume that n 1 is even.

We fix a large bounded domain

B ⊂ C 2 satisfying K ⊂ B, K -∩ B ⊂ f (B), and K + ∩ B ⊂ f -1 (B).
By induction, the inclusions above imply that

(3.1) K ⊂ B, K -∩ B ⊂ f n (B), and K + ∩ B ⊂ f -n (B) for all n ≥ 1.
Because of Lemma 2.7 and the fact that we are only interested in the values of the g j 's on the support of µ, we can assume that all the g j 's are compactly supported in C 2 and satisfy

(3.2) g j C 2 ≤ 1 on C 2 and i∂g j ∧ ∂g j ≤ dd c g j on B.
For simplicity, write

h := g 1 (g 2 • f n 2 -n 1 ) . . . (g κ • f nκ-n 1 ). We need to prove that | µ, g 0 (h • f n 1 ) -µ, g 0 • µ, h | d -min 0≤j≤κ-1 (n j+1 -n j )/2
since this estimate, together with the induction assumption applied to µ, h , would imply the desired statement. In order to obtain the result, we will prove separately the two estimates

(3.3) µ, g 0 (h • f n 1 ) -µ, g 0 • µ, h d -min 0≤j≤κ-1 (n j+1 -n j )/2 and (3.4) -µ, g 0 (h • f n 1 ) + µ, g 0 • µ, h d -min 0≤j≤κ-1 (n j+1 -n j )/2 .
Set M := 10κ and fix a smooth function χ with compact support in C 2 and equal to 1 in a neighbourhood of B. Consider the following four functions, which will later allow us to produce some p.s.h. test functions:

g + 0 := χ • (g 0 + M ) and h + := χ • (g 1 + M )(g 2 • f n 2 -n 1 + M ) . . . (g κ • f nκ-n 1 + M ) and g - 0 := χ•(M -g 0 ) and h -:= χ• (g 1 +M )(g 2 •f n 2 -n 1 +M ) . . . (g κ •f nκ-n 1 +M )-2(M +1) κ .
Recall that n 0 = 0. To prove (3.3) and (3.4), it is enough to show that (3.5) µ, g + 0 (h

+ • f n 1 ) -µ, g + 0 • µ, h + d -n 1 /2 and (3.6) µ, g - 0 (h -• f n 1 ) -µ, g - 0 • µ, h - d -n 1 /2 .
Indeed, we observe that χ does not play any role in (3.5) and (3.6). Hence, the difference between the LHS of (3.5) and the one of (3.3) (resp. of (3.6) and of (3.4)) is a finite combination of expressions involving no more than κ functions among g 0 , . . . , g κ , that we can estimate using the induction hypothesis on the mixing of order up to κ -1. It remains to prove the two inequalities (3.5) and (3.6).

Denote by (z, w) the coordinates on C 4 = C 2 × C 2 and define φ ± (z, w) := g ± 0 (w) h ± (z). We have the following lemma for a fixed domain U 1 as in Section 2.

Lemma 3.1. The functions φ ± satisfy (i) dd c φ ± ∧ T + ≥ 0 on B × B;

(ii) dd c φ ± ∧ T + * ,U 1 ≤ c κ , where c κ is a positive constant depending on κ, but not on the g j 's and the n j 's.

Proof. (i) For simplicity, we set 0 = 1 := 0 and j := n j -n 1 . Define also gj := g j • f j . In what follows g0 depends on w and gj depends on z when j ≥ 1. Observe that by the invariance property of K + ∩ B in (3.1) and the constraints in (3.2), the following inequalities hold in a neighbourhood of K + ∩ B:

(3.7) i∂g j ∧ ∂g j = (f j ) * (i∂g j ∧ ∂g j ) ≤ (f j ) * (dd c g j ) = dd c gj .
In particular, we have dd c gj ≥ 0 in a neighbourhood of K + ∩ B. Note that for g0 = g 0 the properties hold on B, which contains K -∩ B. Now, since T + is closed, positive, and supported by K + = K + × K -, in order to prove the first assertion it is enough to show that dd c φ ± ≥ 0 on a neighbourhood of (K + ∩ B) × (K -∩ B) in C 4 where χ = 1. In what follows, we only work on such a neighbourhood. We have

dd c φ + = κ j=0 dd c gj l =j (g l + M ) + 0≤j =l≤κ i∂g j ∧ ∂g l m =j,l (g m + M ),
where we recall that g0 is g0 (w) and the other gj 's are gj (z) for 1 ≤ j ≤ κ. For the first term in the RHS of the last expression, we have

(g m + M ) ≤ (κ + 1)(M + 1) κ-1 κ j=0 i∂g j ∧ ∂g j ≤ (κ + 1)(M + 1) κ-1 κ j=0 dd c gj . It follows that dd c φ + ≥ (M -1) κ κ j=0 dd c gj -(κ + 1)(M + 1) κ-1 κ j=0 dd c gj = (M -1) κ 1 - (κ + 1) M + 1 1 + 2 M -1 κ κ j=0 dd c gj , which gives dd c φ + ≥ 0 since the choice M = 10κ implies that (1 + 2 M -1 κ < (1 + 1 κ κ < 3.
Similarly, in the same way, we also have

dd c φ -≥ (M -1) κ κ j=0 dd c gj -(M + 1) κ-1 0≤j =l≤κ i∂g j ∧ ∂g l ,
which gives dd c φ -≥ 0. This concludes the proof of the first assertion of the lemma.

(ii) The second assertion of the lemma is a consequence of Lemma 2.6.

End of the proof of Theorem 1.2. Recall that it remains to prove (3.5) and (3.6). Since µ is invariant and n 1 is even, we have

µ, g ± 0 • (h ± • f n 1 ) = µ, (g ± 0 • f -n 1 /2 )(h ± • f n 1 /2
) . We now transform the last integral on C 2 to an integral on C 4 in order to use the dynamical system F = (f, f -1 ) on C 4 introduced above. We are using the coordinates (z, w) on C 4 = C 2 × C 2 . We will also use the diagonal of C 2 × C 2 , which is given by ∆ = {(z, w) : z = w}.

Recall that we have µ = T + ∧ T -and that the currents T ± have local continuous potentials in C 2 . It follows that the intersections of T ± with positive closed currents on C 4 are meaningful. Moreover, the invariance of T ± implies that (F n 1 /2 ) * (T + ) = d -n 1 T + on C 4 . Thanks to the above identities, we have

(3.8) µ, g ± 0 • (h ± • f n 1 ) = T + ∧ T -, (g ± 0 • f -n 1 /2 )(h ± • f n 1 /2 ) = (T + ⊗ T -) ∧ [∆], (g ± 0 • f -n 1 /2 (w))(h ± • f n 1 /2 (z)) = T + ∧ [∆], (F n 1 /2 ) * (φ ± ) = d -n 1 T + ∧ (F n 1 /2 ) * [∆], φ ± = d -n 1 (F n 1 /2 ) * [∆], φ ± T + .
We apply Corollary 2.5 with the functions φ ± = g ± 0 (w) • h ± (z) instead of φ and the current [∆] instead of R. For this purpose, since ∆ ∩ I + = ∅, we can choose a suitable open set V 1 containing ∆. We also fix an open set V 2 as in Section 2. Since B is large enough, Lemma 3.1 implies that dd c φ ± ≥ 0 on a neighbourhood of K + ∩ V 2 . Thus, we obtain from Corollary 2.5 that d -n 1 (F n 1 /2 ) * [∆] -T -, φ ± T + d -n 1 /2 or equivalently (3.9)

d -n 1 (F n 1 /2 ) * [∆], φ ± T + -T -, φ ± T + d -n 1 /2 .
Together, (3.8), (3.9), and the fact that T -, φ ± T + = T + ∧ T -, φ ± = µ ⊗ µ, g ± 0 (w) • h ± (z) = µ, g ± 0 • µ, h ± give the desired estimates (3.5) and (3.6). The proof of Theorem 1.2 is complete. Remark 3.2. Friedland and Milnor [START_REF] Friedland | Dynamical properties of plane polynomial automorphisms[END_REF] proved that any polynomial automorphism of C 2 is either conjugate to an elementary automorphism which preserves a fibration by parallel complex lines or to a Hénon map as above. So, our results apply to all automorphisms of C 2 which are not conjugated to an elementary one. Remark 3.3. When f is a regular automorphism of C k with k even and p = k/2, the map (f, f -1 ) is regular on C 2k . The same proof as above gives us the exponential mixing of all orders and the CLT for f . The results still hold for every regular automorphism but the proof requires some extra technical arguments that we choose to do not present here for simplicity, see for instance de Thélin and Vigny [START_REF] De Thélin | Entropy of meromorphic maps and dynamics of birational maps[END_REF][START_REF] Vigny | Exponential decay of correlations for generic birational maps of P k[END_REF].

Remark 3.4. When f is a horizontal-like map such that the main dynamical degree is larger than the other dynamical degrees, the same strategy gives the exponential mixing of all orders and the CLT, see the papers [START_REF] Dinh | Dynamics of horizontal-like maps in higher dimension[END_REF][START_REF] Dinh | Geometry of currents, intersection theory and dynamics of horizontal-like maps[END_REF] by Nguyen, Sibony and the second author, and in particular [START_REF] Dinh | Dynamics of horizontal-like maps in higher dimension[END_REF] for the necessary estimates. In particular, these results hold for all Hénon-like maps in dimension 2, see also Dujardin [START_REF] Dujardin | Hénon-like mappings in C 2[END_REF].

Remark 3.5. In the companion paper [START_REF] Bianchi | Exponential mixing of all orders and CLT for automorphisms of compact Kähler manifolds[END_REF], we explain how to adapt our strategy to get the exponential mixing of all orders and the Central Limit Theorem for automorphisms of compact Kähler manifolds with simple action on cohomology. As the proof in that case requires the theory of super-potentials, which is not needed for Hénon maps, we choose to do not present it here.

(

  iii) For every n ≥ 1, both F n and F -n are regular automorphisms of C k , of algebraic degrees d + (F ) n and d -(F ) n , and indeterminacy sets I + and I -, respectively. Example 2.3. (Generalized) Hénon maps on C 2 correspond to the case k = 2 in Definition 2.1.In this case, we have p = k -p = 1 and d + = d -= d, the algebraic degree of the map, see[START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Friedland | Dynamical properties of plane polynomial automorphisms[END_REF][START_REF] Sibony | Dynamique des applications rationnelles de P k[END_REF].

  Since g j C 2 ≤ 1 we have |g j | ≤ 1. Denote by ω FS the Fubini-Study form on P k . Then κ j=0 dd c gj l =j gl κ j=0 (F j ) * ω FS and an application of Cauchy-Schwarz inequality gives ) * (i∂g j ∧ ∂g j ) κ j=0 (F j ) * (ω FS ).

  + M ) ≥ (M -1) κ κ j=0 dd c gj . For the second term, an application of Cauchy-Schwarz inequality and (3.7) give 0≤j =l≤κ i∂g j ∧ ∂g l m =j,l

Acknowledgments. We would like to thank the National University of Singapore, the Institut de Mathématiques de Jussieu-Paris Rive Gauche, and Xiaonan Ma for the warm welcome and the excellent work conditions. We also thank Romain Dujardin, Livio Flaminio, and Giulio Tiozzo for very useful remarks and discussions.

This project has received funding from the French government through the Programme Investissement d'Avenir (LabEx CEMPI /ANR-11-LABX-0007-01, ANR QuaSiDy /ANR-21-CE40-0016, ANR PADAWAN /ANR-21-CE40-0012-01) and the NUS and MOE through the grants A-0004285-00-00 and MOE-T2EP20120-0010.

Regular automorphisms of C k and convergence towards Green currents

Let F be a polynomial automorphism of C k . We still denote by F its extension as a birational map of P k . Denote by H ∞ := P k \ C k the hyperplane at infinity and by I + , I -the indeterminacy sets of F and F -1 respectively. They are analytic sets strictly contained in H ∞ . If I + = ∅ or I -= ∅, then both of them are empty and F is given by a linear map and its dynamics is easy to describe. Hence, we assume that I ± = ∅. The following definition is due to Sibony [30].

CNRS, Univ. Lille, UMR 8524 -Laboratoire Paul Painlevé, F-59000 Lille, France Email address: fabrizio.bianchi@univ-lille.fr National University of Singapore, Lower Kent Ridge Road 10, Singapore 119076, Singapore Email address: matdtc@nus.edu.sg


Rodrigo Rill-García

Eva Dokladalova

Petr Dokládal

Jean-François Caron

Romain Mesnil

Pierre Margerit

Malo Charrier

Inline monitoring of 3D concrete printing using computer vision

Keywords:

The detection of anomalies is at the basis of any 3D printing control. In this paper, we propose a methodology for detection of anomalies based on computer vision. This methodology is composed of three modules: (1) image acquisition, (2) interlayer line and layer segmentation and (3) characterization of the local geometry and texture of the layers and detection of anomalies. The image acquisition is performed with a camera fixed to the printing nozzle. The proposed layer segmentation method recognizes and locates the lines separating the printed layers (F-score = 91%). The third module -taking as input the segmentation and the original image -evaluates the geometry of the layers and the texture of the material. The results are used to detect geometry anomalies when the values are outside the expected range. The material texture is classified into four classes of quality (macro-averaged F-score = 94%). We present the results and show the suitability of our methodology for automatic detection and localization of anomalies on images acquired during a printing session.

Introduction

Construction based on additive manufacturing is now heading towards its maturity [START_REF] Barbosa | Reinventing Construction: A Route of Higher Productivity[END_REF], with practical examples of 3D concrete printing (3DCP) at the building scale [START_REF]All3DP Pro, 8 biggest companies building 3D printed houses[END_REF]. However, despite these high-profile applications, there is still work to be done on material development, construction strategies and the process control to make this technology reproducible and certifiable.

In principle, extrusion-based 3D printing is based on an established workflow: (1) design of the part, (2) definition of the print path and (3) extrusion of the material (printing) [START_REF] Hager | 3D printing of buildings and building components as the future of sustainable construction[END_REF]. Nevertheless, many different process parameters influence the 3DCP extrusion e.g. the velocity of the extrusion nozzle, the distance between the nozzle and the surface of deposition, the dosing pump pressure, the extruded concrete's rheological properties, etc. For this reason, various monitoring methods have been proposed at different stages: before printing (determining the feasibility of the process), during printing (process feedback), and at post-production stages (post-printing processing) [START_REF] Buswell | Inspection methods for 3D concrete printing[END_REF].

The primary objective of inline monitoring (during printing) is to allow the operators to carry out corrective actions in case of anomalous behavior. This principle can be extended to automatic control systems, as seen in the literature for manufacturing processes at different size scales. For example, in [START_REF] Wu | Accurate linear and nonlinear modelbased feedforward deposition control for material extrusion additive manufacturing[END_REF], a feedforward control system is tested for fused filament fabrication. In [START_REF] Shojaei Barjuei | Real-time vision-based control of industrial manipulators for layer-width setting in concrete 3D printing applications[END_REF], a closed-loop control system is tested for 3DCP. In both cases, a vision system is used to measure the parameter of interest: the layer width, which can be controlled by adjusting the velocity of the extrusion nozzle.

Nonetheless, the layer width is not the only parameter that can be measured to be able to take corrective action and minimize defects. Some examples of geometrical deformations caused by process defects, that can be identified by inspecting the layers' thickness, are described in [START_REF] Carneau | Layer pressing in concrete extrusion-based 3D-printing: Experiments and analysis[END_REF]:

• Excessive velocity appears when the velocity of the extrusion nozzle exceeds the extrusion rate, producing discontinuous layers also called longitudinal tearing (Fig. 1(a))

• Over-pressing occurs when the pressing force resulting from the extruded layer exceeds the strength of the penultimate layer, resulting in staggered layer patterns with a loss of control of layer width and thickness (Fig. 1

(b))

• Flow-out occurs when the yield stress of the material is not sufficient to hold its own weight. This introduces distance between the Fig. 1. Printing process defects and their influence in the local geometry of printed samples [START_REF] Carneau | Layer pressing in concrete extrusion-based 3D-printing: Experiments and analysis[END_REF]. printed piece and the extrusion nozzle, producing a poor material deposition (e.g. coiling) (Fig. 1(c))

The measurement of the layers' thickness can be performed by vision systems too. Furthermore, vision systems are not limited to the monitoring of geometrical properties, allowing the analysis of the observed surfaces in terms of properties such as the texture. In this work, using a lateral view of the printed piece, we propose a methodology for inline 3D concrete printing monitoring using computer vision methods (see Fig. 2). Our contributions can be summarized as:

1. A method for automatic layer segmentation. We use a model based on Deep Learning (DL) to locate the interlayer lines in an image (i.e. the lines separating two printed layers), and we segment independent printed layers by using these lines as guide. 2. A monitoring method based on geometrical characterization of the layers. We use image processing to determine the local geometry, based on the interlayer lines, in terms of: orientation, curvature, layer thickness and distance to the nozzle. These measurements are used to detect and locate anomalies with respect to the expected values in the observed layers. 3. A monitoring method based on textural characterization of the layers. We use a machine learning approach to classify the extruded concrete, region-wise, based on textural properties. Textures are classified either as good or as one of fluid, dry and tearing. These classes, which are dependent on the water content of the mixture, are highly related to the layer's rheological properties.

The rest of this paper is organized as follows. In Section 2, we make a literature review of approaches used for inline 3DCP monitoring. In Section 3, we describe our experimental setup, focusing on the inline image acquisition. In Section 4, we discuss the properties of the images used in our experiments and the method used for interlayer line segmentation. The methods for layer geometry characterization and texture classification are described in Section 5. In Section 6, we show our results on experimental images acquired in a 3DCP laboratory. In Section 7, we discuss these results and close the paper by proposing future improvements.

Related literature

For inline monitoring of 3DCP, a lot of research has been published on rheological properties of the material. For example, an online yield strength measurement on an uni-axial flow [START_REF] Ducoulombier | The ''Slugs-test'' for extrusion-based additive manufacturing: Protocol, analysis and practical limits[END_REF], where the mass of concrete drops is measured, provides information about the yield strength and heterogeneity of the material before the printing. The material properties can be controlled in real time using approaches such as a two-component strategy -with the addition of accelerator at the nozzle outlet [START_REF] Gosselin | Large-scale 3D printing of ultra-high performance concrete -a new processing route for architects and builders[END_REF]. This type of information is useful because insufficient stiffness or strength of fresh concrete can lead to an unacceptable cumulative error or even failure, as observed in [START_REF] Archez | Strategy to shape, on a half-meter scale, a geopolymer composite structure by additive manufacturing[END_REF]. On the other hand, other parameters may not make the printing process completely fail but would lead to the creation of visible defects leading to the rejection of the 3D printed part afterwards [START_REF] Buswell | 3D printing using concrete extrusion: A roadmap for research[END_REF].

As a non-invasive approach to inline monitoring, methods based on specialized optical sensors have shown promising results to detect such defects. Most published research in this direction focuses on geometrical features. For example, a 1D Time-of-Flight distance sensor was mounted on the printing nozzle in [START_REF] Wolfs | A real-time height measurement and feedback system for 3D concrete printing[END_REF]. The goal was to ensure that the printing nozzle is always at an adequate distance with respect to the last printed layer, so that the extruded material does not get unexpected deformations in its trajectory from the nozzle to the surface of deposition.

This method allows compensating the effect of non-planar print beds, and it forces the robot to extrude material at a proper distance even if lower layers are deformed by the weight of recent layers (flowout). However, if such a flow-out deformation exist, the final result will be lower than the reference model, with excessively wide layers. This problem can be solved with approaches such as the one proposed by [START_REF] Lindemann | Development of a shotcrete 3D-printing (SC3Dp) technology for additive manufacturing of reinforced freeform concrete structures[END_REF]. There, the printing method is shotcrete rather than layer extrusion; for printing monitoring, the authors use laser triangulation. Similarly to [START_REF] Wolfs | A real-time height measurement and feedback system for 3D concrete printing[END_REF], the authors measure the distance from the nozzle to the last printed layer. In this work, however, the used printing method allows compensating the height of the printed piece if the bottom layers suffered deformation.

Although this method allows correcting the height of the piece, any anomaly in the width of the layer (over-pressing) is ignored. With respect to this problem, we can find alternatives like the one discussed in [START_REF] Kazemian | Real-time extrusion quality monitoring techniques for construction 3D printing[END_REF] based on computer vision.

Computer vision

3DCP monitoring can be performed with data retrieved from economic devices such as RGB cameras. For example, the work of [START_REF] Kazemian | Computer vision for realtime extrusion quality monitoring and control in robotic construction[END_REF] proposed a monitoring system fully based on computer vision, using a camera fixed to the extrusion nozzle. In that work, the image is first converted to grayscale and blurred with a Gaussian filter; the filtered image is binarized with the Otsu's method to segment the printed layer from the background. Finally, the authors map the width of the segmented layer from pixels to inches (with a top view). The patent number US 8944799 B2 [START_REF] Khoshnevis | Techniques for sensing material flow rate in automated extrusion[END_REF] was generated based on this method that adopts computer vision techniques for contour crafting. Low extrusion rates are expected to create thinner layers, while wider layers are produced by high extrusion rates. Another similar work is presented in [START_REF] Shojaei Barjuei | Real-time vision-based control of industrial manipulators for layer-width setting in concrete 3D printing applications[END_REF], where the borders of the printed layer are detected as edges. There, the width of the layers is associated with the velocity of the extrusion nozzle, rather than the extrusion rate.

Towards more robust vision systems, DL has gained a significant attention in recent years. For printed concrete monitoring, we can see the example of [START_REF] Davtalab | Automated inspection in robotic additive manufacturing using deep learning for layer deformation detection[END_REF]. There, DL-based segmentation is performed to distinguish a contour-crafting-printed piece from background (similarly to [START_REF] Kazemian | Computer vision for realtime extrusion quality monitoring and control in robotic construction[END_REF]). In that work, the camera acquires images from an external, post-printing, lateral view rather than from an inline top view. This allows simultaneously monitoring the piece's height and the interlayer lines. Once the printed piece is isolated from the background in the captured images, a surface smoothing filter is used along with the Canny edge detector to identify the interlayer lines in the printed piece. Then, the Hough transform is used to estimate, region-wise, the orientation of these lines. Finally, a module for defect detection locates lines with unexpected orientations inside the printed piece.

For 3D geometry assessment, [START_REF] Nair | Mathematical morphology-based point cloud analysis techniques for geometry assessment of 3D printed concrete elements, Addit[END_REF] compares the 3D model used for printing with a point cloud obtained (post-printing) by scanning the printed piece on a rotatory base. Using techniques from mathematical morphology, the authors calculate a distribution of cloud-to-cloud distances; this distribution allows quantifying global (all piece) and local (layer-wise) errors with respect to the model used to design the printing path.

As seen in the reviewed literature, optical sensors are useful for non-invasive 3DCP monitoring methods. However, research on these methods is still limited in the literature. Furthermore, research on methods based on computer vision is even more scarce. While geometrical assessment is a popular trend in methods using computer vision, monitoring of the surface has been ignored in 3DCP. In fact, surface analysis for inline monitoring of small-scale extrusion is an active topic of research [START_REF] Oleff | Process monitoring for material extrusion additive manufacturing: a state-of-the-art review[END_REF]. Among the methods based on computer vision, the ones based on texture analysis for surface quality assessment can be extended for the monitoring of extruded concrete.

In this work, we propose to detect anomalies in both aspectsthe geometry and the texture -with a single camera, during the printing process. Our methodology is compared in Table 1 with the different works based on optical sensors that were discussed in this section. With the experiments presented in this work, we will show that: (1) an analysis of the interlayer lines allows measuring multiple geometrical parameters of the printed piece to locate anomalies, making the segmentation of these lines an important task; (2) the texture of the extruded layers allows locating abnormal regions when the mixture presents anomalies in the water content; (3) the analysis of the geometry and the texture can be performed from the input of a single RGB camera.

Experimental setup

For the results presented in this paper, we created a dataset from images recorded during a printing session run in the 3D printing laboratory of École des Ponts ParisTech. The cell is composed of a 6axes industrial robot ABB IRB 6620, 1 using the xHEAD printing head developed by XTreeE 2 as extruder. It relies on a bi-component strategy, similar to the one described in [START_REF] Gosselin | Large-scale 3D printing of ultra-high performance concrete -a new processing route for architects and builders[END_REF], with an external feeding pump for the concrete, and a dosing pump for concrete and accelerator in the printing head. A mixing paddle guarantees the homogeneity of dispersion of accelerator in a mixing chamber. The material is the NAG3 concrete developed by Lafarge and tested in [START_REF] Esnault | Experience in online modification of rheology and strength acquisition of 3D printable mortars[END_REF]. The nominal speed is 80 mm/s, with a nozzle diameter of 20 mm, while the pipe diameter is 30 mm. The typical flow rate is 1 L/min, but may slightly vary during the experiment.

The robot trajectories are constructed from a 3D model and a simple slicing algorithm in Grasshopper [START_REF] Rutten | Grasshopper[END_REF], which outputs a series of planes (targets). These targets are later converted into RobotStudio files through HAL [START_REF] Schwartz | HAL: Extension of a visual programming language to support teaching and research on robotics applied to construction[END_REF], a general programming language for industrial robotics. Constraints of orientation were added in order to guarantee the perpendicularity between the camera frame and the deposited layer; the inverse kinematics is solved with HAL, which also indicates if a position is not accessible and allows either repositioning the printed piece in a more accessible region or changing its geometry in the 3D modeling environment.

The robot controller controls the position and speed of the robot, the flow rate of the dosing pump and the flow rate of accelerator. The printing parameters used during the session were dynamically varied to purposely create geometrical and textural imperfections (see Fig. 3).

Using this setup, the proposed monitoring methodology consists of 3 sequential modules (see blue dashed rectangles in Fig. 2):

1. Image acquisition 2. Interlayer line and layer segmentation 3. Local layer characterization and anomaly detection

Image acquisition and interlayer line segmentation

Several setups for image acquisition are possible. As illustrated in Fig. 4(a), we chose a camera fixed to the nozzle in order to have a fixed point of view with respect to the extrusion zone.

Image acquisition setup

To avoid blur in the acquired images, caused by regions of the printed piece outside the depth of field of the camera, the camera orientation should remain perpendicular with respect to the printed wall. This normally requires a seven axis in the robot that allows the extrusion nozzle to rotate around its longitudinal axis. This problem of orientation also exists when printing with non-circular nozzles, which are common place in 3DCP; although out of the scope of this paper, there are technical solutions to this issue [START_REF] Liu | Rotation nozzle and numerical simulation of mass distribution at corners in 3D cementitious material printing[END_REF]. Regarding image resolution, it should allow correctly capturing texture i.e. identifying the aggregates in the concrete. In the experiments for this article, the largest aggregate is roughly 0.3 mm. The desired camera resolution is therefore >3.3 px/mm. Due to eventual deformations or lack of control of width, small apertures are preferred in our applications: for a smaller aperture, the depth of field increases, allowing obtaining sharp images even if the distance between the observed concrete and the lens differs from the expected value.

Small aperture is usually used with long exposure times, which implies a trade-off on this parameter, because mitigation of motion blur requires short exposure times. Here, we are mostly in the case of linear motion blur. If we want the blur to be limited to 1 pixel, we get a constraint on exposure time. Using experimental data, we get an upper bound of 𝑡 𝑒𝑥𝑝𝑜𝑠𝑢𝑟𝑒 ∼ 1∕3000 s, which holds for a wide range of printed concretes, where typical aggregate size remains submillimetric and robot speed in the 100 mm/s range.

Properties of the acquired images

The distance of the lens was adjusted to contain approximately 4 layers per frame, with an expected layer thickness of 6 mm. To acquire images, a Raspberry Pi Camera High Quality Module featuring a SONY IMX477 CMOS sensor with a macro lens was fixed on the printing nozzle (more camera specifications available at [START_REF] Ltd | Buy a Raspberry Pi High Quality Camera[END_REF]). This allowed having a fixed distance (∼13 cm) between the lens and the extruded material. Additionally, we mounted a led lamp on the nozzle to allow a short, and constant, exposure time; this is illustrated in Fig. 4(a). The image resolution is 1280 × 960 pixels; therefore, our images contain ≈40 px/mm. The focus, exposure time and lens aperture were manually adjusted prior to printing the piece and fixed during extrusion.

To have uniformity in our final dataset, from the acquired images we chose a subset of pictures meeting 2 conditions: (1) no background is visible at the sides of the printed piece, (2) the center of the nozzle is coplanar with the center of the observed wall (as illustrated in Fig. 4(a)). We preserved 628 raw images in total; a schematic diagram of a typical frame with acceptable printing is shown in Fig. 4(b).

Interlayer line and layer segmentation

Once an image is acquired, the next step in our methodology is to segment its interlayer lines. On dry pieces, the detection of these lines can be based on finding dark edges (see Fig. 5(a)). However, fresh concrete exhibits specular reflectance as well as layer merging and superposition. Consequently, interlayer lines can take different aspects: black, bright or even not visible at all (see Fig. 5(b)). Furthermore, given the constraints imposed by inline image acquisition (as described in Section 4.1), the problem becomes particularly hard: loss of focus, unexpected motion blur, camera vibration, polluted air (presence of particles), etc.

To approach this problem, we trained a neural network originally used for crack segmentation: U-VGG19 [START_REF] Rill-García | Pixel-accurate road crack detection in presence of inaccurate annotations[END_REF]. It is a fully convolutional neural network inspired by U-net which uses VGG19 as backbone model; the architecture of this encoder-decoder network is illustrated in Fig. 6. By using this network, we get as output a binary mask with the location of the interlayer lines in the input image (as illustrated in Fig. 7).

Given an input image 𝐼, let us denote by 𝑆 ⊂ Z 2 the interlayer lines segmented by the model. This segmentation 𝑆 can be used to extract the layers delimited vertically between each two lines. Let us denote by 𝑅𝑂𝐼 ⊂ Z 2 this set of layers (see Fig. 7(c)).

Inline 3D concrete printing monitoring

Once the interlayer lines and the independent layers are segmented, the next step in our methodology is to perform a local layer characterization for anomaly detection.

Geometrical characterization

The proposed geometrical characterization consists of local measures of the geometry of the printed layers and interlayer lines. These measures, and the posterior anomaly detection, are performed using the following methods.

Local orientation of interlayer line

As discussed by [START_REF] Davtalab | Automated inspection in robotic additive manufacturing using deep learning for layer deformation detection[END_REF], studying the orientation of the interlayer lines allows detecting deformations in the printed layers. To do this, we propose a method based on image filtering using rotating structuring elements [START_REF] Soille | Morphological operators with discrete line segments[END_REF].

First, we perform a thinning [START_REF]Morphological operations on binary images (r2020b[END_REF] to approximate the interlayer lines to 1-pixel-width lines. Let 𝑆 𝑇 ⊂ Z 2 denote the thinned image, 𝑅(𝑑, 𝑙) a line oriented by 𝑑 with length 𝑙, and 𝐷 the set of possible orientations 𝑑 ∈ [-90, 90). The line length 𝑙 is chosen to be 1/5 of the image width. This hyperparameter presents a trade off between the available line angles [START_REF] Soille | Morphological operators with discrete line segments[END_REF] and the locality of the measurement. The value (256 pixels) was selected to be similar to the expected layer thickness, but it is not critical for the method. Let 𝑝 ∈ Z 2 be a point in an image. Then, 𝐹 (𝑝, 𝑑) = [𝑆 𝑇 * 𝑅(𝑑, 𝑙)](𝑝) is the convolution of 𝑆 𝑇 with a line oriented in 𝑑. Finally, our local orientation map 𝜙 is defined as:

𝜙(𝑝) = { arg max 𝑑∈𝐷 F(𝑝, 𝑑), 𝑝 ∈ 𝑆 𝑇 undefined, otherwise (1)

Local curvature of interlayer line

Let L𝑛 be a piece-wise approximation of the 𝑛th interlayer line in 𝑆 𝑇 , using 𝑚 cubic splines [START_REF] Boor | A Practical Guide To Splines[END_REF]:

L𝑛 (𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝑠 1 𝑛 (𝑡 -𝑡 0 ), 𝑡 0 ≤ 𝑡 < 𝑡 1 ⋮ 𝑠 𝑚 𝑛 (𝑡 -𝑡 𝑚-1 ), 𝑡 𝑚-1 ≤ 𝑡 < 𝑡 𝑚 (2)
The knots are selected to split the interlayer line of interest into 𝑚 fixed-length segments. This length is a hyperparameter, and it implies a trade-off between the accuracy and the smoothness of the spline approximation. The value was chosen to be 1/10 of the image width (128 pixels), but it is not critical for the method. The independent variable 𝑡 corresponds to a parametric representation of the spline 𝑠 𝑖 𝑛 (𝑡) = (𝑥(𝑡), 𝑦(𝑡)).

Since this smooth approximation is twice differentiable, the local curvature 𝜅 of L𝑛 at point 𝑡 can be computed as:

𝜅(𝑡) = 𝑥 ′ 𝑦 ′′ -𝑦 ′ 𝑥 ′′ (𝑥 ′2 + 𝑦 ′2 ) 3 2 (3)
The primes refer to the derivatives with respect to 𝑡. The magnitude of 𝜅(𝑡) is the inverse of the radius of the osculating circle touching the point (𝑥, 𝑦) defined by L𝑛 (𝑡). The sign indicates the concavity.

Local thickness of layer

The local thickness of a layer is determined by the distance between its two interlayer lines. We calculate this thickness using the fast Euclidean Distance Transform (EDT) [START_REF] Maurer | A linear time algorithm for computing exact Euclidean distance transforms of binary images in arbitrary dimensions[END_REF].

Recall, 𝑆 𝑇 is the image containing the 1-pixel-wide representation of the interlayer lines. The function 𝐷 ∶ Z 2 → R + is the result of applying the EDT to 𝑆 𝑇 (see Fig. 8(a)). The local maxima of 𝐷 indicate the centers of the printed layers (see Fig. 8(b-c)). Since we expect the layers to be nearly horizontal in the input picture, we estimate the locations of these local maxima along the vertical direction. Let 𝑘 1 = (0, -1, 1) 𝑇 and 𝑘 2 = (1, -1, 0) 𝑇 be vertical kernels. The convolution of the distance map 𝐷 with these kernels allows detecting the maxima of 𝐷: For the points located in the middle of two interlayer lines (the layer center), and inside ROI (see Fig. 7(c)), the layer thickness 𝑊 is

𝑃 𝑀 = {𝑝 | [sign(𝐷 * 𝑘 1 ) * 𝑘 2 ](𝑝) = 2} (4)
𝑊 (𝑝) = { 2𝐷(𝑝), 𝑝 ∈ (𝑃 𝑀 ∩ 𝑅𝑂𝐼) undefined, otherwise (5)
The factor 2 allows obtaining the thickness (the double of the shortest distance from 𝑆 𝑇 to the centerline).

Relative nozzle height

Similarly to [START_REF] Wolfs | A real-time height measurement and feedback system for 3D concrete printing[END_REF][START_REF] Lindemann | Development of a shotcrete 3D-printing (SC3Dp) technology for additive manufacturing of reinforced freeform concrete structures[END_REF], we measure the relative height of the printing nozzle with respect to the surface of deposition. From the lateral point of view proposed in this paper, this height corresponds to the vertical distance from the nozzle to the line receiving the currently printed layer.

Since the camera is fixed to the robot, the center of the nozzle is at a constant position 𝑐 = (𝑥 𝑐 , 𝑦 𝑐 ) in the image. Let 𝐿 2 , the second interlayer line from top to bottom in the image, be the set of points 𝑝 = (𝑥 𝑝 , 𝑦 𝑝 ) separating the currently printed layer from the layer of deposition (see Fig. 9). Then, the relative nozzle height is:

𝐻(𝑝) = { 𝑦 𝑝 -𝑦 𝑐 , 𝑝 ∈ 𝐿 2 undefined, otherwise (6)
With the vertical axis pointing downwards in an image, this height is positive as long as 𝑝 is below the nozzle.

Geometrical anomalies

To summarize the information retrieved by Eqs. ( 1), (3), ( 5) and (6) (orientation, curvature, thickness and height, respectively), we obtain a histogram per measure. Under good printing conditions, the distributions generated from these histograms should be centered and concentrated inside a range of admissible values.

The principle is illustrated in Fig. 10. There, we show the interlayer segmentation maps from a frame with acceptable printing and one exhibiting signs of over-pressing (see Fig. 1 for reference). Additionally, we show the distributions of layer thickness measured pixel-wise by our method for each frame. The hypothesis is that a defect can be modeled as a deviation (Fig. 10(d)) from the distribution expected for a defectless printing (Fig. 10(c)). Then it becomes straightforward to identify the position of these deviations and localize the corresponding portion of layer in the frame. Since the fixed pose of the camera with respect to the nozzle is known, and the pose of the nozzle at the time of acquiring the analyzed frame can be retrieved from the robot's control software, it is possible to calculate the physical localization of the anomalies in the printed piece.

In the next section, we present the characterization of the texture of the printed layers.

Texture characterization

The rheological properties of the extruded material are directly related to the composition of the mixture -which impacts the structural properties of the printed piece. In this work, we focus on the water content: lack of water can cause cold joints at the interface, but an excess of water makes the layers more prone to deformation or even collapse. The water content is one of the several factors affecting the texture of the extruded layers (see Fig. 11). When the texture exhibits abnormal properties, it becomes a good visual indicator of anomalies in the printing process.

In our method, the texture of the printed piece is analyzed locally, layer by layer, in small windows as shown in Fig. 12. The height of the windows is adapted as to cover the entire thickness of either layer, excluding adjacent layers. The windows have a fixed width of 200 pixels (∼5 mm) and are adjacent each to another. The value of this hyperparameter (chosen to have near-square windows) presents a trade off between the available area to analyze per window and the locality of the analysis. Each window is then independently classified, either as good or one anomalous class.

Pre-processing

Before extracting the windows, we level the image's grayscale intensities to preserve textural information while reducing the effects of lightning and shadows. Given the input image 𝐼, first we obtain a local-mean-intensity image using a Gaussian filter with 𝜎 = 40. By subtracting the filtered image from 𝐼, we obtain a new image with an intensity distribution centered around 0. We further subtract the minimum to avoid negative values. The result of this whole process is illustrated in Fig. 13.

Texture features

Each window is analyzed individually to provide a label using a machine learning approach. Our feature extraction is based on two visual descriptors: gray-level co-occurrence matrices (GLCMs) and local binary patterns (LBP). These descriptors are obtained from the maximum rectangular box contained within concrete pixels in the window. Let each of these boxes be called a 𝑇 𝑏𝑜𝑥 .

Before feature extraction, each 𝑇 𝑏𝑜𝑥 is quantized to 𝑞 discrete values. To avoid that outlier intensity values have an undesired influence, each 𝑇 𝑏𝑜𝑥 is clipped to values in the range 𝑚𝑒𝑎𝑛(𝑇 𝑏𝑜𝑥 ) ± 3.1 ⋅ 𝑠𝑡𝑑(𝑇 𝑏𝑜𝑥 ) before quantization. Let the quantized version of the image, in the range [0, q-1], be called a 𝑇 𝑞 𝑏𝑜𝑥 . GLCMs. Lets define a GLCM [START_REF] Haralick | Textural features for image classification[END_REF] as:

𝐺𝐿𝐶𝑀 𝐼 𝛥𝑥,𝛥𝑦 (𝑖, 𝑗) = 𝑤 ∑ 𝑥=1 ℎ ∑ 𝑦=1 ⎧ ⎪ ⎨ ⎪ ⎩ 1, if I(x, y)=i and I(x + 𝛥x, y + 𝛥y)=j 0, otherwise (7)
Here, ℎ and 𝑤 are the height and the width of an analyzed image 𝐼, respectively; 𝛥𝑥 and 𝛥𝑦 are a horizontal and a vertical offset, respectively. 𝐼(𝑥, 𝑦) returns the intensity value of the pixel at the position 𝑥, 𝑦 in 𝐼.

GLCMs are distributions of co-occurring pixel values at a given scale (the chosen offset). From preliminary experiments, 𝑞 = 12 was chosen for GLCM calculation. With 12 intensity levels in a 𝑇 𝑞 𝑏𝑜𝑥 , per each offset we obtain a square 12 × 12 matrix (see Fig. 14). Per each GLCM, we obtain scalar features in terms of statistical properties. With 𝑃 𝑖,𝑗 the value of the normalized GLCM at the position 𝑖, 𝑗 (i.e. the probability of the intensity pair (𝑖, 𝑗)), the features are obtained based on 5 properties: Homogeneity: a measure of the closeness of the distribution of elements in the GLCM to its diagonal. It is defined as: [START_REF] Davtalab | Automated inspection in robotic additive manufacturing using deep learning for layer deformation detection[END_REF] Correlation: a measure of how correlated are the pixels to their neighbors with the given offset. It is defined as:

𝑙𝑒𝑣𝑒𝑙𝑠-1 ∑ 𝑖,𝑗=0 √ √ √ √ √ 𝑙𝑒𝑣𝑒𝑙𝑠-1 ∑ 𝑖,𝑗=0 𝑃 2 𝑖,𝑗
𝑙𝑒𝑣𝑒𝑙𝑠-1 ∑ 𝑖,𝑗=0 𝑃 𝑖,𝑗 ⋅ (𝑖 -𝜇 𝑖 )(𝑗 -𝜇 𝑗 ) √ 𝜎 2 𝑖 ⋅ 𝜎 2 𝑗 ( 12
)
Due to the printing process, the texture properties are anisotropic. Because of this, we extract GLCMs with offsets in the horizontal and vertical directions. Per each direction, we use offset distances from 1 to 50 with step 1. The total number of features obtained per 𝑇 𝑞 𝑏𝑜𝑥 is 500 (2 directions × 50 distances × 5 statistical properties).

LBP. This descriptor consists of assigning a 𝑛-bit binary number to each pixel depending on its 𝑛 neighbors with radius 𝑅 [START_REF] Ojala | Performance evaluation of texture measures with classification based on Kullback discrimination of distributions[END_REF]. Following the neighbors along the hypothetical circle of radius R, each bit is assigned a 0 if the center pixel's intensity value is greater than the corresponding neighbor's value. Otherwise, the bit is assigned a 1. The transformed image contains, per pixel, the decimal equivalent of the binary number calculated according to their neighbors. We use 8 neighbors, therefore the transformed images contain values in the range [0, 255].

From preliminary experiments, 𝑞 = 64 was chosen for quantization to extract LBP (see Fig. 15). After getting the LBP transform of a 𝑇 𝑞 𝑏𝑜𝑥 , we calculate the histogram of intensities of the transformed image to extract a fixed-size vector per 𝑇 𝑞 𝑏𝑜𝑥 . To reduce the dimensionality of this vector, the histograms are calculated using 8 bins. These histograms are normalized so that their magnitudes are not influenced by the size of the 𝑇 𝑞 𝑏𝑜𝑥 . Similarly to GLCMs, we use different radii to extract information at different scales. These radii are from 1 to 31 with step 10. The total number of features per 𝑇 𝑞 𝑏𝑜𝑥 is 32 (4 radii × 8 bins).

Classification

With the extracted textural features, we classify each 𝑇 𝑞 𝑏𝑜𝑥 into 1 of 4 classes as illustrated in Fig. 16: excessively fluid, good quality, excessively dry, and tearing.

The fluid class is characterized by a smooth surface. Since the fluidness is associated with high water content, this class is also characterized by a high specular reflectance. The good class corresponds to the desired material properties. The amount of present water is adequate, producing the appearance of more visible grains at the surface with higher homogeneity and a lower reflectance. The dry class is characterized by a rougher texture, caused by an increased amount of visible grains with bigger size. Since the amount of water is reduced, the material exhibits very low to no reflectance. However, this class is also determined by the absence of cracks. The tearing class is the one characterized by the appearance of cracks. It is also identified by the lack of reflectance and often by a rougher texture than the dry class.

From the images acquired during the printing process, we extracted and labeled a total of 111 texture windows. Each of these windows is labeled by the agreement of two annotators. The distribution of classes is: 24 fluid, 27 good, 24 dry and 36 tearing. The base dataset is composed by the 532 features calculated per each of these images (GLCM + LBP features), and their corresponding class labels. To avoid that some features dominate the classifier because of their magnitude, we standardize each one of them during training.

The model used to learn from this domain is a small convolutional neural network (CNN) with one convolutional input layer and one fullyconnected output layer (see Fig. 17). Because of the feature extraction implementation, adjacent elements in the feature vector represent the same property but with different offsets (see Fig. 18). Since these features have a high likelihood of being correlated, we use the convolution to capture meaningful information from the neighborhood. The kernel size is 6, with a ReLu activation to obtain non-linear outputs. To reduce the likelihood of overfitting, we use dropout regularization with 30% probability. The output layer has 4 neurons and uses a softmax activation so that ∑ 𝑐∈𝑐𝑙𝑎𝑠𝑠𝑒𝑠 𝑃 (𝑐|features) = 1.

Textural anomalies

Similarly to the approach used for geometrical anomaly detection, we can summarize an image in terms of the histogram of window labels. In this case, the categorical distribution should be centered and concentrated in the good class. High occurrence of other classes indicates some anomaly in the printing process.

Experiment results

Segmentation

Once an image from the printing process is acquired, the next step is to segment the interlayer lines to obtain the binary segmentation maps.

We used 128 manually annotated images for training and 32 for testing. The labeling was done in Krita 4.4.5 using a circular brush with a diameter of 20 pixels (see Fig. 19 for an example of manual segmentation). During training, we provided 256 × 256 image crops with a batch size of 4. We implemented U-VGG19 using Tensorflow 2.1.0. Similarly to [START_REF] Rill-García | Pixel-accurate road crack detection in presence of inaccurate annotations[END_REF], we initialized the weights of the encoder with those of the convolutional layers of VGG19 (pre-trained on ImageNet). We used the same optimizer and loss function i.e. Adam and a weighted combination of binary cross-entropy loss and dice score loss. The model is trained during 200 epochs, with an early stop if the testing loss did not improve for 20 consecutive epochs. The learning rate was decreased by a factor of 2 on testing loss plateau with 5 epochs tolerance; the initial learning rate was 10 -4 .

We evaluate the performance of U-VGG19 on the test images in terms of precision, recall and F-score. We calculate each score per image and we report the average over all the test images. There is a possibility of offsets between prediction and annotation caused by the inaccuracy of the annotators at the moment of locating the center of the interlayer lines. To address this, we use a tolerance margin similarly to works such as [START_REF] Amhaz | Automatic crack detection on twodimensional pavement images: An algorithm based on minimal path selection[END_REF][START_REF] Shi | Automatic road crack detection using random structured forests[END_REF]. The chosen tolerance is 2 pixels, meaning that predicted positive pixels no more than 2 pixels away from a positive pixel in the manual annotation are considered true positives. This represents approximately 0.5 mm in real life. The scores are shown in Table 2. The average precision and the average recall have very similar values (around 90%). Therefore, U-VGG19 exhibits a balanced ability to detect interlayer lines while discriminating false positives. An example of automatic interlayer line segmentation using U-VGG19 is illustrated in Fig. 19. With an average F-score of 90.63%, we use the predictions of U-VGG19 to perform the characterization methods proposed in this paper.

Texture classification

In this multi-class context, we evaluate the proposed method with 5-folds stratified cross-validation i.e. each fold has the same proportion of observations with a given class. To further reduce the likelihood of overfitting, per iteration we perform another stratified split to divide the training folds into a training and a validation split; the size of the validation split is half the size of the test split. Finally, we perform data augmentation on the images from the resulting training split, obtaining 16 additional images per texture window. The transformations are rotations (with angles -4, -3, -2, 2, 3, 4 degrees), zooms (with scale factors 0.97, 0.98, 0.99, 1.01, 1.02, 1.03) and illumination rescales and shifts defined by 𝛼 ⋅ intensity + 𝛽 (with pairs 𝛼∕𝛽: 0.9/0, 1.1/0, 1.0/-50, 1.0/50). Per iteration, the approximate number of samples per split is 1326 training, 11 validation and 22 test. Validation and testing are performed on real images only, while the training split does not contain modified versions of the images in the validation and test splits. The approximate class distribution is 22% fluid, 24% good, 22% dry and 32% tearing. We train the CNN with a binary cross-entropy loss using the Adam optimizer with default parameters in Tensorflow 2.1.0. We train for a maximum of 2000 epochs with batch size 32, with an early stop if the loss in the validation split did not improve for 100 consecutive epochs.

We evaluate the trained models in terms of the macro-averaged Fscore in the test splits; the macro F-score gives equal importance to all the classes. Table 3 shows the macro F-score per fold and per class; the single F-scores per fold and class are shown too.

Among the four classes, fluid exhibited the best scores (100% in all folds). On the other side, the most difficult class was dry, with a 90.2% average over the 5 folds. The good and tearing classes have slightly better averages: 93.7 and 92.2%, respectively. From an analysis of the resulting confusion matrices per fold (see Fig. 20), we discovered that the main source of these errors was classifying good and tearing images as dry. Nonetheless, this kind of confusion is understandable since the dry class is, physically, a transitory state between good printing and tearing. Furthermore, these errors are infrequent, as denoted by the average 94.1% macro-averaged F-score over the 5 folds (average of the last row).

With good models for interlayer line segmentation and texture classification, in the next section we show the results of using these models along with the proposed methods for characterization and anomaly detection.

Layer characterization and anomaly detection

In this section, we show two study cases to characterize, based on the proposed measurements: interlayer line orientation and curvature; nozzle height; and layer thickness and texture. Besides the measurement of orientation, expressed in degrees, all the other geometrical measurements are based on pixels as metric unit. Pixels are converted to millimeters by taking 40 px/mm in our setup (Section 4.1.1). The locally measured values are plotted using color maps (see Fig. 21). Additionally, we show the distributions resulting from these measurements (see Fig. 22). The ranges of admissible values, delimited with green dotted lines in the distributions, are user-defined; deviations outside those ranges are considered anomalies, and they are plotted in blue if below the minimum threshold and in red if above the maximum one.

For orientation, we expect near-to-horizontal layers; the userdefined range is (-10, 10) degrees. For curvature, we expect near-tostraight layers; the range is (-0.05, 0.05) mm -1 . The expected layer thickness is 6 mm; the range is (4.5, 7.5) mm. The nozzle should keep a height close to the expected layer thickness; the range is (5, 7) mm. Regarding texture, any class but good is considered as anomaly; the texture classes are associated to a color code similar to the color maps used for geometrical measurements: fluid/blue, good/green, dry/yellow, tearing/red. Additionally to the class, we show the probability of the prediction according to the classifier.

The first study case, analyzed in Figs. 21 and22, portrays an overall acceptable printing. As observed in Fig. 21(b-e), there are almost no geometrical anomalies (except a few high curvature segments); consequently, the distributions obtained from these plots are contained inside the ranges of admissible values, as seen in Fig. 22(b-e).

The second study case, analyzed in Figs. 23 and24, depicts a scenario with severe anomalies. As observed in Fig. 24(b-e), the distributions of the geometrical measurements are very different from the ones in Fig. 22: in all the distributions, there are high density regions outside the range of admissible values. The most extreme case is present in Fig. 24(e), where the distribution is concentrated outside the admissible range. This behavior is likely to cause a coiling effect on the material deposition; in fact, as depicted in Figs. 23(c) and 24(c), we see high curvatures that can be directly related to this phenomenon.

Furthermore, with respect to the texture analysis, we see that the second layer from top to bottom exhibits many regions with tearing. This is the extreme case of texture anomalies, since it implies the appearance of cracks that can jeopardize the structural safety of the printed piece.

As shown by this study case, as well as the first one with overall acceptable printing, the proposed methodology is able to provide an inline characterization of the process based on visual inspection of the last printed layers. With this characterization, the proposed methodology detects and locates anomalies in the extruded layers.

Conclusions and future work

In the present article, we proposed a methodology for computervision based, inline monitoring of 3D concrete printing. This methodology consists of three sequential modules: (1) image acquisition,

(2) interlayer line and layer segmentation, and (3) local layer characterization and anomaly detection.

With the presented experiments using the aforementioned methodology, we show that:

• The use of Deep Learning allows an inline segmentation of the interlayer lines separating adjacent extruded layers. The proposed model obtains an F-score of 91%. We demonstrate that locating these lines in the analyzed images allows a further segmentation of the independent layers. • An analysis of the segmented interlayer lines allows measuring multiple geometrical properties of the printed piece simultaneously. We demonstrate that the measurement of these properties allows detecting and locating anomalies in the printed piece. • An analysis of the texture of the independent layers allows detecting when the extruded material exhibits anomalies in the water content. We suggest 4 classes for this analysis, and the model proposed to classify the textures obtains a macro-averaged F-score of 94%.

Our experiments show that the geometrical and the textural analyses can be performed simultaneously from the input of a single RGB camera. With this, we extend the catalogue of possible geometrical measurements using computer vision with respect to the current proposals in the literature. Additionally, we are the first to propose an analysis based on texture for 3D concrete printing, inspired by the approaches already used in small-scale, extruded-based additive manufacturing.

In the two presented study cases, we show that our methodology provides the operator with the exact position and nature of the detected anomalies. If the results are recorded, they can serve as a quality report after the printing process is finished. This report can serve as a proof of acceptability of the printed piece. Inversely, it can also serve, when severe anomalies are detected to reject the piece -and even before the end the printing to avoid material and time waste. Finally, the nature and severity of the detected anomalies can be reported to the operator to make corrective adjustments.

A future research could extend this monitoring towards a closedloop control. The nature and the severity of the anomalies could be used as feedback to an automated decision system issuing a corrective action on the printing system. It could stem from a rule-based expert system, vector-to-action dictionaries, supervised machine learning models, etc.

The images and annotations used for the results presented in this paper are available as a public dataset named I3DCP: https://github. com/Sutadasuto/I3DCP. The code to reproduce our results is available at https://github.com/Sutadasuto/3dcp_cv_monitoring.

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .

 2 Fig. 2. Summary of the proposed methodology.

Fig. 3 .

 3 Fig. 3. (a) Experimental sample printed with purposely created defects. (b) Example of acceptable printing without purposely created defects.

Fig. 4 .

 4 Fig. 4. The printing process with inline image acquisition.

Fig. 5 .

 5 Fig. 5. In contrast with dry concrete, the interlayer lines in fresh concrete can appear black, white, or invisible, because of layer superposition and specular reflectance caused by the material's humidity.

Fig. 6 .

 6 Fig.6. Architecture of U-VGG19 as provided in[START_REF] Rill-García | Pixel-accurate road crack detection in presence of inaccurate annotations[END_REF].

Fig. 7 .

 7 Fig. 7. (a) One image from the collected dataset; the field of view was adjusted to contain mostly concrete regions including the last printed layer. (b) Interlayer lines segmented by U-VGG19. (c) Regions of interest in white; regions within a layer that are not contained between two interlayer lines are ignored because the layer is not seen entirely in the vertical axis.

Fig. 8 .

 8 Fig. 8. Method proposed to measure the local thickness of the layers. (a) Distance transform of S 𝑇 . (b) Maxima from the distance transform in the vertical direction. (c) Local layer centers.

Fig. 9 .

 9 Fig. 9. Relative nozzle height measurement.

Fig. 10 .

 10 Fig. 10. Distributions of layer thickness, measured pixel-wise. (c) In a defect-free process, the distribution of any measured parameter is contained in an expected interval. (d) When printing defects occur, a deviation from this interval is observed.

Fig. 11 . 1 ∑ 1 ∑

 1111 Fig. 11. Texture of concrete with different water contents.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Windows used to analyze the textures of Fig. 7(a). The width is fixed (200 px ≈ 5 mm) and the height is adapted to contain the entire thickness of the layer.

Fig. 14 .

 14 Fig. 14. Example of extracted GLCM.

Fig. 15 .

 15 Fig. 15. Example of extracted LBP.

Fig. 16 .

 16 Fig. 16. Examples of extracted windows for the four proposed classes (with improved contrast for visualization).

Fig. 17 .

 17 Fig. 17. CNN for texture classification.

Fig. 18 .

 18 Fig. 18. Feature vector for texture classification. Adjacent elements represent the same property but with different offsets.

Fig. 19 .Fig. 20 .

 1920 Fig. 19. Comparison between manual and predicted segmentation. The color code is: (Green) True positives, (Red) False positives, (Blue) False negatives. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  Rill-García: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Writing -original draft, Writing -review & editing. Eva Dokladalova: Conceptualization, Methodology, Supervision, Writing -review & editing. Petr Dokládal: Conceptualization, Methodology, Supervision, Writing -review & editing. Jean-François Caron: Funding acquisition, Resources, Supervision. Romain Mesnil: Conceptualization, Methodology, Resources, Writing -review & editing. Pierre Margerit: Conceptualization, Methodology, Resources, Software, Validation. Malo Charrier: Conceptualization, Methodology, Validation.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Plots from the first study case. The black pixels in the plots b-e correspond to undefined values in the 2D space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 23 .Fig. 24 .

 2324 Fig. 23. Plots from the second study case. The black pixels in the plots b-e correspond to undefined values in the 2D space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 1

 1 Summary of methods based on optics for 3DCP inline monitoring.

	Work	Sensor	Type of property	Measurements	Approximate cost a
	[8]	1D Time-of-Flight	Geometrical	Layer to nozzle height	5 USD
				(top view)
	[9]	Laser triangulation	Geometrical	Layer to nozzle height	6500 USD
				(top view)
	[10]	RGB camera	Geometrical	Layer width	40 USD
				(top view)
	[6]	Monochromatic	Geometrical	Layer width	349 USD
		camera		(top view)
	[11]	RGB camera	Geometrical	Layer orientation	40 USD
				(lateral view)
	[12]	3D scanner	Geometrical	Printed piece cloud to 3D model cloud distance	472 USD
				(lateral view)
	Ours	RGB camera	Geometrical,	Layer to nozzle height; layer thickness; interlayer line	57 USD
			textural	orientation and curvature; texture classification
				(lateral view)

a Estimated with base on the sensor models and/or specifications reported in the respective cited works.

Table 2

 2 Interlayer line segmentation testing scores.

	F-score (%)	Precision (%)	Recall (%)
	90.63	92.14	89.22

Table 3 F

 3 -scores (%) of texture classification using 5-folds stratified cross-validation.

	Class	Fold 1 Fold 2 Fold 3 Fold 4 Fold 5 Class
							macro-average
	Fluid	100	100	100	100	100	100
	Good	100	100	88.9	90.9	88.9	93.7
	Dry	90.9	100	83.3	100	76.9	90.2
	Tearing	93.3	100	92.3	92.3	83.3	92.2
	Class macro-average 96.1	100	91.1	95.8	87.3

Acknowledgments

This work has been supported by the project DiXite. Initiated in 2018, DiXite (Digital Construction Site) is a project of the I-SITE FUTURE, a French initiative to answer the challenges of sustainable city. This research project was financially supported by I-SITE FUTURE, France (ANR-16-IDEX-0003).

Data availability

Links to Github repositories containing our data and source code, respectively, are provided in the conclusions of the article.