Corentin Xuyang Chang
email: xuyang.chang@ens-paris-saclay.fr

Le Gourriérec

François Hild

Stéphane Roux

Corentin Le Gourriérec

Brightness and contrast corrections for stereocorrelation: Global and instantaneous formulation with spatial regularization

Keywords: Digital Image Correlation (DIC), High-speed camera, Brightness and contrast corrections, Impact test, Laminated glass, Spatial regularization

come

Introduction

Digital Image Correlation (DIC or 2D DIC) is a full-field method to measure displacement fields by registering images at different time instants. Thanks to its non-contacting and non-destructive nature, DIC is becoming increasingly popular for analyzing dynamic tests [START_REF] Helfrick | 3D digital image correlation methods for full-field vibration measurement[END_REF][START_REF] Reu | Comparison of DIC and LDV for practical vibration and modal measurements[END_REF][START_REF] Hild | Measurement of kinematic fields via dic for impact engineering applications[END_REF]. Stereocorrelation (SC), as one of the most popular extensions of DIC, exploits images acquired by multiple (i.e., two or more) cameras [START_REF] Luo | Accurate measurement of three-dimensional deformations in deformable and rigid bodies using computer vision[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF][START_REF] Hild | Digital image correlation[END_REF][START_REF] Dufour | Shape, Displacement and Mechanical Properties from Isogeometric Multiview Stereocorrelation[END_REF][START_REF] Wang | Hybrid multiview correlation for measuring and monitoring thermomechanical fatigue test[END_REF] to analyze a common region of interest (ROI) of 3D surfaces. SC is very attractive for different reasons. On the one hand, it grants access to 3D shape characterization and deformation. On the other hand, SC is very flexible to set up experimentally compared to classical DIC methods since the cameras are no longer restricted to being positioned perpendicular to flat specimen surfaces.

The principle of SC is to exploit gray level conservation between the intrinsic texture of the reference configurations (i.e., the gray level of speckle patterns of the ROI or the natural texture of the specimen itself) and the corrected deformed image sequence by minimizing, for instance, quadratic differences between them. However, experimental setups sometimes violate gray-level conservation (e.g., changes in reflectivity, shading, or specular reflection due to large local rotations of the specimen). Gray level variations often raise issues related to calibration and kinematic field measurement [START_REF] Tong | An evaluation of digital image correlation criteria for strain mapping applications[END_REF] in terms of convergence and trustworthiness.

In the early developments of SC, "local" approaches [START_REF] Sutton | The effect of out-of-plane motion on 2D and 3D digital image correlation measurements[END_REF] were proposed to perform registrations for measuring the 3D shape and its corresponding deformations. The external specimen geometry is described as a cloud of 3D points, while a set of displacement vectors describes their kinematics. To account for brightness and contrast variations, the so-called Zero-mean Normalized Sum of Squared Differences (ZNSSD) and Parametric sum of squared difference (PSSD) criteria, for instance, are popular in local SC approach [START_REF] Pan | Equivalence of digital image correlation criteria for pattern matching[END_REF]. These criteria are independent of any affine transformation of gray levels, and hence, even if brightness and contrast corrections are not made explicit, the registration criterion is equivalent to a gray level quadratic difference including constant brightness and contrast corrections (BCCs) over each local zone of interest (ZOI), a.k.a. "subsets" in local SC. Thus, BCCs are constant per subset, but show no specific regularity across subsets. Moreover, the same discretization is set both for the kinematics and for the BCCs, a drastically limiting property.

A dense multi-view approach of SC has been proposed later on (whereby external surfaces are described in forms of facets [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF][START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF][START_REF] Dubreuil | Mesh-based shape measurements with stereocorrelation[END_REF][START_REF] Pierré | Finite Element Stereo Digital Image Correlation: Framework and Mechanical Regularization[END_REF] or freefom surfaces [START_REF] Dufour | Shape, Displacement and Mechanical Properties from Isogeometric Multiview Stereocorrelation[END_REF][START_REF] Wang | Hybrid multiview correlation for measuring and monitoring thermomechanical fatigue test[END_REF][START_REF] Beaubier | CAD-based calibration of a 3D-DIC system: Principle and application on test and industrial parts[END_REF]), and referred to as "global approach". The major benefits of such approaches are two fold:

• comparisons between numerical simulations and experimental measurements can be performed over the exact same mesh, allowing for detailed inspections and validation;

• reducing the degrees of freedom to describe displacement field without "betraying" its accuracy, the conditioning of the Hessian matrix is significantly improved; hence measurement uncertainty is reduced and allows for trustworthy registration.

In such a framework, BC fields are easily discretized over FE meshes in a similar way as the displacement fields [START_REF] Charbal | Regularised digital-level corrections for infrared image correlation[END_REF]. Although such a formulation endows BCCs with some regularity properties (i.e., continuity of both fields over the ROI), the increase in the number of degrees of freedom (from 3 to 5 unknowns per node) may be detrimental to uncertainty. BCCs enhanced SC to Analyze Impact Test on Laminated Glass As a major obstacle to the broader usage of global SC in industrial applications, the Hessian matrix (of the cost function) may still be ill-conditioned, thereby resulting in low accuracy and possibly lack of robustness when very fine meshes are considered, even more so when BCCs are included. To mitigate this issue, regularization strategies have been proposed. The first family, referred to as "hard regularization", consists in restricting the minimization problem to a reduced sub-space of displacement or BC fields using additional assumptions or prior knowledge. In such cases, a set of well-suited shape functions are chosen, for instance, to promote smoothness [START_REF] Rouwane | Architecturedriven digital image correlation technique (addict) for the measurement of sub-cellular kinematic fields in speckle-free cellular materials[END_REF] or to describe problems with closed-form solutions [START_REF] Roux-Langlois | DIC identification and X-FEM simulation of fatigue crack growth based on the Williams' series[END_REF]). The second family (i.e., "soft regularization") aims to penalize (rather than forbid) nonphysical variations (e.g., Tikhonov regularization [START_REF] Dufour | Shape, Displacement and Mechanical Properties from Isogeometric Multiview Stereocorrelation[END_REF], Laplacian penalty operator [START_REF] Charbal | Regularised digital-level corrections for infrared image correlation[END_REF]). Laplacian regularization was chosen herein for BCCs for its flexibility and generic applicability.

A very severe test case was chosen to offer a proof of concept for the proposed implementation of BCCs in FE-based SC, namely, an impact test on PVB-laminated glass. The large stereo angles between different cameras, specular reflections, multiple cracking of glass, and the occurrence of large outof-plane motions of the impacted zone induce very large gray level variations both in space and time, thus making the test case extremely discriminating.

The novelty of this paper is summarized as follows:

-In this work, a penalty term built on the square modulus of the gradient of the brightness and contrast fields is proposed. Such a form is a step forward compared to Ref. [START_REF] Charbal | Regularised digital-level corrections for infrared image correlation[END_REF], which was based on a different penalty (norm of Laplacian) introduced to mimic steady heat transfer problems, meant to physically describe temperature fields as captured by an IR camera. In the present case, the L2norm of the gradient is better suited to a triangular mesh, and is the most convenient low-pass filter easily tunable with no remeshing. It is shown that the BCCs with the new regularization strategy successfully restore the gray-level conservation between the intrinsic texture and stereo images.

-Global stereocalibration is built on the assumption of gray-level conservation between the intrinsic texture and each stereo image.

For stereocalibration with no BCCs, the assumption is not strictly correct, which slows down the convergence speed. The introduction of BCCs significantly accelerates the convergence of the calibration step. In the presented case, it reduced the computation time by at least 55%.

-It is also shown that FE-based SC enhanced with regularized BCCs significantly increased the robustness of the kinematic measurements. The high-frequency components in the spatial displacement field were significantly dampened thanks to adequate spatial regularization applied to both kinematic and BC fields. It is worth noting that the workload of the numerical implementation of BCCs is very small within current global SC codes.

The paper is organized as follows. First, the formulation of global stereo-calibration and stereocorrelation with BCCs is introduced in Section 2, where different spatial regularization strategies for kinematic measurements and BCCs are discussed. The impact test on laminated glass is described in Section 3. Results and discussions are given in Section 4, before concluding. BCCs enhanced SC to Analyze Impact Test on Laminated Glass

FE-based stereocorrelation

This section first recalls the basic principles and notations of global formulations of stereo-calibration and SC. Then the implementation of BCCs is introduced. Last, different types of spatial regularization strategies are discussed.

Perspective projection

The calibration of the stereosystem is a crucial step prior to SC [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications[END_REF]. The perspective projection transformation, which is used in the pinhole model, is introduced. The 3D shape of the specimen is linked to its 2D projections for the m-th camera via projection transformations

s{x (m) } = [Π (m) ]{X} (1)
where s is the scale factor, {X} = {X, Y, Z, 1} ⊤ defines the homogeneous coordinates of any point of the specimen surface, {x (m) } the homogeneous coordinates corresponding to the 2D projection onto the m-th camera plane, and [Π (m) ] the projection matrix associated with the m-th camera.

The projection matrix [Π (m) ] is usually decomposed as the product of two separate matrices

[Π (m) ] = [K (m) ] [RT (m) ] (2)
where the intrinsic matrix of the m-th camera [

K (m) ] [K (m) ] =       f (m) x s (m) k C (m) x 0 f (m) y C (m) y 0 0 1       (3)
includes the focal lengths (f

(m) x , f (m) y
) and the optical center coordinates (C

(m) x , C (m)
y ), and the skewness s (m) k .

The extrinsic matrix of the m-th camera [RT] accounts for the geometrical relationship (i.e., rotation matrix [R (m) ] and translation vector t (m) ) between the specimen frame and that of the m-th camera

[RT (m) ] = [R (m) , t (m) ] (4)
The rotation matrix [R (m) ] is either described using unitary quaternions q = (q 1 , q 2 , q 3 , q 4 ) or Euler angles. The translation vector is written as

{t} ⊤ = {t X , t Y , t Z }.
The intrinsic matrix for each camera can be pre-calibrated using a calibration target/board [START_REF] Garrido-Jurado | Automatic generation and detection of highly reliable fiducial markers under occlusion[END_REF]. In the following, it will be assumed that intrinsic calibrations were already performed for each camera. Thus, to complete the full calibration for the m-th camera, the six extrinsic parameters {P (m) } = {q 2 , q 3 , q 4 , t X , t Y , t Z } ⊤ have to be determined1 . This identification step will exploit an "integrated" approach as detailed below.

Global formulation of stereo-calibration

In the absence of brightness and contrast variations, the gray levels of the sample surface should appear identically for all cameras. Thus, one can design an "intrinsic pattern" [START_REF] Dufour | Shape, Displacement and Mechanical Properties from Isogeometric Multiview Stereocorrelation[END_REF] f (X) of the speckle pattern for all physical points, as the average over all cameras of the gray level of the projections, The general philosophy of the global approach is to use the 3D geometry of the sample surface, dressed with its pattern f (X), as the reference, and all camera images are compared to the projection of this reference onto the corresponding image plane. This construction holds both for calibration and stereocorrelation.

f (m) (x (m) ) f (X) = 1 N N m f (m) (x (m) ([Π (m) ], X)) (5
The calibration of each camera m consists in the determination of the projection matrix [Π (m) ], whose coefficients are collected in a column vec-

tor {P} ⊤ = {{P (1) } ⊤ , {P (2) } ⊤ , • • • , {P (N ) } ⊤ }.
This is done through the minimization of the following quadratic differences

ρ 2 ({P}) = N m=1 ROI f (m) (x (m) ([Π (m) ({P (m) })], X)) -f (X) 2 dX ( 6
)
where N denotes the number of cameras.

Let us emphasize that the intrinsic texture is initially approximate because so are the projection parameters of the cameras. Thus, throughout the calibration procedure, the intrinsic texture is recomputed from Equation (5) at each iteration.

Parameter identification is an inverse problem. In the context of stereocalibration, since a large number of parameters are required to be calibrated (i.e., 6 per camera): two preliminary steps are needed:

1. A pre-calibration is performed, thereby allowing for an initial decent guess for the extrinsic parameters.

2. A sensitivity analysis is mandatory to calibrate the sought camera model accurately.

Sensitivity fields

The sensitivity field relative to any given camera model parameter is obtained by computing a finite difference after applying an infinitesimal variation to the model parameter, from {P} to {P} + {δP}. This parameter variation induces a "pseudo-displacement" field of the corresponding projected points from x (m) to x (m) + δx (m) , from which one can define a sensitivity field.

S (m) P = δx (m) δP (7)
which is a 2D vector field defined at each physical evaluation point X and for each considered parameter P . The sensitivity matrix [S P ] gathers the full set of 6 × N sensitivities.

Integrated approach for stereocalibration

Let us recall that the initial guess of the extrinsic camera parameters is provided from the pre-calibration step; they are gradually refined through the minimization of the cost function (Equation ( 6)) using an iterative Gauss-Newton scheme.

First, the camera-parameter-based Hessian matrix [H P ] is generated using a set of image gradient and the pre-computed sensitivity fields

[H P ] =          H (1) 0 0 0 0 H (2) 0 0 0 0 • • • 0 0 0 0 H (N)          (8)
where [H P ] denotes the (6 × N ) × (6 × N ) Hessian matrix. The sub-Hessian matrix [H (m) ], relative to camera m, reads

H (m) P Q = ROI (m) ((S (m) P • ∇f (m) )(x (m) )((S (m) Q • ∇f (m) )(x (m) ))dx (m) (9)
where ∇f (m) (x (m) ) represents the image gradient of each camera. Similarly, the camera-parameter-based second member {h P } is assembled from the contributions of each camera {h (m) }

{h P } ⊤ = {{h (1) } ⊤ , {h (2) } ⊤ , • • • , {h (N) } ⊤ } h (m) P = ROI (m) ((S (m) P • ∇f (m) )(x (m) ))( f (x (m) (X)) -f (m) (x (m) ))dx (m) (10)
For each Gauss-Newton iteration, the correction of the extrinsic parameters is the solution of the linear system

[H P ]{δP} = {h P } (11)
The above corrections are repeated until reaching stationarity of the parameters

c v = {δP} (it) {P} it) ≤ 10 -3 (12)
where it denotes the current iteration of Gauss-Newton scheme.

BCCs during stereocalibration

Due to lightning conditions and perspective differences for each camera in actual experimental setups, the assumption of gray level conservation between the intrinsic texture f (X) and each image f (m) (x (m) (X)) is not fulfilled. The brightness and contrast variations often impact the stability of stereocalibration and stereocorrelation algorithms. Implementing brightness and contrast corrections (BCCs) into the calibration step is a way to restore the "desired" gray-level conservation. BCCs generally consist of an affine transformation of the gray levels [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications[END_REF]. Because of the freedom offered by these transformations, it is convenient for numerical implementation to normalize each camera image f (m) with zero average and unitary standard deviation

f (m) (x (m) ) = f (m) (x (m) ) -µ(f (m) ) σ(f (m) ) ( 13
)
where µ(f (m) ) is the mean gray level of the considered ROI, and σ(f (m) )

the corresponding standard deviation. For the m-th camera, the corrected reference image, built from the intrinsic pattern, reads

f (m) (x (m) ) = (1 + c (m) (x (m) )) f (m) (x (m) ) + b (m) (x (m) ) ( 14
)
where f (m) (x (m) ) denotes the corrected reference image, and

(b (m) (x (m) ), c (m) (x (m)
)) define respectively the brightness and contrast correction fields. They are evaluated by minimizing the following cost function

T (b (m) , c (m) )) = (1 + c (m) (x (m) )) f (m) (x (m) ) + b (m) (x (m) ) -f (X) 2 ROI (15)
Even if the corrections (b (m) , c (m) ) only operate in the image plane of camera m, for ease of notations, one can map them onto the 3D sample surface. The previous b (m) (x (m) (X)) fields are thus denoted as b (m) (X) (and similarly for c (m) ). The same convention also holds for f (m) .

For arbitrary variations of b (m) and c (m) , the problem is ill-posed. Wellposedness is recovered if one assumes that both brightness and contrast fields obey some constraints of regularity. The brightness and contrast correction fields are discretized in the following using a set of spatial shape functions, BCCs enhanced SC to Analyze Impact Test on Laminated Glass

Φ(X), in the 3D space b (m) (X) = i b (m) i Φ i (X) c (m) (X) = i c (m) i Φ i (X) (16)
where {b (m) } and {c (m) } collect in column vectors the nodal BCCs for the m-th camera.

Within the present finite-element formulation, the sensitivities of brightness and contrast fields are now defined as

S (m) b (X) = 1 S (m) c (X) = f (m) (X) (17)
The BCCs Hessian matrix [H

BCCs ] and second member {h

(m) BCCs } are expressed as [H (m) BCCs ] =    H bb H bc H cb H cc    and {h (m) BCCs } =      h b h c      (18)
where the sub-Hessian matrices and sub-second member read

(H αβ ) ij = ROI ((S α Φ i )(X))((S β Φ j )(X))dX, for (α, β) = b, c (19)
and

(h α ) i = ROI ((S α Φ i )(X))( f (X) -(1 + c (m) (X)) f (m) (X) -b (m) (X))dX (20)
The brightness and contrast unknowns are determined in a single step solving a linear system [H

(m) BCCs ]      b (m) c (m)      = {h (m) BCCs } (21)
The complete stereocalibration algorithm with BCCs is summarized in Algorithm 1

Algorithm 1 Stereocalibration with BCCs Initialize extrinsic parameters {P} (pre-calibration) Compute normalized reference images f (m) (X) ▷ Eq. ( 13) Compute intrinsic pattern f (X) ▷ Eq. ( 5) Compute camera-parameters based Hessian matrix ▷ Eq. ( 8) while pre-convergence of parameter not reached (Eq. ( 12)) do Update camera-parameter based second member ▷ Eq. ( 10) Compute incremental parameter correction ▷ Eq. ( 11) Update camera model parameters Update interpolated image and intrinsic texture ▷ Eq. ( 5) if pre-convergence reached (Eq. ( 12)) then

Compute BCCs Hessian matrix and second member ▷ Eq. ( 18) Compute brightness and contrast fields for all f (m) (X) ▷ Eq. ( 21) Update reference images with BCCs ▷ Eq. ( 14) end if end while Several remarks are worth being made:

1. In practice, it is preferred to implement a staggered approach for BCCs during stereocalibration. More precisely, these parameter corrections are performed between the normalized images and intrinsic texture (red arrows in Figure 1), until reaching pre-convergence of the camera parameters (i.e., c v ≤ 5 × 10 -2 ). Then BCCs (blue arrows in Figure 1) and parameter identification (green arrows in Figure 1 ) are alternated until reaching final convergence.

∇ f (m) (X) = (1 + c (m) (X))∇ f (m) (X) + ∇c (m) (X)) f (m) (X) + ∇b (m) (X) ≈ (1 + c (m) (X))∇ f (m) (X) (22)
It is very expensive to compute ∇c (m) (X) and ∇b (m) (X) over all integration points, but more importantly because BCCs are slowly varying functions, the last two terms of the gradient are negligible compared to the rapidly varying ∇ f (m) (X). Thus in this paper, a first-order approximation was considered.

Finite element formulation of stereocorrelation (FE-SC)

The sought displacement field U should minimize the quadratic differences between the intrinsic texture and the corrected deformed image for all cameras

U(X) = argmin U m=N 1 ρ (m) (X) 2 ROI = argmin U m=N 1 ROI ( f (m) (X) -f (m) U (x (m) (X))) 2 dX ( 23
)
where

f (m) U (x (m) (X)
) represents the corrected deformed image the m-th camera recorded at time instant t.

The displacement field U is sought as a linear combination of a set of

spatial shape functions Φ i U(X) = i=Ns 1 a i (t)Φ i (X) (24)
where a i (t) stands for the nodal displacement to be determined at each time instant t, and {a(t)} collects all nodal displacements in a column vector form.

The nodal displacement can be obtained through a Gauss-Newton scheme [START_REF] Pierré | Finite Element Stereo Digital Image Correlation: Framework and Mechanical Regularization[END_REF][START_REF] Berny | Displacement Uncertainty Quantifications in T3-Stereocorrelation[END_REF]. Hence, at each iteration, the stereocorrelation Hessian matrix [H U ] and second member {h U } are computed

(H U ) ij = N m=1 ROI ((Φ i • ∇ f (m) )(X))((Φ j • ∇ f (m) )(X))dX (25)
and [START_REF] Nourry | Laminated glass behaviour under perforating impact[END_REF] BCCs enhanced SC to Analyze Impact Test on Laminated Glass where ∇ f (m) (X) denotes the gradient of the 3D intrinsic pattern recorded by the m-th camera.

(h U ) i = N m=1 ROI ((Φ i • ∇ f (m) )(X))( f (m) (X) -f (m) U (x (m) (X)))dX
Let us note that the intrinsic pattern f (m) was transported onto the sample surface defined in the 3D space. Similarly, the displacement field basis was also defined as the 3D deformation of the monitored surface. However, since any 3D surface is defined in its 2D parametric space [START_REF] Dufour | Shape, Displacement and Mechanical Properties from Isogeometric Multiview Stereocorrelation[END_REF], the gradient of f (m) is computed in the camera plane

∇ f (m) (X) = ∇ x f (m) (x) δx δX (27)
At each iteration, the nodal displacement vector {a(t)} is corrected by {δa} computed from the linear system

[H U ]{δa} = {h U } (28)
such that

{a(t)} = {a(t)} + {δa} (29)
The gray level residuals, whose norm is minimized, are used to estimate the registration quality. To allow for a fair comparison, it is usual to normalize them residual by the dynamic range of gray levels for each individual camera. Hence local residuals ρ (m) (X), camera-wise residuals ρ

(m)
instant , and global residuals ρ instant are defined

ρ (m) (X) = f (m) (X) -f (m) U (X) dyn( f (m) (X)) (30) ρ (m) instant = ⟨(ρ (m) (X)) 2 ⟩ X ( 31
)
ρ instant = 1 N N m=1 (ρ (m) instant ) 2 (32)
where ⟨•⟩ X denotes spatial averages for all integration points.

The Gauss-Newton scheme is iterated until reaching stationarity of the global instantaneous residual

δρ it instant ρ it instant ≤ 10 -5 (33)
where it denotes the current iteration

BCCs for stereocorrelation

In a similar spirit to the BCCs used for calibration, when SC is used along a long time-series, BCCs are crucial to improve the accuracy and robustness of kinematic measurements, and these corrections, as for displacements, are likely to change over time. Brightness and contrast field, (b (m) (X), c (m) (X)) at time t, are thus sought to minimize the cost function

η 2 (b (m) , c (m) ) = (1 + c (m) (X)) f (m) (X) + b (m) (X) -f (m) U (X) 2 ROI (34)
Once the brightness and contrast fields are determined, one could either update the reference or the corrected deformed images. In the present work, the latter is chosen to avoid re-computing stereocorrelation Hessians at each iteration. The corrected deformed image for the m-th camera is expressed as

f (m) (U,BCCs) (X) = f (m) U (X) -b (m) (X) 1 + c (m) (X) (35)
One could determine displacements and perform BCCs in one shot [START_REF] Hild | Digital image correlation[END_REF][START_REF] Charbal | Regularised digital-level corrections for infrared image correlation[END_REF].

However, proper weighting factors for each block matrix calls for proper weighting of the different sensitivity fields (of different nature, namely, S U , S b and S c ). Alternatively, a staggered approach determines the displacements and the BCCs enhanced SC to Analyze Impact Test on Laminated Glass BCCs sequentially (in a so-called 'U-BCCs-U' pathway). The latter is adopted in this work for its simplicity of implementation.

Spatial regularization

Spatial regularization for kinematic fields

For the presented case test, the kinematic field is assumed to be axisymmetric [START_REF] Gourriérec | High speed stereo-vision study of laminated glass fragmentation upon impact[END_REF]. Consequently, the orthoradial displacement component is null. The motion of each physical point (in polar coordinates) only depends on its initial radial distance to the impact center, r, instead of its full 3D coordinates X.

U(X) = U r (r)e r + U z (r)e z (36)
However, the initial mesh is not adapted to such axisymmetric formulation.

One may introduce a radial discretization over which U r and U z are decomposed, and assign to each mesh node the corresponding value. This operation involves a rectangular matrix [T ] that relates each nodal displacement component to the axisymmetric discretization. Then, the axisymmetric formulation is performed with a mere change in Hessian and second member

[H (U,axis) ] = [T ] ⊤ [H U ][T ] {h (U,axis) } = [T ] ⊤ {h U } (37)

Spatial regularization for BCCs

As mentioned earlier, the introduction of BCCs leads to a significant increase in unknowns and hence to a deterioration of the conditioning and uncertainties.

One way to circumvent this difficulty is regularization [START_REF] Tikhonov | Solutions of Ill-posed Problems[END_REF]. However, in contrast with displacements, one cannot use an axisymmetry property of BCCs, because the orientation of the surface with respect to light sources induces local changes. Adding a penalty on high spatial frequencies of BCCs is a way to limit the effective number of degrees of freedom. Let us introduce [D] the discrete gradient operator integrated over all elements of the mesh, so that the penalty term is proportional to

φ 2 b = N m {b (m) } ⊤ [D] ⊤ [D]{b (m) } (38)
φ 2 c = N m {c (m) } ⊤ [D] ⊤ [D]{c (m) } (39)
where one recognizes [D] ⊤ [D] as a discrete Laplace operator.

Let us note that in the case when SC is hybrid, say, using an IR camera, in addition to visible light cameras, then BCCs are related to temperature fields [START_REF] Charbal | Integrated Digital Image Correlation considering gray level and blur variations: Application to distortion measurements of IR camera[END_REF][START_REF] Maynadier | One-shot measurement of thermal and kinematic fields: Infra-red image correlation (iric)[END_REF], and hence penalizing the temperature gradient is a crude way to model the associated heat transfer problem. This analogy can be further enriched using a finer model [START_REF] Wang | Hybrid multiview correlation for measuring and monitoring thermomechanical fatigue test[END_REF][START_REF] Charbal | Regularised digital-level corrections for infrared image correlation[END_REF].

The total cost function to minimize for BCCs thus reads

φ 2 total = η 2 ({b (m) }, {c (m) }) + w b φ 2 b + w c φ 2 c ( 40
)
where w b and w c represent different weight factors attached to brightness and contrast corrections, as discussed below.

Boundary regularization

The above formulation does not constrain the BC fields at boundaries (a harmonic field, i.e., whose Laplacian vanishes, can match any value at its boundaries). An easy (and cheap) way to bypass this issue is to assign the value of the average of their first neighbors to boundary nodes. BCCs enhanced SC to Analyze Impact Test on Laminated Glass

Regularization weights

The question of weights to be given to regularization terms is always a delicate problem. Let us however point out that a characteristic length scale is hidden in the weights w b and w c . A dimensional analysis shows that each of these weights is proportional to the square of a length scale. Thus, if one normalizes each cost function η 2 , φ 2 b and φ 2 c in φ 2 total from their value computed for a trial field with a well-defined scale, then one can easily attach any scale to the regularization, which then acts as a low-pass filter. Here, the chosen reference trial field is

v ref (X) = exp (in.X/ξ) ( 41
)
where ξ is the characteristic scale, and n an arbitrary unit vector. To change the weight w for setting a filter cut-off length to ℓ, then

w = ℓ ξ 2 η 2 (v ref ) φ 2 (v ref ) (42)
To summarize, Laplacian regularization dampens out (ideally spurious) fluctuations of the corresponding field below a given wavelength.

Impact test on a PVB-laminated glass

Inside the family of safety glasses, laminated solutions are good candidates when the user safety requires good impact performance to be achieved with limited debris projection. This layered material consists of hard (and brittle) glass plies alternating with softer interlayers made of polymers such as polyvinyl butyral (PVB). This transparent shield in transportation, buildings, and public facilities achieves various properties, namely, thermal and acoustic insulation, UV reflection [START_REF] Martín | Polymeric interlayer materials for laminated glass: A review[END_REF], or color changing can be endowed with different properties (e.g., thermal and/or acoustic insulation) including user safety. Moreover, Glass-polymer adhesion allows for large strain regimes beyond glass fragmentation [START_REF] Del Linz | Delamination properties of laminated glass windows subject to blast loading[END_REF][START_REF] Elzière | Large strain viscoelastic dissipation during interfacial rupture in laminated glass[END_REF],

without fragment detachment [START_REF] Fourton | Adhesion rupture in laminated glass: influence of adhesion on the energy dissipation mechanisms[END_REF]. during which delamination of interlayers from glass plies and their stretching dissipate a large amount of impact energy [START_REF] Fourton | Adhesion rupture in laminated glass: influence of adhesion on the energy dissipation mechanisms[END_REF]. Laminated glass thus exhibits superior pre-fracture performance and good post-fracture behavior, thereby ensuring customer safety by keeping large glass fragments during impact.

In order to predict the structural behavior of the laminated glass composite under impact, it is crucial to observe the full history of the event, especially after glass cracking where combined delamination and stretching of the polymer layer may occur at newly formed crack positions. The laminated glass performance is assessed and certified in the building industry via the requirements of EN 356 standard, which consists of three successive impacts on a 1 m 2 laminated glass panel. To study all the mechanisms at work during impact, Nourry [START_REF] Nourry | Laminated glass behaviour under perforating impact[END_REF] has developed an experimental device -named cannon -that mimics the same loading conditions as the standard, on a smaller 300 × 300 mm 2 layered glass plate as shown in Figure 2. The plate is tightened between two polished steel rings allowing bending of fragmented samples during the first tens of milliseconds of an impact.

A compressed air system launches a blunt impactor ended by a steel halfsphere. Its apex meets the glass surface at its center at an adjusted and controlled velocity. The laminate starts to bend, and radial cracks rapidly initiate under the impactor dividing the previously intact plate into triangular beam sections [START_REF] Seshadri | Mechanical response of cracked laminated plates[END_REF]. The continuously bending fragments finally break at a certain radial distance from the impactor center, thereby creating circular cracks. After this initial crack network development, interlayer stretching and polymer-glass delamination occur.

Monitoring the laminate deformation on the opposite side of impact is performed with two synchronized high-speed Photron ® SA-5 cameras positioned on each side of the cannon axis (Figure 2 and Table 1). Each camera recorded the test at 30,000 Hz with a definition of 640 × 376 pixels. Lighting was provided by two 400-W spotlights placed above the cannon axis (Figure 2) to compensate for the short integration time of the high-speed cameras. The speckle pattern was applied over one-half of the specimen surface, except for a small inner central disk. The reasons for painting such an ad-hoc region of interest are two fold:

• First, the upper-right half of the laminated glass remained unspeckled, thereby allowing for direct observation of the development of the crack network. A white vinyl-based paint was deposited onto the semi-circular section of the observed surface. Black dots were sprayed over the surface to create a random pattern. This thin layer had good adhesion properties with glass, and could be deformed into finite strain regimes before breakage without interfering with the glass mechanics and providing enough contrast for SC purposes. The laminate deformation could be captured at late stages of glass fragmentation, thus providing precious information for the post-fracture performance of such assemblies [START_REF] Gourriérec | Multi-instrumented dynamic experiments on laminated glass[END_REF].

• Second, the inner radius of the speckle pattern was precisely prescribed and located at the center of the specimen so that it could be used as a trustworthy boundary during the stereo-calibration procedure. This inner central disk also helped detect millimetric indentation that occurred during contact between the projectile and the plate, which was used to validate the synchronization of the cameras and other sensors. BCCs enhanced SC to Analyze Impact Test on Laminated Glass

An FE-mesh with a slightly reduced external radius was used to avoid potential obscuration and shadowing related to camera perspectives for large laminate deformation. The motivation for full-field measurements of the glass surface was to monitor crack initiation and growth during the test and by the 3D kinematic field measurement of the deformed fragmented laminate to capture the extension of the PVB interlayer responsible for energy dissipation, thereby allowing for the validation of global models. Table 1 gathers all parameters that characterize the setup, while Table 2 recalls the SC analysis parameters. For the sake of simplicity, in the present work, the index t refers to the frame number in the image sequence. The corresponding acquisition rate is F = 30 kHz so that t/F is the real time. The initial impactor speed was about 7.53 m/s. By plotting stereo-image sequences at different instants of time (Figures 3 and4), several qualitative observations are made:

1. Before t = 50, the projectile has not yet reached the glass surface. The reference image is thus chosen as that acquired when t = 50; 2. At t = 70 (0.67 ms after impact), namely, soon after the projectile has reached the laminated glass, radial cracks in glass are first observed;

3. At t = 90 (1.33 ms after impact), along with the densification of radial cracks, orthoradial cracks are initiated;

4. Between t = 110 and t = 130, due to large rotation of some glass fragments, specular reflections of lighting are observed leading to local brightness saturation (particularly for the first camera (Figure 3); 5. Between t = 110 and t = 150, a densification of orthoradial cracks is observed;

6. Despite very large out-of-plane motions during projectile impact, adhesion to the laminated PVB layer limits detachments of small glass fragments. Yet some small pieces ejected from the sample center are seen at about t = 250 (i.e., 6.6 ms after impact). The first fragment detachment is observed when t = 247. This phenomenon violates gray-level conservation and cannot be faithfully captured;

7. After t = 250, the assumption of axisymmetry is no longer valid due to the progressive localization of PVB tearing. The SC analysis was thus carried out in the range t ∈ [50 : 250] for these last two reasons.

Pre-calibration

The pre-calibration scene is shown in Figure 5. CharUco boards were positioned randomly in the 3D space, thereby allowing for a one-shot calibration of the intrinsic parameters of each camera. The "Open-book" calibration target was placed at the center of the laminated glass where two CharUco boards were positioned at each surface center. As illustrated in Figure 5, the local frames of these two CharUco boards matched well the edge of the boards, thus validating the accuracy of the calibration of intrinsic camera parameters.

Since a geometrical coincidence between the calibration target and the laminated glass centers was prescribed, by identifying the extrinsic parameter of these two CharUco boards one easily obtains the initial guess for the extrinsic parameters for the laminated glass through a mere translation.

Calibration

At convergence of the calibration step, the projected mesh is superimposed with each reference image (Figure 6). The edge of the FE mesh matches well the speckled region, thereby validating the calibration accuracy. The decrease in global gray level residuals during stereocalibration when BCCs are included is illustrated in Figure 7. When BCCs are introduced, at iteration #10, the residual suddenly drops from about 4% to less than 3%.

Several comments are worth being emphasized:

1. The initial guess of extrinsic parameters is important since it allows a fair estimate of the intrinsic texture. In the case when the initial guess is too poor, convergence of the proposed algorithm becomes much slower, and a multiscale calibration approach [START_REF] Turpin | Coupling tomographic and thermographic measurements for in-situ thermo-mechanical tests[END_REF] becomes necessary.

2. Similarly, it is important for stereocalibration to reach pre-convergence of the sought parameters before launching the BCCs step. Otherwise, the stability of the proposed staggered approach could be affected.

tion (ending 25 iterations earlier than the standard calibration with no BCCs). with BCCs follows a narrow Gaussian distribution (Figure 9). The restoration of gray-level conservation between the intrinsic texture and corrected reference images proves both the robustness of BCCs and stereocalibration.

In Figure 8(d), the localized gray level residuals in stripes (from bottom left to top right) remain. They may be related to shadowing of the speckle pattern and related perspectives for different stereo cameras.

Kinematic measurements

The measured out-of-plane displacements over the deformed mesh is laid over

Gray-level residual

The normalized local gray-level residual is plotted for several snapshots (at early Figure 12 Let us note that the density of radial cracks from the residual field can only be measured thanks to the BCCs. The identification of crack patterns thus becomes much easier and more trustworthy.

Performance Comparison

Two different global SC variants were considered herein, namely, i) standard SC, ii) SC with BCCs. To allow for a fair comparison between these procedures, the exact same FE mesh and convergence criterion were used.

First, the global gray level residual ρ instant using the two variants are compared in Figure 14. The standard SC algorithm yields levels within the range 4.7-7.5%. The SC algorithm with BCCs leads to lower residuals (varying between 0.5% and 5.3%). In the latter case, the temporal history of the resid- The nodal displacements along the radial and out-of-plane directions with a frame step of 20 is plotted in Figure 15 using these two algorithms. The displacement estimates are actually very similar in terms of amplitude. However, when BCCs were enabled, the spatial fluctuations in the kinematic fields are significantly damped, thereby suggesting that they were mostly spurious due to changes in gray levels. In the present case, It is concluded that BCCs significantly improved the kinematic measurements. In terms of algorithm performance2 , Figure 16 shows that even though more unknowns had to be determined with the new implementation, it converged faster than without BCCs. As consequence, the total CPU tile decreased.

Conclusions

This work presented an FE-based stereocorrelation approach enhanced with brightness and contrast corrections in the stereocalibration and kinematic steps. The analysis of an impact test on laminated glass was used to validate the proposed implementation. The large out-of-plane motions, large stereo-angles of the cameras, and significant brightness and contrast variations make this experiment extremely challenging [START_REF] Chang | Proper generalized decomposition stereocorrelation to measure kinematic fields for high speed impact on laminated glass[END_REF][START_REF] Gourriérec | High speed stereo-vision study of laminated glass fragmentation upon impact[END_REF]. It was shown that including BCCs in stereocorrelation was extremely beneficial in the present configuration:

1. The correction of gray level changes revealed much more clearly cracks, now visible in the fields. For the analyzed test, this increased performance is thus very rewarding. Following references [START_REF] Chang | Proper generalized decomposition stereocorrelation to measure kinematic fields for high speed impact on laminated glass[END_REF][START_REF] Gourriérec | High speed stereo-vision study of laminated glass fragmentation upon impact[END_REF], the kinematic field was regularized using an axisymmetric assumption. The reduced spatial basis successfully decreased the condition number of the Hessian matrix. It was shown that the limited number of spatial degrees of freedom did not degrade the accuracy of the sought kinematics. A "soft" regularization penalizing short wavelength variations of the corrections was introduced for the brightness and contrast fields with a tunable cut-off frequency. With the proposed spatial regularization strategy applied to kinematic and BC fields, the specular reflections were successfully mitigated, and, more importantly, the presence of cracks, as revealed from the residuals, remained untouched; actually enhanced.

A PGD formulation of stereocorrelation was proposed [START_REF] Chang | Proper generalized decomposition stereocorrelation to measure kinematic fields for high speed impact on laminated glass[END_REF][START_REF] Passieux | A Space-Time PGD-DIC Algorithm[END_REF]. Thus, including BCCs as extra degrees of freedom in the PGD framework is appealing to boost the computational speed [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF]. The implementation of BCCs with a PGD formulation is a natural outlook of the work.

with Xavier Brajer, René Gy and Keyvan Piroird. CLG was supported by the ANRT and Saint-Gobain Research Paris through a CIFRE grant.

)

  BCCs enhanced SC to Analyze Impact Test on Laminated Glass

Fig. 1

 1 Fig. 1 Schematic representation of stereocalibration. Red arrows: calibration until reaching pre-convergence of parameters; blue arrows: BCCs ; green arrows: calibration between f and BC corrected normalized image f (m) until reaching convergence of parameters

Fig. 2

 2 Fig. 2 Stereovision set-up [27]. 1-Left-side camera. 2-Right-side camera. 3-400 W spotlights. 4-Laminated glass or, presently shown, the calibration dihedral through a 280 mm diameter circular window.

Fig. 3 Fig. 4

 34 Fig. 3 Image sequence monitored by the first camera from t = 50 to 250

Fig. 5

 5 Fig. 5 CharUco local reference frame plotted in RGB triplet over different stereo cameras: (a) camera 1, (b) camera 2

Fig. 6

 6 Fig. 6 Reference images superimposed with the projected mesh at convergence of the calibration procedure. (a) camera 1, (b) camera 2

Fig. 7

 7 Fig. 7 Global gray level residual as a function of iteration number. Calibration without (resp. with) BCCs is plotted in blue (resp. red)

Fig. 8

 8 Fig. 8 Result of stereocalibration. (a) Normalized reference image for the first camera at all evaluation points, f (1) (X) (b) Intrinsic texture, f (X) (c) Residual field with no BCCs, f (1) (X) -f (X), and (d) with BCCs, f (1) (X) -f (X).

5 .

 5 As a more discriminative validation, the local residual field with or without BCCs are compared in Figure8(c,d) and Figure9. The one with no BCCs reveals a significant background (i.e., it appears in light blue in the bottom and light red in the top parts). Conversely, the local residual field

Fig. 9

 9 Fig. 9 Histograms of residuals at convergence of stereocalibration for cameras 1 (a) and 2 (b) with and without BCCs.

6 .

 6 To highlight the impact of spatial regularization in BCCs, the nodal-wise BCCs are compared in Figure10. With the distribution of BC with no spatial regularization fluctuates more. With spatial regularization (with a length equal to three times the element size), the distribution of B vs. C corrections follow an affine relationship, which highlights the robustness of the implemented spatial regularization.

Fig. 10

 10 Fig. 10 Nodal brightness vs. contrast corrections. Camera 1, with no spatial regularization (a) and with spatial regularization (b). Camera 2, with no spatial regularization (c) and with spatial regularization (d).

  the deformed image acquired by the first camera at several time snapshots t = [50, 90, 130, 170, 210, 250]. The edge of the deformed mesh matches precisely the speckled region along the image sequence. It proves that the SC analysis successfully delivered accurate results in spite of the aforementioned experimental challenges.

Fig. 11

 11 Fig. 11 Deformed mesh with color-coded Uz displacement fields laid over different frames acquired by camera 1.

  and later Figure13 stages). In the early stage, the initial graylevel residual after BCCs is very low for the most part of both cameras in the interval t = [50, 90]. With the increase of impact load and densification of orthoradial cracks t = [90, 130], the speckle pattern itself begins to crack, leading to locally higher residuals in the form of scattered spots but the background of residual remains almost zeros. A major radial crack close to the inner radius in the middle part of specimen is first observed at t = 130.

Fig. 12

 12 Fig. 12 Normalized local residual fields for different snapshots t = [50, 70, 90, 110, 130, 150] of cameras 1 and 2.

Fig. 13

 13 Fig. 13 Normalized local residual fields for different snapshots t = [150, 170, 190, 210, 230, 250] of cameras 1 and 2.

  ual is remarkably smooth. The sharp spikes in global residuals for the former variants are due to fleeting specular reflections and local brightness and contrast variations (t ∈ [70, 80] ∪ [130, 150] ∪ [240, 250]) and are mostly canceled by BCCs. Such improvement in terms of convergence and robustness of the SC code leads to more trustworthy results while increasing the number of frames that could be analyzed with the present method.

Fig. 14

 14 Fig. 14 Global residual as a function of image number with standard SC and SC with BCCs.

Fig. 15

 15 Fig. 15 Nodal-wise displacement fields along (a) in-plane and (b) out-of-plane directions with a frame step of 20 with standard SC and SC with BCCs.

Fig. 16

 16 Fig. 16 Performance comparison. (a) Total CPU time and (b) average iteration per image with and with no BCCs.

2 .

 2 The competition between kinematic, brightness, and contrast fields were greatly reduced by including the BCCs between the reference image and the corrected deformed image for stereocorrelation BCCs enhanced SC to Analyze Impact Test on Laminated Glass purposes. The displacement fields were less noisy compared to those obtained with no BCCs.3. After restoring the gray-level conservation between each reference image and the intrinsic texture, the stereocalibration step converged much faster.[START_REF] Luo | Accurate measurement of three-dimensional deformations in deformable and rigid bodies using computer vision[END_REF]. With an existing FE-based SC code, the numerical implementation of BCCs was easy.

Table 1

 1 DIC hardware parameters of the stereoscopic system[START_REF] Chang | Proper generalized decomposition stereocorrelation to measure kinematic fields for high speed impact on laminated glass[END_REF]

	Cameras	Photron ® SA-5
	Definition	640 × 376 pixels
	Color filter	none
	Gray Levels rendering	8 bits
	Lens	Nikon AF Nikkor 24mm f/2.8D
	Aperture	f/11
	Field of view	500 mm × 300 mm = 0.15 m 2
	Image scale	40 mm / 56 px≈ 0.7 mm/px (in the center)
	Stereo-angle	right camera: 33°, left camera: 37°S
	tand-off distance	between 70 and 80 cm
	Image acquisition rate	30,000 fps
	Patterning technique	B/W paints
	Pattern feature size	in the range of 3 to 10 pixels

Table 2

 2 Stereocorrelation analysis parameters

	DIC software	Correli 3.0 [30]
	Image filtering	None
	Element sizes	1.5-1.7 mm (≈ 25 -26 px, see Figure 6)
	Shape functions	Linear (T3 elements)
	Evaluation points (per element)	528 (i.e., N IP ≈ 23)
	Matching criterion	Sum of squared differences (23)
	Interpolant	Linear
	Displacement noise-floor	[σ(Ux), σ(Uy), σ(Uz)] = [8.2, 8.5, 4.7] µm
	4 Results and discussions
	4.1 Qualitative observations

For unitary quaternion q, it is sufficient to determine three components (i.e., [q2, q3, q4]) to characterize rotation matrices[R]

BCCs enhanced SC to Analyze Impact Test on Laminated Glass

BCCs enhanced SC to Analyze Impact Test on Laminated Glass 3. The Gauss-Newton algorithm used in this study rests on the assumption that the cost function minimum is close to 0. For stereocalibration with no BCCs, this assumption is not correct and hence convergence tends to be slow (Figure7). Conversely, BCCs enhanced the parameter-based Hessian matrix, which significantly boosts the convergence speed of stereocalibra-

The computations were run on workstation with a 2x Intel(R) Xeon(R) CPU E5-2630 v4 (Broadwell) 2.2 GHz and 20 cores.

Mechanics 61, 771-790 (2021)

Acknowledgments

The "cannon" experiment was made possible with different pieces of equipment provided by Saint-Gobain Research Compiègne, Saint-Gobain Research Paris and ENS Paris-Saclay, who are gratefully acknowledged. Many thanks to Quentin Guy de Chamisso, Grégorie Remeur and Richard Villey for conducting the cannon tests. We acknowledge useful suggestions and discussions