i Cette dernière année de thèse aura été riche en évènements. Entre la crise sanitaire, la recherche d'un premier poste à la DGA, la rédaction et soutenance du projet PreDICT, la recherche d'un nouveau logement et la candidature au Corps Technique de la DGA, je n'aurais pas réussi à tout mener de front sans l'aide et le soutien que nombre de personnes m'ont apporté. Je souhaite tous les remercier ici, à commencer par ceux que j'oublierai sans doute dans les lignes qui suivent... Je souhaiterais commencer ces remerciements par les doctorants, ingénieurs de recherche et jeunes docteurs que j'ai croisés, à la fois pour les discussions scientifiques et les moments de convivialité partagés.

Ensuite, un grand merci à l'ensemble des membres de mon Jury de thèse d'avoir accepté de participer à ma soutenance. Merci aux rapporteurs, Julien Réthoré et Stéphane Roux, pour leurs relectures attentives, d'autant que le rendu du rapport précédait de peu la conférence PM-iDICs 2021 organisée à Nantes, à Jean-Denis Durou pour son intervention éclairante au laboratoire, dans le cadre des "séminaires image", qui m'a permis de mieux appréhender les thématiques issues de la vision par ordinateur, à Michel Grédiac pour la présidence bienveillante de mon Jury de thèse ainsi que ses retours concernant mon manuscrit, also special thanks to Elizabeth Jones for kindly accepting to participate in my defense in French.

Je souhaite remercier la Fondation Jean-Jacques et Felicia Lopez-Loreta et son Conseil d'avoir primé le projet PreDICT et de donner à ces travaux de thèse la possibilité d'aller plus loin en visant des développements plus appliqués. Un grand merci Nghia au passage pour ton aide dans la structuration de ce projet. Merci aussi à Florent, Renaud et Nicolas pour les discussions et leur aide.

Je tenais également à remercier la DGA pour la qualité de la formation que j'ai reçue depuis l'École d'ingénieur jusqu'au Doctorat, pour la définition de mon premier poste en lien avec mes travaux de thèse et mes appétences, et pour m'avoir permis de réaliser mon rêve d'enfant de devenir pilote.

Merci Damien pour la grande liberté laissée à cette thèse dans ton rôle de référent technique, ainsi que pour ton aide dans la définition de mon premier poste aux Techniques Aéronautiques. Merci Christophe pour les discussions enrichissantes autour de la corrélation, ta curiosité et aussi ton engagement.

Évidemment, un grand merci à mes directeurs de thèse. Merci Robin pour tes relectures attentives de l'ensemble des productions de ces trois dernières années, merci pour tes questions qui ont souvent soulevé les points chauds et permis de mieux comprendre ce que l'on faisait. Merci Jean-Charles pour ta pédagogie et ton aide à comprendre la portée du travail réalisé et à l'expliquer dans les productions. Et last but not least, merci Jean-Noël d'avoir endossé le rôle de directeur de thèse et d'avoir fixé un cap à ces travaux, sans me contraindre à suivre une route fixée, mais en m'orientant avec beaucoup de subtilité. J'ai eu beaucoup de chance de tomber sur cette équipe, très complémentaire (mathématiques, mécanique, expérimental), le tout dans une super ambiance. Également membres de l'équipe, merci aux co-bureaux Ali et Guillaume de m'avoir supporté pendant presque 3 ans. Merci Guillaume pour le temps que tu as pris au début pour répondre aux questions que je me posais, et m'avoir permis de mettre plus vite le pied à l'étrier. Merci Ali pour les randos, ton aide pour beaucoup trop de choses et les discussions.

Un grand merci aux amis, à toute la famille, la [CHK+], aux frangins et à Papa et Maman d'avoir été là, en particulier ces derniers temps, qui ont été indéniablement difficiles et exigeants. Merci Papa et Maman de nous avoir permis de faire des études dans de bonnes conditions à Cri-Cri, Flo et moi. Merci Maman de nous avoir suivis dans les matières littéraires et Papa dans les matières scientifiques, merci de nous avoir donné la chance d'en être là où nous en sommes aujourd'hui.

Enfin merci Morgane de partager ma vie

Notations

General notations

n th -order tensors are underlined with n lines.

Camera model and photogrammetric calibration

For all point M of the three-dimensional space we will denote m the point defined as the intersection of the considered camera retinal plane and the ray going from the camera optical centre O c to M . As m is a projection of M , it will be considered as a point of a two-dimensional space, and thus described thanks to two coordinates. Other notations are inspired to a large extent from previous works [Garcia 2001].

M

Vector of the three-dimensional space, e.g. M w m

Vector of the two-dimensional space, e.g. m pix M w M coordinates in the world coordinate system R w :

M w = O w M = (X w Y w Z w ) T M c
M coordinates in the camera coordinate system R c :

M c = O c M = (X c Y c Z c ) T m r
m coordinates in the retinal coordinate system R r : m r = O r m = (x r y r ) T m pix m coordinates in the picture coordinate system R pix :

m pix = O pix m = (u v) T m pix
Homogeneous coordinates of m in the picture coordinate system R pix :

m pix = O pix m = (

Photometric Digital Image Correlation (PhDIC) formulation

Most quantities in this section may depend on the considered camera (often denoted with a subscript c) and time step (superscript t). Cameras may be indexed with a subscript i to avoid misleading notations (e.g. Z t c,i is the unit Z vector associated to the camera reference frame of camera i at time t).

L e

Emitted radiance L f

Falling radiance Φ r or E Irradiance α Camera throughput Ω f Solid angle of the falling light θ f

Colatitude associated with the direction of the falling light ϕ f Azimuth associated with the direction of the falling light r(X)

Unit vector going from X to the considered camera optical centre n Normal vector to Ω n e (θ e , ϕ e ) Unit vector pointing the direction (θ e , ϕ e ) (light emitted) n f (θ f , ϕ f ) Unit vector pointing the direction (θ f , ϕ f ) (light falling) f (θ f , ϕ f , X, θ e , ϕ e ) Bidirectional Reflectance-Distribution Function at point X and for direc-

Application on a real test case

As camera positions are completely different between the shape measurement step and the displacement measurement one, pictures associated with the first step are denoted I i and deformed state ones J j . This is to avoid misleading the reader with previous notations, as they might assume that I 0 i and I 1 i correspond to the camera i lying at the same place between t = 0 and t = 1. This is not the case here.

For the same reason, quantities such as visibility, camera projection model... are denoted with a (e.g. V j , P j ).

General introduction Context

Taking all sectors into consideration, the global aircraft industry generates about 2-3% of all carbon gas emissions [START_REF] Becken | What role for offsetting aviation greenhouse gas emissions in a deep-cut carbon world[END_REF]. Concurrently, and despite the ongoing sanitary and economic crisis, the long term aircraft industry growth is expected to average 5% a year until 2030 [European Commission 2020]. This trend has been observed over the last decades and, as shown in Figure 1, air traffic has proven itself to be relatively resilient throughout the crises that it has faced. Based on these considerations, the question is not whether air traffic will recover from COVID 19 pandemic but when it will recover. To mitigate the impact of aircraft industry on carbon gas emissions associated with this exponential growth, public authorities defined ambitious road maps. For instance, by 2050, a 75% reduction in CO 2 emissions compared to 2000 levels stands for one of the most challenging goals [European Commission 2015, Part 2 "The Rationale for CleanSky 2"]. To achieve these goals, ground-breaking innovations in aircraft design cycles are needed, including in terms of weight and structural mechanics [ACARE 2017, "Fostering implementation"]. Currently, aircraft design cycles are based on the pyramid of tests [Rouchon 1990]. Such a pyramid is sketched in Figure 2. It relies on the matching of data from tests (left hand side) together with data from models (right hand side). This matching is performed at different scales, from coupon tests (bottom of the pyramid) up to full scale tests (top). It should be noted that as scale increases (from bottom to top), the number of tests performed and models available shrinks, resulting in the pyramidal shape. The reasons for that are simple: as scale increases, costs associated with tests soar and, typically, only one full-scale test will be performed for a Figure 2: Testing pyramid in aeronautics. Increasing test scale goes hand in hand with increasing costs, thus the number of tests decreases as the test scale increases. Regarding modelling, the greater the scale, the wider the gap between predictions and tests. As a result, the number of available models decreases as scale increases. Based on [Rouchon 1990] and [Passieux et al. 2015a, Figure 1]. Pictures are courtesy of DGA Aeronautical Systems.

given aircraft. At the same time, the number of available models decreases at an even faster rate, because of the lack of predictability at the largest scales.

Most, if not all, mechanical models rely on a set of parameters describing the behaviour of a material under given conditions (geometry, loads, boundary conditions...) and at a given scale. The knowledge of these parameters and of the associated uncertainties allows one to predict the response of a system within a margin of error. When this margin becomes too wide, tests need to be performed in order to adjust the model. That is to say identify more precisely the model parameters. However, uncertainty on the identified parameters is closely related to the amount of available relevant data -without any data, there is no chance to retrieve the associated parameters. As scale increases the number of parameters to identify increases as well, because geometries, loadings, anisotropy and/or inhomogeneities, and boundary conditions become more complex to represent. For instance, the geometry, or shape, of a part in an assembly may substantially differ from its geometry when considered alone, because of the mechanical loading resulting from the weight of other parts in the assembly. Also, this shape in a complex environment may depend on the manufacturing process or on the way in which the system was assembled. For two full-scale tests, this may result in two different initial geometries for a same model. But as already explained, the greater the scale, the smaller the number of tests. And here arises our main concern: the number of parameters to identify is greater at the largest scales, while the amount of potentially available data shrinks because of the scarcer tests. There is also another issue which is linked to the measurement point density. Instrumenting tests accounts indeed for a significant portion of costs, because of the unit cost of sensors together with the time spent by operators to set them up. Hence, the measurement point density may be reduced at the top of the testing pyramid (see Figure 3). This results in an imbalance at the largest scales between the amounts of available data and parameters to identify. This is a major issue in the aircraft design cycles because increasingly large uncertainties tend to propagate in the upscaling process. On top of that, aircraft structures are based on thin, lightweight, slender parts which are much more prone to local non-linear phenomena (see Figure 4). Unfortunately, non-linearities appear to be particularly sensitive to uncertainties and associated structures mechanical responses may vary dramatically [Stanley 1985]. These two aspects, uncertainty propagation together with non-linearities, explain the current lack of predictability at the largest scales provided by the models. Hence, empirical manufacturing knockdown factors are defined to account for the uncertainty on the structures response. This can result in the best-case scenario in oversized designs (which is detrimental to the environment), and in the worst-case scenario in unpredicted mechanical failures. Thus, although the testing pyramid has been the source of countless successful designs, it can be seen, in its current form, as a barrier to innovative and more efficient concepts [START_REF] Certest | Certification for Design: Reshaping the Testing Pyramid[END_REF]]. As explained earlier, ground-breaking innovations are needed to meet the targets defined in terms of CO 2 emissions for the aircraft industry, making the reshaping of the testing pyramid necessary. Indubitably, what comes out of previous considerations is that a better understanding of aircraft structures mechanical behaviour is needed to be able to keep on designing ever more efficient concepts. In order to achieve this goal, two main issues should be tackled. First, tests should not be considered only as a way to validate predictions made by "perfect" models anymore. Instead, a true test-simulation dialogue should be aimed for. A dialogue where models provide mid-scale tests with data such as boundary conditions and mechanical loading stemming from full-scale loads. But most importantly, a dialogue where similar data acquired during tests are, in turn, used to enrich models by accounting for aforementioned uncertainties. Second, a way to reduce costs associated with data acquisition on large and complex geometries should be found so as to increase the measurement point density.

Of course, the first identified issue is closely related to the second one, as extracting enough data from tests is a prerequisite to the test-simulation dialogue. However there is also an organisational reluctance which is beyond the scope of the present manuscript but that we wish to briefly discuss. As explained earlier, people are used to a one-way test-simulation matching, where tests are performed according to the assumption that, for instance, the "perfect" CAD model is a good approximation of the true geometry. Thus, there is a lack of feedback loop and associated processes to include information coming from tests in simulations. Also, relying on sparse information stemming from tests is actually somewhat comfortable for (some) people in charge of explaining discrepancies between tests and simulations. Even though the model used does not account for the observed mechanical behaviour, it is indeed possible, when relying on sparse measurements, to adjust model parameters to fit test data.

Measurement method specifications and choice

In order to increase measurement point density, using full-field measurement methods [Grédiac & Hild 2013] stand for an interesting strategy. Compared to punctual measurements, it is indeed possible to increase the total number of measurement points by up to a factor 100 while decreasing costs associated with measurement device setup and consumable items. From previous considerations, a flexible measurement technique capable of operating at a broad range of scales, on complex geometries while remaining as inexpensive as possible is desirable. Also, the measurement method should be able to facilitate test-simulation dialogue, as it is one of the two main identified issues of the previous section. Of course, the measurement data should provide relevant information with respect to the models, in order to reduce uncertainties listed above.

Regarding mechanical tests, a wide variety of full-field measurement methods exist [START_REF] Grédiac | The use of full-field measurement methods in composite material characterization: interest and limitations[END_REF], Grédiac & Hild 2013]. Some allow to retrieve an in-plane displacement field (Digital Image Correlation (DIC) [START_REF] Horn | [END_REF], Lucas & Kanade 1981[START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF], Localized Spectrum Analysis (LSA) [Grédiac et al. 2019] which is a grid method improvement [Grédiac et al. 2016, Parks 1969], Moiré interferometry [START_REF] Post | High sensitivity moiré: experimental analysis for mechanics and materials[END_REF]), an out-of-plane displacement field (Moiré interferometry [Asundi et al. 1989]), a three-dimensional displacement field over a surface (Stereo Digital Image Correlation (SDIC) [Lucas & Kanade 1981, Sutton et al. 2009]), or even over a whole 3D volume (Digital Volume Correlation (DVC) [START_REF] Bay | [END_REF]). Shapes of specimens can also be extracted (X-ray tomography [START_REF] Desrues | [END_REF], Fringe projection [Zhang 2010], SDIC [Sutton et al. 2009], Deflectometry [Nguyen et al. 2019]). Others provide a strain field as Photoelasticity [START_REF] Ramesh | [END_REF], Germaneau et al. 2008] or (commonly out-of-plane) displacement gradients as Shearography [Hung 1982]. Even temperature can be obtained by the mean of thermal paints [Lempereur et al. 2008] or InfraRed Thermography (IRT) [Gilblas et al. 2014].

In this context, DIC and its three-dimensional 1 extension, SDIC, stand for interesting candidates 2 . These are indeed contact-free, popular, fast, and easy-to-setup measurement techniques, which require only lights and cameras together with an optical path to the Region Of Interest (ROI) (which may have an arbitrary complex shape) to operate. Moreover and as presented in the last paragraph, two quantities relevant in the test-simulation dialogue perspective can be obtained: displacements as well as specimen shapes (SDIC only) over the whole ROI. These measurements are performed by relying on the matching of features contained in different pictures. On top of that, the DIC variational formulation [Fedele et al. 2013] allows the identification of shapes and displacements on an arbitrary kinematic basis. This is of great interest for test-simulation dialogue as the same basis as the one used for solving the mechanical problem associated with the model can be used (see Figure 5). The comparison between measurements and simulations is then straightforward, as defining an error requires only to subtract corresponding degrees of freedom. This naturally leads to frameworks allowing to retrieve quantities of interest either by minimising such an error (e.g. Finite Element Method Updating (FEMU) [Kavanagh & Clough 1971, Kim et al. 2007, Lecompte et al. 2007]) or by direct identification (e.g. Integrated Digital Image Correlation (IDIC) [Leclerc et al. 2009, Réthoré 1 By three-dimensional, we mean that the displacement field identified with SDIC is a 3D one and that it is measured over a (possibly) non-planar surface. It should not be mistaken with Digital Volume Correlation which extracts a displacement field over a whole volume. In the remainder of the manuscript, we thus consider surface measurements only.

2 We wish to make clear here that we do not pretend that DIC should not be coupled with other full-field measurement techniques or even punctual measurements to address the data issue in mechanical tests. We even less pretend that it is the only possible way to better understand phenomena occurring during tests. This overview should rather be taken as an introduction to the motivations for keeping on developing tools in this research field. et al. 2009, Roux & Hild 2006]). Eventually, by its optical nature, DIC appears to be well suited for performing measurements over large areas and at different resolutions. This kind of measurements should be able to provide enough data to obtain a digital representation of the particular test which is considered, a digital twin. By the expression "digital twin", we adopt a broad definition. That is, the extraction from measurements of information relevant either for the modelling (e.g. geometry, displacements, model parameters) or for the measurement itself (e.g. in DIC camera positions, focals), that allows in turn to generate data to compare simulations with.

Research outline

In this work, we would like to present and focus on some of the remaining technical obstacles preventing DIC from fully operating at larger scales, on complex geometries, and to provide more data in the test-simulation dialogue perspective:

• Obtaining more data from tests is in fact something easy to achieve by simply increasing the number of cameras in the experimental setup. However, ensuring consistency of the whole dataset in a test-simulation dialogue perspective is an open problem. Here arises the (sometimes hackneyed) topic of data assimilation or data fusion.

• A related issue is to perform measurements over a wide range of scales (or multiscale measurements). Merging measurements from different scales over the same regions and taking advantage of all available information is indeed not straightforward.

• Features (speckle patterns) used to perform image registration in DIC are often adapted to a single scale. This can lead to inhomogeneous measurement precision across cameras.

• The wider the structure, the larger the displacements for a given strain magnitude. Hence, DIC in its current form is not perfectly suitable to address multiscale measurement, as features used to perform image registration may disappear from view between two time steps, especially at the most revolved scales.

• Ensuring a homogeneous scene lighting, as required by DIC implicit assumptions, becomes harder as scale increases. This is even more true on complex geometries and/or when considering large deformations. In this context, lighting issues may result in poor measurement accuracy, or may even incapacitate DIC from operating.

This manuscript is broken down into 7 chapters. Chapter 1 presents a very condensed introduction to ill-posedness and inverse problems. As these topics will be encountered all along the present work, we wanted the reader to be familiar with them beforehand. Then, Chapters 2 and 3 deal with the photogrammetric camera calibration process and frameworks used in DIC and SDIC, especially for global approaches, and even more specifically, for Finite-Element (FE) ones. This allows to show very precisely where current limitations come from, while introducing some tools that are used later in the manuscript. In Chapter 4, multiscale approaches in DIC are presented. It encompasses two different but strongly related views: initialisation of the measurement problems and multiscale measurements, that is cameras observing the structure at different resolutions. As explained earlier, there is a current lack of patterns able to provide homogeneous measurement accuracy over a wide range of scale. For this reason, we investigate the use of fractal (or self-affine) patterns which deal for interesting candidates thanks to their statistical scale invariance properties. We show for one of them that it does result in a measurement accuracy independent of the scale [Fouque et al. 2021c]. This was an important point, before starting the development of numerical methods related to multiscale measurements, to ensure that such patterns existed. Chapter 5 stands for the main contribution of the present work. Based on Computer Vision (CV) research works, a formulation allowing to address the different identified issues such as lighting and visibility is proposed. It relies on a physical modelling of the scene and the irradiance equation as in [Fouque et al. 2021a]. However, we go a step further as we provide here a way to make use of all available information. This is a clear breakthrough compared to usual frameworks which are based on the matching of features in some of the available data and the grey level conservation equation. Then, a theoretical viewpoint is adopted, and justifications of the proposed formulation are presented in Chapter 6. There are three main arguments. The two firsts justify a posteriori the use of the formulation, as a link with usual frameworks is established. It is shown that the proposed formalism encompasses usual ones while being able to perform on a wider variety of tests [Fouque et al. 2021b]. Also, a sensitivity analysis is performed and significant improvements in terms of measurement precision compared to classical approaches is demonstrated. Finally a third a priori argument is provided based on estimation theory. We show that this formulation is optimal, provided that some assumptions, that we make explicitly, are met. Finally, a practical test case is proposed in Chapter 7 [Fouque et al. 2021a]. This brings us to the Conclusion and Perspectives.

Chapter 1. A brief introduction to ill-posedness and inverse problems

Well-posed problems

A model relies on different elements of which we will try to provide a general definition [Allaire 2005]. x will denote the data (e.g. right-hand side, geometry, initial/boundary conditions), y the solution sought (or output) and A, the operator acting on y. In a partial differential equation, for instance, A would refer to the equation itself. Hence, it may depend on some model parameters as well (denoted p). With these notations (summarised in Figure 1.1), the problem is to find y such that:

A(y, p) = x (1.1)
In the Hadamard sense, a problem is well-posed when each of the three following conditions are met:

1. for all x a solution y exists 2. for all x the existing solution y is unique 3. the solution y depends continuously on the data x Conversely, a problem is ill-posed if one (or more) of the three previous points is not satisfied.

Inverse problems

An inverse problem assumes a direct well-posed problem [Isakov 2006]. This direct well-posed problem defines the data x, output y and model A (see Section 1.1). An inverse problem associated to this direct one would not consider the solution y as an unknown anymore. Instead, y would be considered as conditions to meet to obtain the data x (first kind inverse problem, see Figure 1.2a) or to identify the set of parameters p (second kind inverse problem, see Figure 1.2b) that would have generated the solution y in the direct problem [Weisz-Patrault 2012, 1.4. Capteurs et méthodes inverses].

Remark It is convenient to distinguish two kinds of inverse problems but note that hybrid problems exist (see Chapter 2 for instance).

Most inverse problems are ill-posed. As a result, it may be difficult to obtain a solution, as the existence is not guaranteed. Thus, it is often obtained by relying on a regularisation step 1 that may be based on optimisation methods. If a set of parameters p (or data x) is obtained, it might not be unique, but, most importantly, it might not be continuous with respect to the measured "solution" y. Heuristically, it means that a small error in the measurement of y can lead to tremendous discrepancies between the "true" p (or x) and the estimated one.

After these brief theoretical prerequisites, we can move on to the first practical inverse problem of this manuscript, that is camera calibration. Chapter 2

Camera model and photogrammetric calibration

Camera models are an elementary tool in DIC and SDIC allowing to describe how a point of the three-dimensional space is captured by a camera in a two-dimensional picture. In this chapter, we provide a general overview regarding this topic, as the different steps to fully implement a photogrammetric stereo calibration process were not available in a single document in the literature reviewed. Practically, by the end of the chapter, the interested reader should be able to implement their own photogrammetric calibration process. To do so, we will first focus on a simple, distortion-free model, allowing us to go through the important steps of the transition from the world reference frame to the picture one in an easy way. Then, more evolved parameters accounting for optical system flaws (e.g. optical axis misalignment, non-coincident focal points) will be presented. Eventually, the photogrammetric calibration process is introduced and details regarding its application on real images are provided for better understanding. This allows to show that some parameters (position of the rig with respect to the target) are usually discarded, whereas we found a way to make use of them in a practical test case in Chapter 7.

Camera modelling

Because, among others, it allows to write translations, rotations, scalings and projections as matrix operations, homogeneous coordinates will be of great help in this chapter [Garcia 2001].

For instance, the affine transformation M c = R M w + t, can simply be written as M c = T M w (see Equation (2.1)) which makes it more compact and reduces the amount of notations in what follows. A significant part of notations used in this chapter are presented in Figure 2

.1. v u θ Y c X c Z c y r x r Y w X w Z w T P f A m f O R c R r R pix R w O w c M O r O pix

The pinhole model

From the world coordinate system to the camera coordinate system Let M be a point of the three-dimensional space. We will denote M w = (X w Y w Z w ) T the coordinates of M relative to R w coordinate system. Typically, in our calibration process, R w will be the same as the target or the ROI one. In order to easily do the calculations in the following subsection, we will express M coordinates in the camera coordinate system

R c (M c = (X c Y c Z c ) T ) with respect to M w : M c = T M w with M c = (X c Y c Z c 1) T , M w = (X w Y w Z w 1) T and T = R t 0 T 1 . (2.1)
Remark From a more practical viewpoint t stands for the translation vector going from

R c reference point O c to R w reference point O w expressed in the camera coordinate system R c (t = (O c O w ) c ).
In the following, we will note t = (t x t y t z ) T and for R, rotation matrix of T :

R = (R ij ) (i,j)∈[[1,3]] 2 = R z (γ) R y (β) R x (α), (2.2)
where:

R x (α) =    1 0 0 0 cos(α) -sin(α) 0 sin(α) cos(α)    ; R y (β) =    cos(β) 0 sin(β) 0 1 0 -sin(β) 0 cos(β)    ; R z (γ) =    cos(γ) -sin(γ) 0 sin(γ) cos(γ) 0 0 0 1    . (2.3) v u Y c X c Z c y r x r f X c Y c x r y r M m Z c O c R c Figure 2
.2: Schematic drawing showing the two green triangles in which Thales' theorem will be used.

From the camera coordinate system to the retinal plane Once M c is known, it is possible to express its projection on the retinal plane by assuming a perfect optical system. This process is shown in Figure 2.2.

Remark

Contrary to what is suggested in Figure 2.2, we will assume f |Z c |. As a result, we will make no difference between focal and retinal planes, and assume that the projection m of a point M is always located at a distance f from R c reference point O c along Z c .

Applying Thales' theorem in the two green triangles in Figure 2.2, we get:

         x r f = - X c Z c y r f = - Y c Z c , (2.4)
with the convention f > 0, Z c < 0. As a result, we get in homogeneous coordinates:

m r =    -f 0 0 0 0 -f 0 0 0 0 1 0    M c . P f (2.5)
Remark Note that, here, homogeneous coordinates require a renormalisation to actually get x r and y r from Equation (2.5).

From the retinal plane to the picture coordinate system This step describes the way the camera photosensor samples the image located in the retinal plane and makes it a picture. Note that, in the retinal plane, m r coordinates are in meters, whereas in the picture they are in pixels. Therefore we need to introduce two sampling coefficients, α u and α v , allowing to get image coordinates along u and v directions respectively (see Figure 2.1 again). Their dimension is in pixels per meter, and α u samples the x r direction while α v samples the y r direction. With these conventions it is quite straightforward to write the relation between m r and m pix . But first two other parameters (whose dimension is in pixels) have to be introduced to the reader: u 0 and v 0 , the optical centre pixel coordinates. They stand for the coordinates of the two-dimensional translation vector going from R pix reference point O pix to R r one O r , expressed in R pix :

m pix =    α u 0 u 0 0 α v v 0 0 0 1    m r . A (2.6)
The projection matrix As explained in [Garcia 2001], the product A P f can be rewritten:

A P f =    α u 0 u 0 0 α v v 0 0 0 1       -f 0 0 0 0 -f 0 0 0 0 1 0    =    -α u f 0 u 0 0 0 -α v f v 0 0 0 0 1 0    =    -α u f 0 u 0 0 -α v f v 0 0 0 1       1 0 0 0 0 1 0 0 0 0 1 0    A P . (2.7)
This algebra shows that it is not possible to differentiate a camera with a focal length f and sampling parameters α u and α v from a camera of focal length 1 and sampling parameters α u f and α v f . Since the triplet (1/f, α u , α v ) is defined up to a scale factor, we are only able to identify two parameters out of the three initial ones. We choose to keep the products α u f and α v f that we will denote respectively f x and f y . Also, since there is an infinity of pairs (A, P f ) for a same camera, it does not make sense to keep on distinguishing them. Therefore, we consider in the following

K = A P f = A P =    -f x 0 u 0 0 0 -f y v 0 0 0 0 1 0   , the projection matrix.
Finally, the product of K and T allows to write:

       u = -f x X c Z c + u 0 = -f x R 11 X w + R 12 Y w + R 13 Z w + t x R 31 X w + R 32 Y w + R 33 Z w + t z + u 0 v = -f y Y c Z c + v 0 = -f y R 21 X w + R 22 Y w + R 23 Z w + t y R 31 X w + R 32 Y w + R 33 Z w + t z + v 0 , (2.8)
where the components of the rotation matrix R (see Equation (2.2)) are (R ij ) (i,j) . Figure 2.3 sums up previous notations and steps as a block diagram.

Remark From Equation (2.8), a physical meaning can be given to f x (resp. f y ): it stands for the number of pixels taken in the u direction (resp. v direction) by an object measuring one meter in the X c direction (resp. Y c direction) located at a distance of one meter of the camera centre

O c along Z c . M w M c m r m pix T K t x t y t z ¡ β γ ( ) ( f ) u 0 v 0 α u α v ( )

Distortion models

In this section two different types of distortions will be presented: linear and non-linear ones. By distortions, we mean sensor geometrical imperfections resulting in discrepancies with respect to the pinhole model.

Linear distortion or skew

u = u v = cos(θ)u + sin(θ)v
.

By substituting u and v in m pix expression and identifying u and v, we get:

u = u + v cos(θ) v = v sin(θ) .
Finally, expressing u and v with respect to u and v:

u = u -v cot(θ) v = v/ sin(θ) .
(2.9)

Using Equation (2.8) in Equation (2.9), we can now obtain the new expression of the projection matrix K:

       u = -f x X c Z c + f y cot(θ) Y c Z c + u 0 -v 0 cot(θ) v = -f y / sin(θ) Y c Z c + v 0 / sin(θ)
.

(2.10) It holds that:

K =    -f x f y cot(θ) u 0 -v 0 cot(θ) 0 0 -f y / sin(θ) v 0 / sin(θ) 0 0 0 1 0    .
(2.11)

Practically, we often have θ ∼ π 2 . Therefore it is relevant to linearise the expression of K. With θ = π 2 + δθ, δθ 1, we then get:

K ∼    -f x -f y δθ u 0 + v 0 δθ 0 0 -f y v 0 0 0 0 1 0    .
(2.12)

Remark δθ is often called the skew parameter.

Non-linear distortions This section will be based to a large extent on the corresponding section of [Garcia 2001]. Here non-linear distortions are meant to describe the different types of geometrical aberrations arising from camera optical system flaws. Therefore the modification in the process described in Section 2.1.1 will take place in the retinal coordinate system, before the camera sensor sampling. As a result, we define

m d = m r + δ(m r ) given that m pix = A m d . δ x (m r ) = x r r 1 ρ 2 + r 2 ρ 4 + r 3 ρ 6 + 2d 1 x r y r + d 2 3x 2 r + y 2 r + p 1 ρ 2 δ y (m r ) = y r r 1 ρ 2 + r 2 ρ 4 + r 3 ρ 6 + 2d 2 x r y r + d 1 x 2 r + 3y 2 r + p 2 ρ 2 Radial Decentering Prismatic with ρ 2 = x 2 r + y 2 r .
(2.13)

The effect on the pictures of previous parameters can be seen in Figure 2.5. Non-parametric approaches allowing to account for distortions also exist [Cornille 2005]. An interpolation of remaining reprojection errors is performed to recover the distortion shape. However, as we will emphasise in the present manuscript, we believe that relying on physical approaches, such as the parametric one presented previously, is more satisfying. Admittedly, this viewpoint is not very fashionable at the moment, but we believe that relying on a modelling step, when possible, makes explicit model validity domains, allows to criticise the assumptions made, and enables us to derive uncertainties. To qualify previous statement, it may be hard in some cases (e.g. nonoptical methods, complex image formation process) to derive a parametric model. However, it is important to take into account the image formation process and include the non-parametric model at a relevant position in the block model.

Remark

Even though this is a non-linear operation, we can write this step:

m d = ∆(m r ) m r with ∆(m r ) =    1 0 δ x (m r ) 0 1 δ y (m r ) 0 0 1    . (2.14)
With such notations, m pix can be retrieved:

m pix = A ∆ P f T M w ,
where the projection matrix K can be identified. Noting indeed that

P f P T 1 f =    1 0 0 0 1 0 0 0 1    = 1, (2.15)
we have: In the end, we built a camera model in which different kind of parameters need to be identified:

m pix = K P T 1 f ∆ P f T M w ∆ with ∆ =      1 0 -δ x /f 0 0 1 -δ y /f 0 0 0 1 0 0 0 0 0      . ( 2
• Six extrinsic parameters: p ext = (t x t y t z α β γ) T .

• Four intrinsic parameters:

p int = (f x f y u 0 v 0 ) T . • Up to eight distortion parameters: d = (δθ r 1 r 2 r 3 d 1 d 2 p 1 p 2 ) T .
These parameters are collected in a camera model parameter vector p = p ext T p int T d T T , and the whole process described above, can be summed up by the camera model P :

P : W × P → Π (M w , p) → m pix ,
where W and Π denote respectively the physical space and the image plane. P stands for the camera parameter space and 10 ≤ dim P ≤ 18 (dim P = 10 with no distortion and dim P = 18 for the full distortion model).

Remark

In what follows, making the distinction between all these reference frames is not relevant anymore, we will thus simply denote M w and m pix respectively by X and x.

Photogrammetric calibration process

As explained in Section 2.1, a camera model relies on different kind of parameters that we need to identify. This identification problem is an inverse problem (which may not be classified as first nor second kind). If we know the intrinsic parameters (and possible distortions) of a camera (p int , d) as well as its position (p ext ) with respect to a set of points, it is straightforward, thanks to Section 2.1, to compute, for each position (X), the associated projection (x) in the corresponding picture. It is much more difficult though, from a set of pictures of the same set of points, to get the parameters (extrinsics, intrinsics and possible distortions) that generated these pictures. Given its inverse nature, this problem is an ill-posed one (for different reasons, e.g. inability of the camera model to account for the whole physics, picture noise or deformation of the set of points between picture captures, point detection precision). As a result, the existence of a parameter set allowing to obtain the exact same projections is not guaranteed.

To circumvent these issues, the problem will be formulated as a minimisation. That is, instead of looking for the parameter set allowing to obtain the exact same projections, we will look for the set generating the projections which are as close as possible to the detected ones. A functional standing for the distance between the detected points and the generated ones will be minimised using an iterative, gradient-based, optimisation algorithm. It will require an initialisation (first set of parameters) leading to a low value of the functional (in order to benefit from the local convergence of the algorithm).

What is rather surprising is that parameters of interest (intrinsics and distortions) are quite easy to initialise in a first approach. On the contrary, it is rather difficult to estimate extrinsic parameters from a target picture, while we do not particularly want to know them precisely at this stage.

Parameters initialisation

Intrinsic parameters and distortions

As mentioned before, we can easily estimate intrinsic parameters, with a single picture:

• u 0 and v 0 , respectively, with the half-height and half-width of the image (in pixels).

• Knowing approximately the target size and its distance with respect to the camera allows to get f x from a cross-multiplication (see Equation (2.4)) ; then, one can assume in addition f x = f y . Alternatively, it can be obtained by dividing the focal length f of the optical system by the size of a pixel (on the image sensor).

• Distortion vector d = 0 is often a good first guess.

Extrinsic parameters

To estimate the extrinsic parameters, our approach will be based on the method developed by [Ravn et al. 1994]. As we need to establish the equations of this paper for our choice of reference frames, we will use the same steps as in [Horaud & Monga 1995].

Equations establishment Denoting, for the rotation matrix, the line vectors

R i = (R ij ) j∈[[1,3]] = (R i1 , R i2 , R i3 ), we have thanks to Equation (2.8):    su sv s    =    -f x R 1 + u 0 R 3 -f x t x + u 0 t z -f y R 2 + v 0 R 3 -f y t y + v 0 t z R 3 t z         X w Y w Z w 1      , (2.17)
which is consistent with d = 0 as first guess. Since our target can be considered as a plane in (X w , Y w ), we will assume Z w = 0 in the following (and we change T accordingly). As a result, Equation (2.17) can be rewritten (considering s = s tz ):

   s u s v s    =    -fxR 11 +u 0 R 31 tz -fxR 12 +u 0 R 32 tz -fxtx+u 0 tz tz -fyR 21 +v 0 R 31 tz -fyR 22 +v 0 R 32 tz -fyty+v 0 tz tz R 31 tz R 32 tz 1       X w Y w 1    .
(2.18)

From the last line of Equation (2.18), we get s = R 31 tz X w + R 32 tz Y w + 1 which can be first multiplied by u and also by v. This yields two equations:

       u = s u - R 31 t z X w u - R 32 t z Y w u v = s v - R 31 t z X w v - R 32 t z Y w v . (2.19)
Finally, thanks to Equations (2.18) and (2.19) and by denoting

κ = (κ 1 κ 2 κ 3 κ 4 κ 5 κ 6 κ 7 κ 8 ) T =                 -fxtx+u 0 tz tz -fyty+v 0 tz tz -fxR 11 +u 0 R 31 tz -fxR 12 +u 0 R 32 tz -fyR 21 +v 0 R 31 tz -fyR 22 +v 0 R 32 tz R 31 tz R 32 tz                
, we can write in the same way as in [Ravn et al. 1994]:

u v = 1 0 X w Y w 0 0 -X w u -Y w u 0 1 0 0 X w Y w -X w v -Y w v κ.
(2.20)

From at least four points of which both positions in the world coordinate system (X w , Y w ) and in the picture (u, v) one are known, we can get all eight κ components. With four wellchosen points it is possible to invert the resulting 8 × 8 matrix. With more than four points, different methods for solving overdetermined linear systems can be used (e.g. Singular Value Decomposition, least squares).

Remark As the method with only four points produced satisfactory results as an initialisation, we did not implement more sophisticated strategies in the present work.

Remark Because close points or nearly-aligned points can lead to ill-conditioned 8×8 matrices, the four points which we chose were the four target corner points.

Parameters estimation

Knowing approximately intrinsic parameters (typically from an a priori estimation as in Section 2.2.1), [Garcia 2001] provides an effective manner to obtain T from κ. Using R orthonormality properties, we can write:

R 2 11 + R 2 21 + R 2 31 + R 2 12 + R 2 22 + R 2 32 = 2.
Substituting κ components and intrinsic parameters for R ij :

t z = - 2 κ 3 -u 0 κ 7 fx 2 + κ 5 -v 0 κ 7 fy 2 + κ 4 -u 0 κ 8 fx 2 + κ 6 -v 0 κ 8 fy 2 + κ 2 7 + κ 2 8 .
(2.21)

Remark Writing t z < 0, we implicitly made the assumption the object was in front of the camera, according to Figure 2.1, which seems reasonable.

Then, all other values can be computed from κ and Equation (2.21):

                                                         t x = u 0 t z -κ 1 t z f x t y = v 0 t z -κ 2 t z f y R 11 = t z u 0 κ 7 -κ 3 f x R 12 = t z u 0 κ 8 -κ 4 f x R 21 = t z v 0 κ 7 -κ 5 f y R 22 = t z v 0 κ 8 -κ 6 f y R 31 = t z κ 7 R 32 = t z κ 8 . (2.22)
Denoting, for the rotation matrix, the row vectors

R i = (R 1i R 2i R 3i ) T , i ∈ [[1, 3]], we finally get R 3 = R 1 ∧ R 2 .
But as explained in [Garcia 2001], Ravn's method does not guarantee the orthonormality of R. Therefore, it can be a good idea to compute its projection on the set of matrices of M 3,3 (R) of which the determinant equals one. This can easily be done by computing the Singular Value Decomposition (SVD): it allows to write R = U Λ V T . Substituting the identity matrix 1 to Λ is equivalent to computing the aforementioned projection. Finally, we need to compute α, β and γ from R. Denoting, for instance, sin(α) = s α and cos(β) = c β , from Equations (2.2) and (2.3), we can write:

R (α, β, γ) =    c β c γ s α s β c γ -c α s γ c α s β c γ + s α s γ c β s γ s α s β s γ + c α c γ c α s β s γ -s α c γ -s β s α c β c α c β    . (2.23) Since for every (α, β, γ) , R (α, β, γ) = R (π + α, π -β, π + γ)
, in order to guarantee the existence of a unique set α, β, γ for every rotation, we adopt the convention (α, -π, π]. Thus, from Equation (2.23), we have:

β, γ) ∈ [-π, π] × [0, π] × [
         β = -arcsin(R 31 ) s β = -R 31 c β = 1 -R 2 31 .
(2.24)

From Equations (2.23) and (2.24), we then get:

         c γ = R 11 c β s γ = R 21 c β , (2.25)
and

         c α = R 33 c β s α = R 32 c β . (2.26)
Finally, from Equations (2.25) and (2.26) it is easy to compute α and γ.

Point positions

Considering point space positions as unknowns will allow us to use rather imprecise targets, printed ones for instance. This kind of targets is much more flexible and cheaper than those where the point relative positions must be precisely known.

As soon as the target can be considered as rigid the whole time between the first and the last picture, the implemented algorithm will be able to reevaluate point relative positions from an initial one. In order to initialise these positions, we can simply measure them approximately in the chosen target reference frame. Now we will present the framework allowing to refine all previous initialisations at once, the photogrammetric calibration process. This may also be called bundle adjustment 1 .

Process for a single camera

As explained before, the calibration process will be based on the minimisation of a functional. More precisely, the idea is to minimise the error between the measured pixel coordinates of 1 "Bundle adjustment is the problem of refining a visual reconstruction to produce jointly optimal 3D structure and viewing parameter (camera pose and/or calibration) estimates. Optimal means that the parameter estimates are found by minimizing some cost function that quantifies the model fitting error, and jointly that the solution is simultaneously optimal with respect to both structure and camera variations. The name refers to the 'bundles' of light rays leaving each 3D feature and converging on each camera centre, which are 'adjusted' optimally with respect to both feature and camera positions. Equivalently -unlike independent model methods, which merge partial reconstructions without updating their internal structure -all of the structure and camera parameters are adjusted together 'in one bundle'. Bundle adjustment is really just a large sparse geometric parameter estimation problem, the parameters being the combined 3D feature coordinates, camera poses and calibrations." Citation from [Triggs et al. 2000] points in a picture and the estimated pixel coordinates through the camera model. In practice, a target with N p points will be used. N i pictures of this target will be taken by the camera. As a result, as introduced in [Garcia 2001], 2N i N p equations will be available. At the same time, there are 4 intrinsic unknowns, 8 distortion parameters (at most), 6N i extrinsic unknowns, and 3N p unknown tridimensional coordinates.

However, in order to reduce functional kernel dimension, rigid body motions of the {camera + target} system should be avoided. To do so, the easiest way is to prescribe the tridimensional position of 2 points in the world coordinate system as well as one of the three coordinates of a third point. Without loss of generality, we can assume that the first point coordinates are (0 0 0) T . Rigid translations of the system are now forbidden. Then, a second point coordinates can be enforced as follows: (L 0 0). Obviously if there are points in the target of which the distance uncertainty is minimal, it would make sense to choose among these points, as the precision on L measurement drives the entire calibration process precision. The only remaining rigid body motion is the rotation about the axis defined by the two previous points. Finally, and again without loss of generality, since three points are always coplanar, assuming Z w = 0 for a third point allows to prevent any rigid body motion of the {camera + target} system. Obviously, this point should not belong to the line defined by the first two points. Eventually, this leads to 12 unknown intrinsic parameters, 6N i extrinsic unknowns and 3N p -7 unknown tridimensional coordinates together with 2N i N p equations, N i ≥ 1, N p ≥ 3. In order to have at least an overdetermined set of equations, we will need:

2N i N p ≥ 5 + 6N i + 3N p ,
or equivalently:

N p (2N i -3) ≥ 5 + 6N i .
(2.27)

From Equation (2.27), we get that:

N i ≥ 2 ⇔ ∃ N p ∈ N, 3 ≤ N p ≤ 17, 2N i N p ≥ 5 + 6N i + 3N p .

Image points detection

In our framework, we used a target made of a coding system (allowing to know the orientation of the target) and circular holes. The library opencv (or more precisely its Python binding) provides the method findCirclesGrid which allowed us to detect the subpixel coordinates of target holes in the picture. Using then the coding system to make sure that for every point j ∈ [[0, N p -1]], the index j referred always to the same physical point in the target for every picture, we obtained xij

(i,j)∈[[0,N i -1]]×[[0,Np-1]]
, the table of the measured coordinates for every point and every picture in the picture coordinate system.

Functional construction

The residual for picture i and point j reads:

r ij = xij -P (X j Y j Z j ) T , p ext i T p int T d T T .
(2.28) with, for instance, X = (X j ) j∈ [[0,] the vector containing the coordinates in the X w direction of every point. p ext i is the vector of extrinsic parameters of the world coordinate system with respect to the camera for picture i. Then, the picture i residual can be recast as:

F i = r T i0 r T i1 . . . r T i(Np-1) T .
(2.29)

Finally the total residual reads:

F = F T 0 F T 1 . . . F T N i -1 T .
(2.30)

In the end the functional minimised by Levenberg-Marquardt algorithm is [Garcia 2001]:

J = F T F . (2.31)
Remark Because of the implementation of the optimisation algorithm (Levenberg-Marquardt), we considered a functional whose physical meaning is not obvious. However if we consider:

σ = 1 N i N p -1 F T F = 1 N i N p -1 N i -1 i=0 Np-1 j=0 r T ij r ij .
(2.32)

σ can be seen as the error standard deviation, in pixel.

Functional argument

We wish to formulate our problem as follows:

Find θ such that, θ = arg min θ F (θ) T F (θ) .

(2.33)

Genuinely, we would define

θ = t x T t y T t z T α T β T γ T p int T d T X T Y T Z T T (2.34)
with, for instance, t x = (t x,i ) i∈[[0,N i -1]] the vector collecting the N extrinsic parameters of translation in the X c direction corresponding to the N i pictures. Now, as mentioned in the introduction of this section, in order to avoid rigid body motion, we will prescribe 7 coordinates among the 3N p coordinates of the N p target points. We will denote respectively X , Y and Z the set of points whose X, Y or Z coordinate is constrained. And we will recall:

θ = t x T t y T t z T α T β T γ T p int T d T X T Y T Z T T (2.35) with X = (X j ) j∈[[0,Np-1]]\X , Y = (Y j ) j∈[[0,Np-1]]\Y and Z = (Z j ) j∈[[0,Np-1]]\Z .
To make it more practical, using a 9 × 12-hole target, we chose to have X = {0, 11}, Y = {0, 11} and Z = {0, 11, 99}, and we looked for:

θ = arg min θ F T F θ , X 0 = Y 0 = Z 0 = Y 11 = Z 11 = Z 99 = 0, X 11 = L T .
(2.36)

Multiple-camera setup

In the following, we will denote N c the number of cameras of the considered setup.

When dealing with several cameras, with a stereo-rig for instance, there are a few modifications to implement compared to the previous framework. First, when counting the number of unknowns and number of available equations, we have, at most 12N c intrinsic unknowns (with distortions) and still 3N p unknown tridimensional coordinates. However, if we consider a rigid setup, there are 6N i extrinsic unknowns for the position of the world reference frame with respect to a reference camera and 6(N c -1) extrinsic unknowns (relative position of other cameras with respect to a reference one). With notations of Figure 2.7, it means that if we know T (for a given picture or, equivalently, target position) and T s (which does not depend on the considered picture), we can get T = T s T . It should not be forgotten that 7 points coordinates must be constrained to avoid rigid body motions of the {camera+target} system. In total, there are thus 12N c + 6N i + 3N p + 6(N c -1) -7 = 18N c + 6N i + 3N p -13 unknowns. At the same time, the number of equations is 2N i N p N c . It means that with a given number of cameras, an additional camera increases the number of unknowns by 18 (12 intrinsic and 6 extrinsic parameters) while providing at the same time 2N i N p new equations. Therefore, as soon as N i N p > 9, it helps regularising the problem.

Remark Even in the case N i = 1, the condition N i N p > 9 is easily met and thus adding cameras generally provides a strong regularisation.

ij (i,j)∈[[0,N i -1]]×[[0,Np-1]]
the table of the measured coordinates for every point and every picture in the picture coordinate system of camera c.

Analogously to Section 2.2.2, we define r c ij the vectorial pixel residual of camera c, point j, image i:

r c ij = xc ij -P p ext c i T p int c T (d c ) T T , (X j Y j Z j ) T .
(2.37)

We will denote F c the vectorial residual of camera c (see Equation (2.30)). With this convention the total residual reads:

F = F 0 T F 1 T . . . F Nc T T . (2.38)
And the functional minimised by the optimisation algorithm is:

J = F T F (2.39)
Remark As in Section 2.2.2, we can then define:

σ = 1 N i N p N c -1 F T F = 1 N i N p N c -1 N i -1 i=0 Np-1 j=0 Nc-1 c=0 r c ij T r c ij (2.40)
Functional argument In the following, we will denote the functional argument Θ. Genuinely, we would define Θ as θ 0 θ 1 . . . θ Nc , with θ c the argument of F c for camera c as defined in Equation (2.35). But as mentioned before, we do not need T transformation for every single picture as soon as c > 0. That is the reason why we define:

θ c = t x 0 T t y 0 T t z 0 T α 0 T β 0 T γ 0 T p int 0 T d 0 T T if c = 0 t c x,s t c y,s t c z,s α c s β c s γ c s p int c T (d c ) T T if c > 0 (2.41)
Remark Note that extrinsic parameters are vectors (containing parameter values for every picture) only if c = 0. As a result, for every camera residual F c , c > 0, we need to compute θ c from θ c . In order to reconstruct extrinsic vectors for every picture i, we can use the identity T i = T s T i , with T s parameters contained in θ c extrinsic parameters and T i parameters contained in θ 0 = θ 0 i th component of extrinsic parameters. It can also be written:

R i = R s R i and t i = R s t i + t s (2.42)
Remark The initialisation process described in Section 2.2.1 does not apply for t c x,s t c y,s t c z,s α c s β c s γ c s ) , c > 0 because it represents extrinsic parameters for the transformation between the reference camera and camera number c (not between the target and camera number c). In order to estimate these parameters, one can derive from Equation (2.42) that R s = R i R T i and t s = t i -R s t i . With i = 0 for instance, we can estimate T 0 and T 0 with the same procedure as in Section 2.2.1 and then get an initialisation for T s .

And the optimisation parameter is:

Θ = θ 0 θ 1 . . . θ Nc X Y Z .
(2.43)

It allows to write the final form of the calibration problem:

Θ = arg min Θ F T F (Θ, X 0 = Y 0 = Z 0 = Y 11 = Z 11 = Z 99 = 0, X 11 = L) .
(2.44)

Remark As explained in the beginning of this section, from this optimisation process, we get useful optimised values for each camera intrinsics as well as relative camera positions. The relative position of the target with respect to the camera rig (or the camera) is usually discarded.

Example on an actual stereo rig

This section aims at illustrating and validating the process described above. To do so, two five-megapixel cameras and a 9 × 12-hole target were used. A set of 16 target positions was captured by both cameras.

Initialisation

Using the method introduced in Section 2.2.1 and in Section 2.2.3, we were able to initialise extrinsic and intrinsic parameters for both cameras. The first step consists in detecting points in every picture and associate for every physical point an integer p. One can see in Figure 2.8 the detection of the target reference frame that allows to know the target orientation and to implement this association. In Figures 2.9 and 2.10, the initialisation is shown. It means that with the parameters estimated with Section 2.2.1 and Section 2.2.3, we were able to evaluate points positions in the picture as the yellow dots. As a result, we got an error standard deviation a bit above one pixel for camera 0 (see Figure 2.9) as well as about 4 pixels for camera 1 (see Figure 2.10). These results show that the initialisation method is very efficient.

Remark

We expected σ 1 > σ 0 since the extrinsics initialisation spanning all pictures is performed only for camera 0 and since only one picture was used to estimate T s .

Results

In this section, from previous initialisation, we apply the photogrammetric calibration process in two different ways:

• process applied to each camera independently, that is two different optimisations,

• multicamera calibration process, that is only one optimisation.

As expected, Table 2.1 shows that error standard deviations are greater when less degrees of freedom are available (with N i pictures and N c cameras, for a method with N c independent calibrations there are 6(N c -1)(N i -1) additional unknowns compared to a method with a single optimisation). However, as already discussed, it shows the benefit to regularise the calibration process, and thus reduces the uncertainties associated with the identified parameters [Garcia 2001].

Remark 17 pictures were actually available for the calibration procedure. However, the target detection step failed for one of them and 2 pictures were excluded because they exhibited large residuals compared to the others. This explains why in Figure 2.11, no point is associated with picture indices 2, 12 and 13 and why 16 pictures are shown in Figures 2.9 and 2.10.

Method

σ 0 σ 1 σ Two optimisations 0.06891 0.08572 0.07778 One optimisation 0.07837 0.09332 0.08617 Table 2.1: Error standard deviation for two different methods : two independent optimisation processes (one for each camera) as described in Section 2.2.2 and a single process as in Section 2.2.3. Error standard deviation is presented for each camera (σ 0 and σ 1 ) as well as for the whole system (σ).

Conclusion

This chapter has provided the opportunity to detail how we implemented an in-house calibration software. We have carefully gone through the different steps of stereo photogrammetric calibration and have presented it as fully as possible, so that a similar implementation may be performed. These developments, made in an early phase of this thesis, turned out to be extremely useful in the practical application proposed in Chapter 7. Of course, more sophisticated processes could be derived, accounting for target point visibility or allowing to calibrate a virtual rig made of a single camera at different positions for instance. i (error standard deviation in pixels for picture i, camera c) at two different stage of the calibration process, before and after optimisation. Error standard deviations are also presented for two independent optimisation processes (one for each camera) as described in Section 2.2.2 and for a single optimisation process as described in Section 2.2.3.

Eventually, it was important to present camera models, before diving into full-field measurements in the next chapter. Cameras stand indeed for a key element in the DIC and SDIC frameworks, as they are used to generate data stemming from tests.

Chapter 3

Full-field measurements in Digital Image Correlation

The DIC and SDIC frameworks are presented. First, the two-dimensional one (DIC) is introduced. It is explained to which extent retrieving a displacement field from pictures at various levels of deformation stands for an ill-posed inverse problem. The Ritz-Galerkin method used to circumvent this issue is detailed and the use of camera models from Chapter 2 is included. Details regarding the implementation are provided. Then, the extension of this framework to the SDIC one, which allows to obtain a three-dimensional surface displacement field, is presented. This allows to start discussions regarding the soundness of this extension and the use of SDIC on large-scale applications. In the context of mechanical tests, DIC is a method allowing to identify a two-dimensional displacement field. This identification is performed over a planar surface, the ROI, from frontoparallel pictures. The knowledge of this displacement field allows then to derive quantities of interest (e.g. strains, displacement maximum amplitude or stress if a model is available) or to adjust for material parameters in a second step by relying on FEMU [Kavanagh & Clough 1971, Kim et al. 2007, Lecompte et al. 2007], or even directly, thanks to IDIC [Leclerc et al. 2009, Réthoré et al. 2009, Roux & Hild 2006]. A more exhaustive list of identification methods is presented in [Colantonio 2020[START_REF] Grédiac | The use of full-field measurement methods in composite material characterization: interest and limitations[END_REF]].

Typically, one would consider two pictures of a same object at different levels of deformation: a reference state one and a deformed state one. They are often respectively denoted by f and g. Based on the matching of features in the images, the idea is then to find the displacement field u allowing to recover f from g (see Figure 3.1). Because of homogeneous reflectance properties of materials used in structural mechanics (e.g. aluminium, steel) patterns are often deposited on test sample surfaces to facilitate this matching.

Compared to standard measurement devices, such as strain gauges for instance, DIC exhibits features that makes it belong to a separate sort of measurement techniques. The main drawbacks associated with this technique is that, depending on a wide variety of causes (e.g. operator, test setup, patterns) measurement uncertainty can be strongly affected. For this reason, guidelines have been defined to help standardise the measurement process [International DIC Society et al. 2018]. However DIC has also some strong advantages provided that an optical path to the ROI can be guaranteed all along the experiments. The first advantages can be expressed in terms of cost, as strain gauges are consumable items, and time, as each gauge requires to be precisely glued and connected to the acquisition system. In addition, strain gauges obviously give only access to sparse information. Also, depending on the targeted measurement point density, and test setups, the additional mass associated with strain gauges may affect the structural dynamic response (see for instance the number of cables in Figure 3). Eventually, a typical length scale is introduced by the size of strain gauges, which averages the strain measurement over a small surface patch.

In contrast, DIC makes use of cameras and associated optical systems together with lights, which can be depreciated over the periods they are used and require a shorter setup time. This optical nature provides a much greater measurement point density and does not interfere with test setups. Also, the measurement uncertainty is strongly related to the physical length corresponding to one pixel and, depending on the optical system magnification, a range of scales as wide as imaginable may be instrumented. This makes DIC a scale-free measurement method. For this reason, DIC is used, for instance, from sub-grain measurements in polycrystalline materials [START_REF] Stinville | [END_REF] to earth surface measurements [Bickel et al. 2018].

Two-dimensional Digital Image Correlation

Grey level conservation or strong formulation

DIC is based on the grey level conservation hypothesis [START_REF] Horn | [END_REF]. Over a domain S in the image f , the problem writes: where the sought displacement field is u = (u v) T and S is a subset of the image plane I of f . This is an ill-posed inverse problem [Besnard et al. 2006]. The main reason is that for each x, there is only one equation for two unknowns (u components). Therefore, only the projection of u along the image gradient may be recovered locally. But the problem is even more pathological as pictures f and g are taken with a digital sensor. This results in a discrete set of data, the pixels, when the sought field u is defined at each point of the ROI. Moreover, these data do not vary continuously as they are encoded (e.g. 8-bit encoding corresponding to 256 values). Eventually, the pictures are corrupted by noise during the acquisition process. Hence, given f and g, the existence of a displacement field u satisfying Equation (3.1) is not guaranteed.

Variational formulation

To relax the problem, it is possible to reformulate it as a least-squares problem. That is, minimising the integral of squared residuals of Equation (3.1):

Find u ∈ L 2 (S) 2 , u = arg min u F (u), F (u) = S f (x) -g (x + u(x)) 2 dx. (3.2)
This allows to address noise issues, as Equation (3.1) is not enforced in a strong way anymore. However, Equation (3.2) is a non-linear least-squares minimisation. To solve it, [Fedele et al. 2013] suggests proceeding iteratively, namely by defining

u (i+1) (x) = u (i) (x) + du (i) (x)
and developing F around u (i) . In order to reduce the amount of notation du (i) , for instance, will stand for du (i) (x) in the following. Also, we define:

φ u (i) : x ∈ S → x + u (i) (x). (3.3)
In this case, we can write:

F (u (i+1) ) = F (u (i) ) -2 S f -g • φ u (i) du (i) T ∇g • φ u (i) + S du (i) T ∇g • φ u (i) ∇g T • φ u (i) -f -g • φ u (i) H g • φ u (i) du (i) +o du (i) 2 , (3.4)
where, for instance, ∇g and H g respectively denotes the gradient and the Hessian matrix of g, and • denotes the usual function composition.

As we want to minimise F , the stationarity condition around du (i) at the second order writes dF (u (i) + du (i) + tdv) dt t=0 = 0 :

∀ dv ∈ L 2 (S) 2 , S dv T ∇g • φ u (i) ∇g T • φ u (i) -f -g • φ u (i) H g • φ u (i) du (i) -f -g • φ u (i) dv T ∇g • φ u (i) = 0.
(3.5)

Remark

The exact same formulation as Equation (3.5) can be obtained by writing a Newton scheme for ∇F = 0 [Passieux & Bouclier 2019].

Because the Hessian term in the second order development is multiplied by the residual f -

g • φ u (i)
, which is assumed to be small, one can neglect it in the iteration scheme. In this case, the problem can be rewritten as a linear least-squares minimisation:

Find

du (i) ∈ L 2 (S) 2 , du (i) = arg min du (i) F (i) lin (du (i) ),
with

F (i) lin (du (i) ) = S f -g • φ u (i) -du (i) T ∇g • φ u (i) 2 .
(3.6)

The stationarity condition around du (i) writes dF

(i) lin (du (i) + tdv) dt t=0 = 0 : ∀ dv ∈ L 2 (S) 2 , S -2 f -g • φ u (i) dv T ∇g • φ u (i) +2dv T ∇g • φ u (i) ∇g T • φ u (i) du (i) = 0. (3.7)
Finally the DIC variational formulation can be expressed as follows, which exactly yields to neglecting the Hessian term in the second order development of Equation (3.5):

Find

du (i) ∈ L 2 (S) 2 , ∀ dv ∈ L 2 (S) 2 , a (i) S (du (i) , dv) = l (i) S (dv) (3.8) with      a (i) S (du (i) , dv) = S dv T ∇g • φ u (i) ∇g T • φ u (i) du (i) l (i) S (dv) = S f -g • φ u (i) dv T ∇g • φ u (i)
.

(3.9)

Note that a

(i)
S is a positive symmetric bilinear form, provided that the images possess a regular texture (i.e. the gradient may vanish exclusively over a null measure subset [Fedele et al. 2013]), and l (i) S a linear one.

Remark Previous equation may be obtained from the usual Gauss-Newton scheme which is almost always used in case of non-linear least-squares minimisation (see Equation (3.2) for instance) [Passieux & Bouclier 2019].

Ritz-Galerkin method

General formulation

As mentioned earlier, Equation (3.1) locally provides information about u only along the image gradient direction. Previous developments did not tackle this issue.

Heuristically, to solve Problem (3.1), u cannot be defined at each position x ∈ S but it can be for subdomains of S, aggregating multiple positions. Mathematically, this idea is achieved by looking for an approximation du (i) h of the displacement field sought du (i) in a finite dimensional subspace V h ⊂ L 2 (S)

2 such that u

(i+1) h = u (i) h + du (i) h . In V h , starting from u (0) = u (0) h ∈ V h , equation (3.8) reads: Find du (i) h ∈ V h , ∀ dv h ∈ V h , a (i) S (du (i) h , dv h ) = l (i) S (dv h ) (3.10) Let N j j∈[[1;N dof ]] be a basis of V h . du (i)
h can then be expressed as follows:

du (i) h = N dof k=1 du (i) S,k N k = N du (i) S , with du (i) S = du (i) S,k k∈[[1;N dof ]]
.

(3.11)

Then, writing Equation (3.10) for dv h = N j yields:

∀j ∈ [[1; N dof ]], N dof k=1 du (i) S,k a (i) S (N k , N j ) = l (i) S (N j )
Previous equation can be rewritten as a linear system:

M (i) S du (i) S = b (i) S , (3.12) with ∀i ∈ [[0; i max ]],      M (i) S = a (i) S (N j , N k ) (j,k)∈[[1;N dof ]] 2 b (i) S = l (i) S (N j ) j∈[[1;N dof ]]
.

(3.13)

One should note that at each step of the iteration process the problem is changed. Indeed both

M (i)
S and b

(i)
S depend on the current estimation u (i) . This Gauss-Newton scheme is thus computationally expensive since matrix M (i) S has to be reassembled and refactored at each step i. In order to reduce computational costs, ∇g • φ u (i) is often replaced by ∇f in (3.8) [Besnard et al. 2006] and as a result in (3.12):

M S du (i) S = b (i) S,f , (3.14) with ∀i ∈ [[0; i max ]],      M S = a S,f (N j , N k ) (j,k)∈[[1;N dof ]] 2 = S N T ∇f ∇f T N b (i) S,f = l (i) S,f (N j ) j∈[[1;N dof ]] = S f -g • φ u (i) N T ∇f , (3.15)
which is the usual Quasi-Gauss-Newton (QGN) scheme implemented in DIC.

Remark By taking u

(0) ∈ V h , ∀i ∈ [[0; i max ]], u (i) ∈ V h .
Remark Because M S involves two shape functions in the integral, it can be seen as a (weighted) mass matrix.

Remark [Passieux & Bouclier 2019] explain that the condition to be met so that the iteration process still converges after ∇f was substituted to ∇g • φ u (i) is:

-b (i) S T M S -1 b (i) S,f < 0 (3.16)
which is a more general case than the small deformation condition ( ∇u 1) usually raised to justify this QGN implementation [Neggers et al. 2016, Passieux et al. 2018]. Typically, it is shown that the QGN scheme converges for large rotations below 90 • , provided that an appropriate step size is used. This QGN implementation is particularly interesting from an algorithmic viewpoint, as the correlation matrix M S is assembled and factorised once and for all in the beginning of the optimisation process. Also, it should be noted that vector ∇f T N , once computed, may be used for both matrix M S assembly and vector b (i) S,f , and only the residual f -g • φ u (i) needs to be updated.

The considered domain S has not been clearly defined yet. Its definition actually depends on whether a local or global approach is implemented.

We defined the ROI as the planar surface over which the displacement measurement is performed. It should be noted that this definition may raise some semantic issues with the community of local DIC. In this community, the ROI is indeed understood as the domain in the reference image corresponding to the planar surface (see for instance [Pan 2009, Introduction]), not as the planar surface itself. In local approaches, this is a very subtle difference. In global approaches, where a model of the surface is available, it actually makes sense to distinguish these two concepts.

Local approaches

In the first case, S stands for a particular subset of pixels: the Zone Of Interest (ZOI). It is included in a wider area, the domain corresponding to the ROI in the reference image (see Figure 3.2). In that case, shape functions are the same for all ZOIs, and their support is included in the associated ZOI. Hence, the iteration process (3.14) is performed for each ZOI independently. ZOIs can be contiguous, overlapping, separated [START_REF] Hild | Comparison of Local and Global Approaches to Digital Image Correlation[END_REF], they may even contain discontinuities [START_REF] Bourdin | [END_REF].

ZOIs are typically square groups of pixels, called subsets, from the reference picture f and shape functions are low-order polynomials. The first coefficients of the polynomial allow to retrieve the displacement field as well as the strain over a set of points corresponding to the ZOI centres.

As ZOIs are independent from one another, local approaches exhibit high performance through parallelisation and seed techniques [START_REF] Pan | Reliability-guided digital image correlation for image deformation measurement[END_REF]].

Global approaches

In the second case, S stands for the whole domain corresponding to the ROI in the reference image. It means that the iteration process has to be performed once and for all, allowing to recover the displacement field directly. Shape functions N j j∈ [[1;N dof ]] support can also be the whole domain S.

An important asset of global approaches compared to the local ones is the opportunity to enforce regularisation on the displacement field based on numerical and/or experimental measurements [Pierré et al. 2017] because shape functions may be the exact same ones as in numerical simulations. However, these methods require a model of the planar surface (e.g. a FE mesh).

Remark The distinction between local and global approaches is somewhat artificial as one could argue that global approaches encompass local ones when choosing the same discontinuous shape functions.

Including camera models

In this section, we assume that a model of the surface is available, so that we can adopt a global approach.

Until now, the domain S has been defined in the reference picture f and the displacement field u as well. This results in a displacement field which is expressed in terms of pixels over a domain parameterised with pixel coordinates. When considering a perfect pinhole model (see Section 2.1.1), it is straightforward to convert pixels in terms of physical length (simple scaling in each direction of the planar surface). However adopting such an approach may result in a displacement misevaluation because of geometrical flaws in the camera optical system (see Section 2.1.2). To circumvent these issues, the camera model P from Chapter 2 may be introduced in the grey level conservation equation as follows [START_REF] Pierré | [END_REF]:

∀X ∈ Ω, f • P (X, p) = g • P (X + U (X), p),
(3.17)

where U = (U V ) T is the sought displacement field defined over Ω which denotes the model corresponding to the ROI (S = P (Ω)). Note that, unlike u from Equation (3.1) which is expressed in pixels, the dimension of U is a physical length (e.g. expressed in meters). In the following, to reduce the amount of notation, we will simply denote P (X) instead of P (X, p).

With such notations the variational formulation associated with Equation (3.17) is:

Find U ∈ L 2 (Ω) 2 , U = arg min U F (U ), F (U ) = Ω f •P (X)-g•P (X+U (X)) 2 dX. (3.18)
To derive the QGN scheme associated with Equation (3.18), the same steps as in Sections 3.1 and 3.2 can be followed and it all comes down to replace f (resp. g) by f • P (resp. g • P ). Noting that, for instance, ∇ (f • P ) = ∇P ∇f • P , the Gauss-Newton scheme writes:

M (i) Ω dU (i) = b (i) Ω , (3.19) with ∀i ∈ [[0; i max ]],        M (i) Ω = Ω N T ∇P ∇g • P ∇P ∇g • P T • φ U (i) N b (i) Ω = Ω f • P -g • P • φ U (i) N T ∇P ∇g • P • φ U (i) , (3.20) where φ U (i) : Ω ⊂ W → W X → X + U (i) (X) . (3.21) Eventually, using ∇P ∇f • P instead of ∇P ∇g • P • φ U (i) (that is substituting ∇ (f • P ) to (∇ (g • P )) • φ U (i)
) yields the associated QGN scheme:

M Ω dU (i) = b (i) Ω,f , (3.22) with ∀i ∈ [[0; i max ]],      M Ω = Ω N T ∇P ∇f • P ∇P ∇f • P T N b (i) Ω,f = Ω f • P -g • P • φ U (i) N T ∇P ∇f • P . (3.23)
This last QGN scheme is used in Section 4.2.

Remark

We would like to stress that there is something unsatisfying with Equations (3.2) and (3.18). It seems indeed equivalent to integrate the residual squared of Equation (3.1) over S or to integrate the one of Equation (3.17) over Ω, which is not. This issue will be addressed later on in Chapter 5, especially in Section 5.2.

Numerical implementation

In previous developments, only mathematical and algorithmic aspects have been discussed. However to solve the DIC problem, we will implement the algorithms on a computer. This requires to address two different topics. First, strategies to numerically compute integrals should be defined. Also, images are acquired and stored in a digital way, meaning that, images correspond to pixel arrays and associated grey level values are defined at pixel centres only. Hence, evaluating, for instance, g • φ u for subpixel values of the displacement is not straightforward. This requires interpolation procedures.

Integration schemes

The algorithms (3.14) and (3.22) will be numerically implemented. That is the reason why we need a quadrature rule in order to respectively estimate S • dx and Ω • dX.

Pixel centres as integration points Classically, a sum over the pixels in the reference image f is used to evaluate S • dx from Equation (3.14). This exhibits an interesting feature: raw data from f may be used. However, the projection of an element of the FE mesh (when considering FE-DIC for instance) in f may not correspond exactly to a set of pixels and thus may intersect some (see Figure 3.3). Hence even constant functions cannot be integrated exactly (meaning that the area of an element can be miscalculated) and it implicitly assumes that the projection of a line of the physical space in the image is a line, which is not necessarily the case as soon as distortions occur in the optical system (see Figure 3.

3).

Remark There is also an unsatisfactory dissymmetry between f and g with this quadrature rule: every data from f is a raw grey level, while almost every data from g is an interpolated one (see next section).

FE integration scheme

We consider that a FE model of the surface is available and that N collects the shape functions associated to the mesh. In this case, the aforementioned issues may be tackled by relying on the approach described in Section 3.3. By denoting N e the set of mesh elements and

X k e k∈[[1,N k e ]]
the integration points of element e, we can approximate the integral of a function h defined over Ω:

Ω h(X) dX ≈ e∈Ne N k e k=1 ω k e h(X k e ), (3.24)
where ω k e is the weight associated to the integration point X k e . There are various ways to define the integrations points and their respective weights [Pierré 2016]. One efficient way is Figure 3.3: With the classical method, the computed area, in red, is not exact. Extracted and modified from [START_REF] Pierré | [END_REF]]. described in Figure 3.4. For a triangular element, the idea is to look for the number of pixels intersected by each edge of the triangle projection in the image f . Then, these numbers of pixels along each direction are used to define the integration points in the isoparametric domain together with the elementary area that they represent. This allows to assign the weight values ω k e based on the area of each elementary area in the physical space. Compared to classical pixel integration, this integration scheme reduces biases in displacement measurements [START_REF] Pierré | [END_REF]. Besides some specific features discussed later in Chapter 7 on a practical application, this is the integration scheme used in this work (in particular in Section 4.2).

Remark This process is easily extended to quadrangles [Pierré 2016]. It may also be applied to other global approaches.

Subpixel image evaluation and image gradient

There are two final aspects that have not been tackled yet to allow for a numerical implementation of the discussed algorithms. First, the need for a method allowing to evaluate images at subpixel positions. Even though f may be evaluated only at pixel centres, which is straightforward, g may not if a subpixel accuracy on the displacement is aimed for. Second, the image gradient should be evaluated as well. To address these two points, a bicubic B-Spline interpolation of images is used (the reader can refer to [Bornert et al. 2009, Schreier et al. 2000, Zhou et al. 2015] for more precisions on the effects of interpolation methods).

Correlation criteria

Practically, there is another issue that we will discuss in detail in this manuscript, particularly in the SDIC case. It is linked to differences in lighting between reference state image f and deformed state one g which result in measurement errors. These differences may be due to Projection in the image

Isoparametric coordinate system

Figure 3.4: Integration points definition for a triangular mesh element. They are positioned so that a density of one point per pixel is reached. The weights ω k e are assigned so that they correspond to each elementary area defined in the top schematic drawing. Hence, their sum equals the triangle area. Extracted from [Pierré 2016]. limitation of the imaging hardware (e.g. scanning electron microscope) or surface degradation [Tong 2005]. To account for these effects, corrections of f and g may be introduced. This results in a correlation criterion (or functional) different from a simple L 2 -norm of the residual. Many different DIC criteria exist and have been evaluated. In [Tong 2005], it is shown that the most robust and reliable criterion is the Zero-mean Normalised Sum of Squared Difference (ZNSSD), with:

C ZN SSD = 1 |S| S   f -f ∆f - g • φ u -g • φ u ∆ g • φ u   2 , (3.25)
where f and g • φ u denote, respectively, the mean value of f and g • φ u , ∆f and ∆ g • φ u denote f and g • φ u standard deviation, and |S| stands for the ROI projection area. It is also expressed in [Tong 2005] the relation between this criterion and the Zero-mean Normalised Cross-Correlation (ZNCC) defined as:

C ZN CC = 1 |S| S f -f ∆f   g • φ u -g • φ u ∆ g • φ u   . (3.26)
Since we can write:

C ZN SSD = 2(1 -C ZN CC ). (3.27)
It is worth pointing out that, for whole-pixel translations, computing C ZN CC or the pattern auto-correlation function at the corresponding point would yield the same results. Hence, it is consistent to assess pattern quality by means of criteria based on the auto-correlation function for the ZNCC (or ZNSSD) cost function, as done in Section 4.2.

Remark

The ZNCC and ZNSSD correlation criteria have been presented for Equation (3.2) (without camera model) but similar criteria may be written for Equation (3.18) (with camera model).

In Section 3.5.2 is discussed the use of such criteria in a SDIC context.

Extension to Stereo Digital Image Correlation

Pieces of this section are adapted from [Fouque et al. 2021a, Introduction].

The main features of SDIC compared to DIC are its ability to perform a three-dimensional displacement field measurement on a possibly non-planar surface Ω. To do so, at least two cameras are needed (see Figure 3.5). In this section, we present two different frameworks used in SDIC. First, the subset-based and then the FE-based ones are introduced. This allows to show the main differences between these two approaches as well as their common characteristics. One of their common traits is the need to calibrate the stereo rig. This can be done (or refined) at various steps but, for the sake of simplicity, we will assume it is the first step of each framework.

For each approach, its extension to multiview pictures of the ROI is introduced (see Figure 3.10 that is presented in the next paragraphs). By multiview pictures we mean that for each state more than two pictures are available, allowing to observe a large part of the structure surface. Note that, with the approaches presented in this section, there are as many reference state images as deformed state ones, which naturally leads to setups where there are as many cameras as there are pictures at each state.

Subset-based Stereo Digital Image Correlation

Case of a single camera pair Considering only two cameras, in Subset-based SDIC [START_REF] Jones | [END_REF],Lucas & Kanade 1981,Luo et al. 1993,Sutton et al. 2009[START_REF] Synnergren | Per Synnergren and Mikael Sjödahl. A stereoscopic digital speckle photography system for 3-D displacement field measurements[END_REF], a master camera is chosen (the left one in Figure 3.6) and the subsets are defined in the reference state picture of this camera (f L ). Basically, each subset of the left image is sought in the right image thanks to standard 2D DIC technique. The knowledge of the positions of the cameras with respect to each other allows then to obtain the position of each subset centre in the world reference frame (triangulation). A 3D representation of the reference state of the whole ROI as a point cloud is then generated. We denote it S. Regarding the deformed state, the procedure is the same, except that the subsets remain defined in the reference state left image f L . Thus additional 2D DIC procedures are used to obtain the position of each subset in both left and right deformed state images g L and g R . Then, thanks again to triangulation, the position of the subset centres in the deformed state S are retrieved. Finally, the displacement U is obtained as a difference between the two point clouds: U = S -S. Let us stress that, consequently, U is not a field strictly speaking, and is the result of several (five) nested optimisation problems.

R c R' c M 1 M 2 m 1 m' 1 m' 2 O c O' c

Remark

In the previous paragraph and in Figure 3.6, a specific strategy regarding temporal matching has been presented but others exist (see Figure 3.7).

Extension to multiview Stereo Digital Image Correlation

The process associated with multiview subset-based SDIC is mapped in Figure 3.10a. Cameras are paired and, for each reference state image pair (f i,L , f i,R ), the point cloud S i corresponding to the projection of the ROI in the image f i,L is retrieved according to the same steps as in the previous section. The different point clouds are then stitched together to obtain the whole shape S of the ROI [START_REF] Solav | [END_REF]]. To perform displacement measurements, the deformed state images (g i,L , g i,R ) are considered and similarly allow, after stitching the results S i , to recover the whole deformed shape S of the ROI. Again, the displacement is evaluated as U = S -S.

Finite-Element-based Stereo Digital Image Correlation

In contrast, in the case of FE-SDIC (see Figure 3.10b), no so-called master camera is defined. Instead, pictures are interrogated at pixel positions corresponding to physical points defined on a FE mesh of the ROI [Réthoré et al. 2013] (see Figure 3.10b), as described in Section 3.4.1. Note that, apart from specific points that will be discussed later on, the extension of FE-SDIC from only two cameras to multiview setups is straightforward. We will thus present the multiview case directly. To do so, we consider a set of N c cameras, each of which took a reference state image of the ROI Ω. The associated set of reference state images and camera models are respectively denoted I 0 c c and (P c ) c . Classically, global SDIC frameworks rely on the steps described in Figure 3.8. As already explained, cameras are first calibrated. Then the relative position of the cameras with respect to the model is searched, which corresponds to the extrinsic calibration. And, before displacements can eventually be measured, a shape correction is introduced to account for small defects of the true geometry with respect to the model. Calibration procedure In the extrinsics calibration and shape measurement steps based on reference state images (I 0 c ) c , the position of each physical point is adjusted so that the different grey levels associated to the projection of this point in each picture match. The best camera parameters p ext c c

t 0 + Δt f L f R g L g R U=S'-S S' S t 0
(extrinsics only) and shape S minimising the discrepancies over all reference state image pairs is sought. This allows to ensure another kind of grey level conservation assuming that the grey level associated to a physical point should be the same for all cameras:

∀c ∈ [[1, N c ]], ∀i ∈ [[1, c -1]], ∀X ∈ Ω, I 0 c • P c X + S(X), p ext c = I 0 i • P i X + S(X), p ext i
(3.28) Based on the same reasoning as the one presented in Section 3.1, the calibration procedure relies on a variational formulation. The associated functional reads:

Find S , p ext ∈ L 2 (Ω) 3 × R 6Nc , S p ext = arg min (Sp ext ) F S, p ext F S, p ext = Nc c=1 c-1 i=1 Ω I 0 i • P i (X + S(X), p ext ) -I 0 c • P c (X + S(X), p ext ) 2 dX, (3.29)
where, for instance, P c is the camera model associated to camera c. p ext is the vector collecting extrinsic parameters of all cameras. Note that we generalised the definition of the projector P to reduce the amount of notation in what follows (we omitted the intrinsic parameters which were calibrated in a prior step and P c is a function of the vector of all extrinsics p ext instead of simply p ext c ).

Remark Note that, as the camera rig is assumed to be calibrated, the relative position of each camera with respect to c = 0 is known (see Chapter 2). Thus it is also possible to formulate a problem where p ext = p ext 0 , that is the rig moves as a whole in the extrinsic calibration step. This is what is done in Chapter 7 (with a virtual rig).

However, Equation (3.29) happens to be ill-posed because of the solution non-uniqueness. For instance, grey level discrepancies resulting from rigid-body motion of the camera rig may be accounted for either by a shape correction field if the space in which the shape is sought contains rigid-body modes, or by a correction in the extrinsic parameters. As it is difficult to ensure that shape space does not encompass (even infinitesimal) rigid-body modes, extrinsics are usually calibrated before shape is measured. Similarly to what is presented in Section 3.1 a Gauss-Newton algorithm may be implemented to minimise Equation (3.29).

• Extrinsics calibration: The Gauss-Newton algorithm corresponding to the minimisation of Equation (3.29) with respect to p ext is:

H (i),ext dp (i),ext = b (i),ext , (3.30) with                                    H (i),ext = Nc c=1 c-1 i=1 Ω ∇ p P i ∇I 0 i • P i -∇ p P c ∇I 0 c • P c × ∇ p P i ∇I 0 i • P i -∇ p P c ∇I 0 c • P c T X + S(X), p (i),ext dX b (i),ext = - Nc c=1 c-1 i=1 Ω I 0 i • P i -I 0 c • P c × ∇ p P i ∇I 0 i • P i -∇ p P c ∇I 0 c • P c X + S(X), p (i),ext dX , (3.31)
where p (i+1),ext = p (i),ext + dp (i),ext and, for instance,

∇ p P c = ∂P c ∂p ext .
• Shape measurement: The Gauss-Newton algorithm corresponding to shape measurement is:

H (i),shape dS (i) = b (i),shape , (3.32) with                                        H (i),shape = Nc c=1 c-1 i=1 Ω N T (X) ∇P i ∇I 0 i • P i -∇P c ∇I 0 c • P c × ∇P i ∇I 0 i • P i -∇P c ∇I 0 c • P c T X + S (i) (X), p ext ×N (X)dX b (i),shape = - Nc c=1 c-1 i=1 Ω N T (X) I 0 i • P i -I 0 c • P c × ∇P i ∇I 0 i • P i -∇P c ∇I 0 c • P c X + S (i) (X), p ext dX ,
(3.33) where S (i+1) = S (i) + dS (i) and S (i) (X) = N (X) S (i) .

At this point, shape measurement remains extremely ill-posed, as the mesh may slide on the object (see Figure 3.9). The cause is that the texture and the mesh are not attached to one another. To regularise the problem, a strategy consists in constraining the shape measurement along the normal to the model. Practically, the normal at each point of the FE mesh is measured by averaging the normals of neighbouring elements. Then a matrix R shape is assembled such that the normal component s(X) = S(X), n(X) of the shape at the nodes of the mesh, denoted s, gives S = R s. This allows to write a Gauss-Newton minimisation scheme for the reduced variable s:

R T shape H (i),shape R shape ds (i) = R T shape b (i),shape . (3.34)
It simply consists of projecting System (3.32) onto the normal directions to the surface. Remark Iterations between these two steps (extrinsics calibration and shape measurement) may be performed so as to obtain both faster and better convergence [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF], Colantonio 2020].

Remark When cameras c and i do not observe the exact same zones in the ROI (e.g. strongly curved surface), a weight may be assigned to the associated residual so that only the area shared by both cameras is considered in the functional (3.29) [START_REF] Chapelier | [END_REF]. However this remains unclear how to consistently define this weight and details are provided regarding this matter in Section 6.1.

Remark Another strategy to perform the shape and extrinsic measurement relies on the definition of an 'intrinsic texture' [Dufour et al. 2015b]. Details regarding this concept are provided in Chapter 5 and Section 6.1.1.

Remark Defining normals at mesh nodes is not straightforward, in particular in the case of sharp edges. A practical way to proceed, more satisfying than averaging the normals of neighbouring elements, is proposed in Section 7.1.4.

Displacement measurement

Then, the displacement measurement step is performed by adjusting the new position of physical points so that the grey levels in the deformed state pictures (I c ) c associated to this point match the reference state ones (I 0 c ) c . However, and unlike the previous step, the minimisation of grey level discrepancy is performed on a camera by camera basis only (i.e. registration between images I c and I 0 c ):

∀c ∈ [[1, N c ]], ∀X ∈ Ω, I 0 c • P c (X) = I c • P c (X + U (X)) . (3.35)
Again, a variational formulation is used to solve this ill-posed problem:

Find U ∈ L 2 (Ω) 3 , U = arg min U F (U ) F (U ) = Nc c=1 Ω I 0 c • P c (X) -I c • P c (X + U (X)) 2 dX.
(3.36)

And the associated QGN scheme writes:

H dU (i) = b (i) , (3.37) with            H = Nc c=1 Ω N T ∇P c ∇I 0 c • P c ∇P c ∇I 0 c • P c T N b (i) = Nc c=1 Ω I 0 c • P c -I c • P c • φ U (i) N T ∇P c ∇I 0 c • P c . (3.38)
Note that H and b (i) are very close to simple additions of the respective two-dimensional cor-

relation terms M Ω and b (i) Ω,f of Equation (3.

22). Remark

The grey level conservation is not enforced over all camera pairs any longer (compared to the calibration step). This allows to reduce computational costs and ill-posedness. Also, as the point of view remains the same a lot of issues associated with light are simply avoided. However it comes at the cost that stereo correspondence may be lost

(i.e. ∃(c, i) ∈ [[1, N c ]] 2 , ∃X ∈ Ω, I c • P c (X) = I i • P i (X)
even in the least-squares sense). This also implies that the same area should be followed, and thus visible, by a given camera at all times and it may become an issue when considering large displacements, especially in a multiscale context for near-field cameras.

When considering many time steps (N t > 2 time steps for instance) and associated pictures

I t c (c,t)∈[[1,Nc]]×[[0,Nt-1]] ,
the displacement measurement procedure is exactly the same except that I c from Equation (3.37) is replaced by I t c . This allows to measure a displacement with respect to the reference state (captured by I 0 c c ) at all times.

Correlation criteria

In subset-based SDIC, the measurement process relies on the correlation of independent small patches in pictures. Correlation criteria introduced in Equations (3.25) and (3.26) are thus equivalent to the application of corrections to each subset. This makes this approach very robust to local illumination differences in images which may be caused by surface normal orientation variations for instance.

In FE-SDIC, and generally in global approaches, the minimisation is performed over the whole ROI Ω. Thus it makes it harder to account for these effects, even when relying on a ZNSSD correlation criterion for instance. However, it is possible to mitigate the impact of illumination changes by distorting the functional and applying a ZNSSD criterion on an element basis for instance [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF],Colantonio 2020]. Note that in this case, a continuous formulation becomes harder, if not impossible, to write. Another approach (Brightness and Contrast Corrections (BCC)) consists in using two low-order polynomials (a c , b c ) and replacing I 0 c by a c I 0 c + b c [Charbal et al. 2020]. Remark Note that these three methods, namely ZNSSD in subset-based SDIC, ZNSSD on an element basis in FE-SDIC and BCC, are a brightness and contrast adjustment. The difference being that when using a ZNSSD criterion, one implicitly assumes that the subset (in the case of subset-based SDIC) or the element (in the case of FE-SDIC) stands for a larger surface than the representative elementary surface. Hence coefficients corresponding to a c and b c in BCC are directly computed from the mean and standard deviation of the considered image portion.

Discussions

Subset-based SDIC is a very robust and fast measurement method as it relies on well-studied two-dimensional image registration techniques. This is made possible by the independent subsets which (a) do not interfere with one another in the case of mismatch and (b) allow high performance through parallelisation. Yet, the measurement process is not focused on the displacement field itself, as it is computed in a post-processing step using outputs of several minimisation problems. In comparison, FE-SDIC (and global approaches in general) is more expensive computationally speaking and requires a model prior to any measurements. However, in a test-simulation dialogue perspective, as the same kinematic basis as simulations may be chosen, the comparison between data from tests and simulations is made much easier. Also, the displacement field U is directly the solution of a minimisation problem, which allows to make use of a priori knowledge of the sought displacement field [START_REF] Bouclier | [END_REF], Passieux et al. 2018,Pierré et al. 2017,Réthoré et al. 2009,Rouwane et al. 2021a]. With such data, a wide variety of powerful algorithms allows to adjust physical parameters associated with numerical models (e.g. FEMU, IDIC) [Hild & Roux 2006, Lecompte et al. 2007, Passieux et al. 2015a[START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF]], densely measure shapes [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF], Dufour et al. 2015b, Etievant et al. 2020] or even identify mode shapes during dynamic tests [Berny et al. 2018b, Passieux et al. 2018].

Stereo Digital Image Correlation and large-scale applications

Performing multiscale measurements using SDIC is not straightforward [Passieux et al. 2015a, Serra et al. 2017]. The formulation (3.36) does not provide a satisfactory way to merge results from pictures taken at various scales. Currently, only data from the most resolved pictures of a zone are considered to perform the associated displacement field measurements [Passieux et al. 2015a, Pierré 2016].

A major issue that we identified with these SDIC methods is that the objective to provide algorithms as efficient as possible might have taken precedence over more essential topics, especially concerning global approaches. As shown above in this chapter, both subset-based and global SDIC displacement measurement algorithms are based on their respective twodimensional counterparts. In other words, assumptions which were justifiable in the case of two-dimensional DIC became questionable. For instance, it is relevant in case of a planar surface, as done in Equations (3.1) and (3.17), to assume that the grey level of a physical point remains constant in pictures taken at different levels of deformation (as long as the texture remains attached to the surface and does not "fade" because of high strain levels and as long as the mechanical deformation does not change the orientation of the local normal vector with respect to the lighting too much). The normal at each point of the ROI can indeed be taken equal to Z c at all times, resulting in the same amount of light energy falling on each point of the surface, at all times. Conversely, if the object geometry is complex, the grey level value may strongly depend on the position on the object as shown in Figure 3.11.

Another example is the same weight assigned in Equation (3.17) for all point of the surface, which is consistent for a planar object located on a plane of normal Z c . However, as shown in Figure 3.12, the same level of confidence cannot be assigned to all unit physical surface patch. The confidence should depend, among others, on the surface normal n orientation.

Eventually, points of the ROI are assumed to be visible, or at least remain visible by a same camera all the time. Of course, fulfilling these assumptions in the case of SDIC is a complex task. It is even more complex in the case of multiview pictures and let us not speak of large structural tests. For these reasons, use cases of SDIC are seriously limited. To clarify this point, we will distinguish the case of small deformations (displacements, rotations and strains are assumed to be small) and the small-strain or large-deformation one (only strains are small). The matching of surface areas that were not visible in previous pictures from the exact same camera is not possible in current SDIC frameworks, thus it is limited to a certain class of large deformations. Also, effects associated with light and surface pixel sampling (which are not accounted for) are strongly related to the local normal vector to the surface [START_REF] Delaunoy | Photometric bundle adjustment for dense multi-view 3d modeling[END_REF],Tsiminaki et al. 2019], hence the current framework is rather limited to small-strain contexts. Overall, this results in the use of SDIC only in small-deformation cases. We believe this is related to the fact that grey level conservation equations (3.1) and (3.17) can be seen as a way to inverse a forward problem consisting in warping a flat picture with a given displacement field. Yet, in the Figure 3.11: Depending on the light position and the surface normal vector, an object of uniform color may not appear uniform in a picture, as illustrated by this sphere. Extracted from [Wehrwein 2015].

Figure 3.12: The sampling performed by a camera of a surface is strongly related to the orientation of the normal vector n. current state-of-the-art of SDIC, there is no so-called 'forward problem' allowing to generate virtual images at various deformation states under given lighting conditions of structural tests integrated in a framework (yet methods to generate virtual images to compare measures to a so-called 'ground truth' do exist [Balcaen et al. 2017, Garcia et al. 2013, Lava et al. 2020, Rohe & Jones 2021]). Thus, we cannot compare these predictions with the observations in actual pictures. As explained by [Tarantola 2005], "the comparison of the predicted outcome and the observed outcome allows us to ameliorate the theory".

Deposited patterns

As suggested by the different Gauss-Newton algorithms ((3.12), (3.14), (3.19), (3.22), (3.30), (3.32) and (3.37)) measurement accuracy in DIC greatly relies on the image gradient which is included in both the matrix that should be inversed in the left-hand side (tensor product with itself) and in the vector of the right-hand side. Note that we already mentioned that, in the strong sense, only the component along the image gradient of the displacement field may be recovered. Hence, the image gradient should exhibit properties allowing to perform measurements as reliable as possible. For that, in the context of mechanical tests, patterns are deposited on measurement surfaces, often by the mean of a paint spray.

As discussed later on, a lot of criteria exist to assess the quality of a pattern [Fouque et al. 2021c], however they are often adapted to a single scale or a discrete set of scales. Hence, before developing numerical methods to extend the use of SDIC to larger scales (which we do in Chapters 5 and 6), the possibility to rely on patterns which could be used in a multiscale context should be assessed. We thus present in the next chapter (namely Chapter 4) multiscale approaches in DIC. Then the use of a fractal pattern for multiscale measurements is introduced.

In the remainder of the manuscript, unless otherwise stated, we consider that we are using a FE-SDIC framework, or at least a global approach.

Chapter 4

Multiscale approaches in Digital Image Correlation

The multiscale aspects in DIC are introduced. They encompass two different but strongly intertwined topics. The first one concerns the displacement measurement initialisation, as multigrid schemes are often used. These schemes rely on a coarse-graining process which generates low-resolution pictures, at various scales, from an initial one. It shows the benefit to both reduce computational costs associated with one iteration of the optimisation algorithms, and smooth the DIC functional to avoid convergence to local minima.

The second topic is related to multiscale measurements. In this context, cameras with different levels of resolution are used to observe the ROI. The pattern deposited on the surface should thus allow measurements as accurate as possible at all scales, and facilitate the initialisation at a scale as coarse as possible. In both cases, it comes down to define a pattern with multiscale properties. However, this is not straightforward as shown by the literature review. The method that we developed to achieve this goal is presented. It stands for one of the contributions made during our PhD thesis [Fouque et al. 2021c].

Contents

Pyramidal approach for initialisation

The algorithms used in DIC have been presented in Chapter 3. They rely on the premise that a satisfactory initial displacement guess is available, so as to benefit from the local convergence of the used minimisation schemes. Practically, it may be hard to directly find a good initialisation for the displacement field and, to circumvent this issue, a pyramidal approach is often used.

Coarse-graining process

Basically, the idea of coarse graining is to aggregate pixels 4 by 4 by averaging their associated grey levels (subsampling), in order to produce a coarser image with a smaller number of pixels. This process can be repeated N s times and produces an N s -level multigrid or pyramidal scheme (see Figure 4.1). A DIC problem is solved at each level by a top-down approach, resulting in a coarse-to-fine search strategy [Hild et al. 2002, Réthoré et al. 2007]. This strategy exhibits many advantages [Burt 1984]. First, iteration costs of QGN algorithms are drastically reduced at the coarsest scales which will then provide a good initialisation for the finest ones. Then, most iterations are performed at the coarsest levels which are the cheapest ones computationally speaking. Eventually, a much greater robustness is achieved when relying on such a pyramidal scheme. The reason for this last asset is the low-pass filtering provided by the grey level averaging that allows to avoid undesirable convergence to local minima. It should be noted that the amount of available data decreases as the level in the multigrid scheme increases (at the level s, the number of pixels is divided by a factor 4 s compared to the initial pictures). At the same time, the number of unknowns (degrees of freedom) remains the same. As a result, the DIC problem becomes more and more ill-posed. Yet, it is very convenient to keep the same degrees of freedom from one level of the pyramid to another, as the expression of the sought field also remains the same. To prevent the ill-posedness to grow with the level, regularisation methods may be useful [START_REF] Luettgen | [END_REF].

Remark

The use of such methods does not limit to pyramidal approaches. They are often implemented when the strong regularisation allowed by the choice of the initial subspace V h (see Section 3.2) is not sufficient.

Regularisation methods in Digital Image Correlation

One may split regularisation methods into two classes, as evoked in [Dufour et al. 2015a]. An instance of the first class has already been presented in Section 3.2, that is a strong regularisation. In this case, the field is sought in a smaller space (e.g. V h ) embedded in the first one [START_REF] Chapelier | [END_REF],Colantonio et al. 2020,Etievant et al. 2020,Passieux et al. 2018,Réthoré et al. 2009, Serra et al. 2017]. The second class may be referred to as weak regularisation methods. It consists in "twisting" the variational formulation by adding a term to the functional. This term penalises the distance to an a priori knowledge on the solution. Yet, at least in some cases, it may be seen as a restriction of the space to a compact set in some topology [Benning & Burger 2018]. This approach is also widely known as Tikhonov regularisation [Benning & Burger 2018, Dufour et al. 2015a].

Strong regularisation

As mentioned above, we have already presented the outline of a strong regularisation when considering the displacement field u in V h ⊂ (L 2 (S)) 2 in Section 3.2. But this process can be repeated if the problem remains ill-posed as done in Section 3.5.2 concerning the shape measurement step for instance. In the following, we give two other examples of such strong regularisations.

In [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF], the shape measurement problem remained extremely ill-posed even when relying on measurements only along the normal to the FE mesh. As a remedy, a subspace of the shape functions associated to the FE mesh is chosen. This subspace is composed of Non-Uniform Rational B-Spline (NURBS) functions which smoothly vary in space, making the problem much less ill-posed. A pyramidal approach is adopted, and the strong regularisation of the shape measurement problem, based on the NURBS functions, relies on a multilevel strategy as well. That is, the number of degrees of freedom decreases as the level in the pyramid increases (see Figure 4.2). The interesting feature of this being that shape measurement is regularised, but more importantly, it is possible to compute the shape correction at the node of the FE mesh exactly.

However, this approach requires an analysis-suitable Computer-Aided Design (CAD) model, and such a model may not be available. Another similar strategy is to immerse the FE mesh in a morphing box [START_REF] Chapelier | [END_REF]. This box may be refined, following a pyramidal approach in the same way as previously presented, but this methodology applies to any mesh. One should also note that both approaches may also be applied, besides shape measurements, to regularise displacement measurements. Weak regularisation In the case of weak regularisation, the number of degrees of freedom of the DIC problem remains the same, instead, as explain earlier, a distance between the sought field and a point in the subspace, symbolising a prior knowledge, is penalised. The functional to be minimised writes:

F reg = F + w R F R R, (4.1)
where F is one of the functionals of Chapter 3, R is the regularisation functional, F and R are typical values taken respectively by F and R, and w R is a regularisation parameter. Many different kinds of regularisation functionals exist. We introduce only a few of them below in the two-dimensional case, for displacement measurements, but they may be generalised to three dimensions, and other quantities of interest. Matrix M S and Vector b

(i) S,f -w R F R Γ ν.
Depending on the chosen R, we will give the corresponding Γ and ν:

• Distance to a prior displacement field u 0 :

R(u) = 1 2 S u -u 0 2 2 , (4.2) yields:    Γ = S N T N ν = u (i) -u 0 (4.3)
Remark This regularisation is often used for material property identification.

• Gradient L 2 -norm squared: It corresponds to the elastic energy minimisation (in the case of plane-stress state and if Poisson ratio equals -1, ∆u = 0, see Appendix A):

R(u) = 1 2 S ∇u 2 F = 1 2 S ∇u T : ∇u, (4.4)
leads to following values for Γ and ν:

   Γ = S G T G ν = u (i) , (4.5)
where G u is the representation of ∇u consistent with the choice of shape functions N (see Section 3.2). Note that ∇u is simply ∇u where columns are appended to one another as a vector.

Remark As

S ∇u T : ∇u = u T S G T G u = - S u T ∆u, L = S G T G is sometimes
referred to as the discrete Laplacian operator.

Remark This gradient-based regularisation functional is often used even in plane-strain state and even though the Poisson's ratio differs from -1.

• Laplacian L 2 -norm squared: In the context of shape measurement, the gradient-based regularisation tends to flatten the solution too much [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF][START_REF] Horn | [END_REF]]. Another choice for the regularisation may be a functional based on the Laplacian [START_REF] Horn | [END_REF]]. Physically, it corresponds to the bending energy of a plate:

R(u) = 1 2 S ∆u 2 2 , (4.6)
which is often approximated by R(u) 1 2 u T L T L u [Rouwane et al. 2021b] for the practical reason that shape functions N are often polynomials of degree 1 [Mendoza et al. 2019]. This results in: Γ = L T L ν = u (i) .

(4.7)

Remark A strong regularisation may be 'weakened' simply by defining R as the distance between the identified field and its projection on the subspace associated with the considered strong regularisation.

Remark Of course, several regularisation functionals may be used at once. For instance in [Leclerc et al. 2011,Mendoza et al. 2019,Rouwane et al. 2021a], an elastic regularisation is used in the bulk while a Laplacian-based one is used at the non-free domain boundaries.

Fractal pattern for Multiscale Digital Image Correlation

A significant amount of the work presented in this section has been published in [Fouque et al. 2021c].

As explained earlier, one of the most important elements in Digital Image Correlation (DIC) is the pattern deposited on test sample surfaces, as measurement accuracy depends strongly on specific features of this pattern [START_REF] Crammond | [END_REF], Lecompte et al. 2006, Pan et al. 2008[START_REF] Yaofeng | Study of optimal subset size in digital image correlation of speckle pattern images[END_REF]. Depending on the application, suitable speckles can be generated on specimen surfaces using a wide variety of experimental methods. For instance, for large scale applications, an airbrush or marker pen may be useful [LePage et al. 2017, Lionello & Cristofolini 2014] while, for small scale applications, the focused ion beam (FIB) technique or spin coating can be used [Kammers & Daly 2011, Wang et al. 2012, Winiarski et al. 2012]. However, for many of these methods, the operator's experience may greatly influence the measurement accuracy obtained [START_REF] Pan | Digital image correlation for surface deformation measurement: historical developments, recent advances and future goals[END_REF]]. Concurrently, test standardisation and robust measurement methods are needed in the industrial context. This need raises demand for methods that are able to generate patterns with high reproducibility [Chen et al. 2015,Zhang et al. 2018].

Consequently, research has sought to define relevant pattern quality criteria. Some researchers have investigated the effect of speckle characteristics in detail [Lecompte et al. 2006, Reu 2014a,Reu 2014b,Reu 2015a,Reu 2015b,Reu 2015c], leading to the definition of an optimal speckle size of 3 to 5 pixels [Lionello & Cristofolini 2014]. We wish to point out here that this optimality criterion relates only to certain pattern generation techniques (e.g. airbrush and marker pen). For example, it relies on the assumption that a typical speckle size can be defined.

Increasing interest is concomitantly being shown in full-field measurement techniques that are able to operate at various scales [Bomarito et al. 2018,Mathew et al. 2018,Passieux et al. 2015a, Serra et al. 2017, Tanaka et al. 2011, Wang et al. 2008]. These techniques dramatically reduce both computational time and hardware requirements. However, in these conditions, it becomes questionable to define an optimal speckle size in pixels (since cameras with different magnification levels are considered). Hence, in some works, zones with different speckle sizes have been defined [Pierré et al. 2017] (see Figure 4.3). However, zones with smaller speckle sizes may appear uniformly grey for far-field cameras or at least generate very noisy results as the speckle size may be much less than 3 pixels. Figure 4.3: Plate of the VERTEX project [Serra et al. 2017]. Two different speckle sizes have been defined. Most of the speckle pattern is adapted to far-field cameras which observe the whole plate. For near-field cameras, observing crack propagation, the speckle size is much smaller which results in a better measurement accuracy. Some other works have adopted a different approach and assessed pattern quality on the basis of global parameters [Bomarito et al. 2017, Bossuyt 2013, Pan et al. 2010]. To date, patterns described as multiscale have been designed for only two different scales [Bomarito et al. 2018, Mathew et al. 2018] and it is therefore more appropriate to call them biscale patterns (see Figure 4.4). Unfortunately, this kind of pattern does not completely fulfil industrial needs. There is little chance that only two scales would be enough and that these two particular scales would be known precisely before the start of experiments. In this context, a truly multiscale pattern would prove particularly useful. Also, for DIC displacement measurements, its multiscale features would make it easier to use initialisation processes based on a pyramidal (or multigrid) approach [Anandan 1989, Hild et al. 2002, Pierré et al. 2017] (see Section 4.1).

In this section, for the sake of simplicity, neither the process of image capture nor associated problems, such as scene lighting, will be investigated. These topics will be treated carefully later on (see Chapter 5). We will focus rather on the generation of a scale-free pattern for texturing objects. This work thus appears as a necessary prerequisite to the development of numerical methods in DIC for large scale applications in a multiscale context, that we present in the remainder of this manuscript.

Pattern quality criteria

We have seen in Section 3.4.3 that the C ZN SSD criterion (see Equation (3.25)) is widely used in DIC. Also, a link has been established with the C ZN CC one (see Equations (3.26) and (3.27)). Eventually, we pointed out that, for whole-pixel translations, computing C ZN CC or the pattern auto-correlation function at the corresponding point would yield the same results. For this reason, we investigate pattern quality criteria based on the pattern auto-correlation function. A supplementary criterion, the Mean Intensity Gradient (MIG) is also introduced. [Bossuyt 2013] suggested quality criteria based on pattern auto-correlation. From physical considerations on properties that an ideal pattern should have (Sensitivity and Robustness), he deduced the corresponding features of the pattern auto-correlation function. In the following, the interesting features of auto-correlation are reviewed.

Criteria based on auto-correlation function

• Main auto-correlation peak sharpness radius: First, the pattern sensitivity needs to be assessed. For that purpose, [Triconnet et al. 2009] defined the main auto-correlation peak sharpness radii as the principal axes of the ellipse formed by the intersection of the osculating elliptic paraboloid in (0, 0) and the zero-height plane. The half-sum of these radii was then used as a quality criterion. [Bossuyt 2013] assumes that this criterion is closely related to the pattern displacement sensitivity. It may appear logical that, the smaller the radius, the greater the change in the auto-correlation function for a given subpixel translation.

• Watershed radius: Then, a good pattern should allow initialisation as far as possible from the displacement to be measured. To this end, [Bossuyt 2013] introduced a criterion based on the primary auto-correlation peak width to take the pattern robustness into account.

The broader the peak the further the initialisation may be from the actual displacement and still avoid undesirable convergence towards local minima of the cost function. This criterion is the watershed radius and it represents the radius of a circle with the same area as the projection of the main auto-correlation peak on a plane. The interested reader is invited to refer to [Bomarito et al. 2017] for more information.

• Secondary auto-correlation peak height: Later on, [Bomarito et al. 2017] suggested taking possible mismatches or convergence in local minima into account by considering displacement fields associated with a stretch deformation. As these displacement fields can lead to primary correlation peaks that are actually lower than secondary correlation peaks in a ZNCC criterion (see Figure 4.5), [Bomarito et al. 2017] introduced a metric based on the height of the secondary auto-correlation peak. et al. 2017]. Some works followed this path [Grédiac et al. 2019] and it led to a method called LSA. However, in a general DIC context, according to [Bomarito et al. 2017], these parameters alone cannot be considered as proper pattern quality criteria because of the non-uniqueness of the identified displacement (up to a translation of a whole number of squares). These patterns cannot be used in a multiscale DIC framework either as the dynamic range of checkerboard pattern pictures taken by far-field cameras is drastically reduced.

For these reasons, both MIG-based and auto-correlation-function-based criteria will be considered here to generate a suitable multiscale pattern for DIC, as described below.

Multiscale pattern generation from auto-correlation function

Based on the auto-correlation criteria mentioned above, [Bomarito et al. 2017] developed an optimisation metric and a framework to generate an optimal pattern. A Boolean parameterisation of the pattern was considered for this, i.e., each pixel could be given the value 0 or N gl -1 by the optimisation algorithm (where N gl is the number of quantisation levels). For each pattern, the auto-correlation function was computed in order to evaluate the metric. Thanks to this first step, [Bomarito et al. 2018] were able to create a method to generate an optimal biscale pattern in which two pixel sizes were defined, one for each camera magnification. However, the generalisation of this process to more than two different scales does not appear straightforward, particularly when a reasonable dynamic range must be obtained at all scales.

In this work, an alternative approach is suggested. We do not parameterise a pattern and then optimise the associated auto-correlation function quality criteria, nor do we define a set of magnification scales. Instead, a suitable auto-correlation function with appropriate multiscale properties is taken as an input. Then the associated pattern is generated directly. This section presents the method used to generate a fractal speckle pattern with respect to a priori criteria. The very same method could be employed to define patterns with user-defined auto-correlation functions.

Auto-correlation function and Fourier Transform

Let us consider a sampled and quantised pattern:

h : x ∈ [[0; 2n + 1]] 2 → h(x) ∈ [[0; N gl -1]], (4.8)
where N gl is the number of quantisation levels (256 levels with 8 bits), 2n + 2 the number of pixels in both dimensions and x the pixel centre. h auto-correlation function A h is defined as:

A h : τ ∈ [[-n; n + 1]] 2 → A h (τ ) ∈ [-1; 1], (4.9)
where τ = (τ x , τ y ) is the shift (or translation) vector. F T (h) denotes h Fourier transform. The inverse Fourier transform will be denoted by F T -1 (•). With such notations, we can compute the auto-correlation, up to a multiplicative constant, using the Wiener-Kinchin theorem [Scargle 1989] as follows:

A h = F T -1 (|F T (h)| 2 ), (4.10)
where |F T (h)| 2 is called h power spectrum. The pattern Fourier Transform modulus can easily be derived from Equation (4.10):

|F T (h)| = F T (A h ) 1 2 . (4.11)
At this point, it may be noted from Equation (4.11) that the amplitude of the pattern Fourier transform is known at each point of the frequency domain. However, the phase information cannot be recovered in the process. In the generation of fractal surfaces, a uniform [0; 2π] random phase ψ is often used at each point of the frequency domain [Feder 1988,Saupe 1988,Voss 1988].

For this reason, we decided to use such a random phase ψ:

F T (h) = F T (A h )
1 2 e jψ with j 2 = -1.

(4.12)

For a thorough analysis on the influence of phase on the pattern obtained, we refer to [de Castro et al. 2017]. Ultimately, the pattern can be generated from the auto-correlation function by applying the inverse Fourier transform and discarding the imaginary part in the previous equation:

h = Re F T -1 F T (A h ) 1 2 e jψ . (4.13)
In the remainder of this chapter, unless otherwise stated, h will stand for the redistributed (in [0; N gl ]) and quantised version of h of Equation (4.13).

Similarities can be found between this numerical way of generating patterns and the method defined in [Funamizu et al. 2013], which optically generates fractal laser speckles. Very similar algorithms for generating fractal surfaces can also be found in the corresponding literature [Feder 1988, Saupe 1988, Voss 1988]. Additionally, a random modulus can also be used in Equation (4.13). That is, F T (A h ) 1 2 can be multiplied in Equation (4.13) by a Gaussian random variable of zero mean (and possibly, unit variance) as in [Saupe 1988, Voss 1988].

In previous works [Bomarito et al. 2017,Bomarito et al. 2018,Bossuyt 2013,Stoilov et al. 2012], the pattern (or its Fourier Transform) was parameterised in order to optimise criteria based on the auto-correlation function. Hence, a computationally expensive inverse problem had to be solved to generate a pattern from these criteria. In contrast, this method is a direct way to obtain a pattern from its auto-correlation features thanks to Fourier Transform properties.

Pattern generation

• Auto-correlation function choice: First, it has to be noted that there are necessary conditions for a function A to be an auto-correlation function. A non-exhaustive list of required properties is given below:

• the image of A should be included in [-1; 1];

• A(0) = 1; • ∀τ , A(τ ) = A(-τ ).
According to these conditions and the criteria defined in Section 4.2.1, we choose a correlation function with circular symmetry, which can thus be captured by a unidimensional function. Some kind of power law (e.g. an nth root function) seems to be an interesting candidate in terms of primary peak sharpness radius, watershed radius and secondary auto-correlation peak height. It is indeed ideally sharp (infinite derivative at the origin) and has no secondary peak. In this work, we choose to use a power type distribution:

A(τ ) = 1 - τ n 2H , (4.14)
where H < 1/2 in order to keep an ideally sharp main auto-correlation peak and 0 < H in order to have A(0) = 1.

This choice for the auto-correlation function has a direct implication on the nature of the pattern that we will obtain. By considering a fractal surface (more specifically a self-affine surface), it is shown in [START_REF] Shepard | [END_REF]] that the auto-correlation function is given by the same kind of function as Equation (4.14), where H denotes the Hurst exponent and 0 < H < 1 [Feder 1988, Shepard et al. 1995]. Reciprocally, the power spectrum of a surface is given by the Fourier transform of its auto-correlation function Equation (4.10).

Hence, the surface associated with A exhibits the same power spectrum as that of a fractal surface. Thus, by using Equation (4.13) with A h = A we would obtain a fractal pattern of Hurst exponent H [Feder 1988, Saupe 1988, Voss 1988]. Such self-affine patterns show some kind of statistical scale invariance: magnifying space coordinates x and y by a factor r (x and y become respectively rx and ry) requires the grey level value to be scaled by a factor r H to remain statistically identical [START_REF] Shepard | [END_REF][START_REF] Voss | [END_REF]].

Concurrently, these surfaces are not stationary [START_REF] Shepard | [END_REF]. Nonstationary surfaces have auto-correlation lengths (usually defined as the distance required for the auto-correlation function to drop from 1 to 1/e 0.37) that depend on the profile or surface area considered. It is therefore irrelevant to associate a correlation length with such patterns, unless it is, for instance, the characteristic length of the profile or surface used to compute the auto-correlation (this characteristic length was introduced in Equation (4.14) via the parameter n, which ensures that the auto-correlation function belongs to [-1; 1]). We expect a scale-invariant, correlation-length-free pattern to exhibit interesting properties in the context of multiscale DIC. Let us finally point out to the interested reader that some of these multiscale properties for such patterns have already been evoked in [START_REF] Wagne | [END_REF]. Also, from [START_REF] Shepard | [END_REF] and Equation (4.13), it is possible to show that the Root Mean Square (RMS) grey level difference between consecutive pixels is proportional to n -H . This demonstrates, in this particular case, that the local gradient is linked to the main auto-correlation peak sharpness radius defined in Section 4.2.1. As H decreases, the sharpness radius decreases and the local gradient increases, which results in better measurement accuracy. Conversely, as H decreases, the slope of the auto-correlation function far from the origin also decreases, making it more difficult to converge far away from the optimum. For 0 < H < 1/2, the fractal field is said to be antipersistent. For 1/2 < H < 1 it is said to be persistent. In the latter case, this means that an increase of grey level value over an arbitrary distance along a profile of h is likely to continue to increase for a comparable distance. In contrast, in the former case, an increase of grey level value is expected to be followed by a decrease [Feder 1988]. Therefore, there is an interesting trade-off between measurement accuracy and robustness depending on the value that is chosen for H. Then, from A(τ ), with H = 1/22, the pattern h is generated using Equation (4.13), and its actual auto-correlation function A h is computed with Equation (4.10) after having subtracted h mean grey level value. In the following, h only stands for the pattern generated for H = 1/22. A h (0, τ y ), symbolised by the orange dashed curve in Figure 4.6, can be compared to the solid blue line of the desired auto-correlation A.

• MIG improvement method: Since the focus has been placed on the auto-correlation function, MIG improvements have not been investigated so far. To remedy this, Figure 4.7 shows the histogram of the grey level distribution for different patterns and the MIG value computed for each of them. More specifically, Figure 4.7a shows the grey level distribution of h (the pattern that we generated in Section 4.2.2), for which a classical bell shape can be observed. This kind of grey level distribution is not optimal in terms of MIG. Heuristically, a way of increasing the MIG value would be to change the shape of the grey level distribution so that a greater number of pixels reached extreme values.

To do this in practice, we make use of the inverse method [START_REF] Devroye | General principles in random variate generation[END_REF]]. Starting with the distribution plotted on Figure 4.7a, a normal cumulative distribution function:

F µ,σ : x → N gl 2 1 + erf x -µ √ 2σ , (4.15)
was applied to obtain the uniformly distributed pattern of Figure 4.7b. In all that follows, unless otherwise stated, µ is set to the grey level mean value of image h (before quantisation), σ is a parameter to be chosen, and the grey level standard deviation of image h (before quantisation) is denoted by σ gl . More precisely, it can be seen in Figure 4.7b that a nearly uniform distribution is obtained by applying F µ,σ to h (before quantisation) for σ = σ gl , and a significant increase in the MIG value, from 23.3 to 58 is achieved. For the sake of simplicity, the redistributed (in [0; N gl -1]) and quantised version of F µ,σ applied to h (before h has been quantised) will be denoted F µ,σ (h). The Python script allowing both h and F µ,σ (h) to be generated is included as supplementary material of [Fouque et al. 2021c] for better understanding.

So as to spread grey levels even more towards extreme values and improve the MIG of the pattern, we now set σ equal to 0.3σ gl when applying F µ,σ to the pattern h (this choice for the value of σ will be explained in Section 4.2.3). Thus, as σ < σ gl , extreme values are more represented than values in the middle of the distribution. It can be seen, by observing Figures 4.6 and 4.7c, that this value for σ gives an even better MIG value, and sharpens the main auto-correlation peak. At the same time -although the possibility existed of this not being the case -the auto-correlation function keeps the same global shape.

Remark Another path could be followed in order to obtain any grey level distribution function. From the uniform grey level distribution of Figure 4.7b and by making use of the inverse method [START_REF] Devroye | General principles in random variate generation[END_REF]] a second time, the desired grey level distribution could be generated from the pattern h as long as its inverse cumulative distribution function was known. However, in this case, it cannot be guaranteed that the pattern auto-correlation function would remain practically unchanged.

Pattern performance evaluation

This section evaluates the suitability of the developed patterns a posteriori. To do this, we compare one of the proposed patterns with two others. As in [Bomarito et al. 2017], a random black and white pattern is taken as a reference. In addition, some numerical experiments are performed on a Perlin noise pattern for a deeper analysis.

Multiscale properties -Coarse graining

As expected from previous discussions, the patterns generated exhibit interesting multiscale properties. Figure 4.8 illustrates this statement by comparing one of the patterns developed (F µ,0.3σ gl (h)) with two other patterns, namely a randomly generated one (where each pixel can be equal to either 0 or 255 with a 50% chance) and a Perlin noise based one with speckles of 3 pixels. For these three patterns, a 32×32 subsampling was computed from the original 16384 × 16384 pictures. This meant that the subsampling required a 9-step coarse graining process (much greater than usual) of the corresponding original picture. As a result, each pixel in the coarse image was the integer part of a mean over more than 260,000 pixels from the original. Unsurprisingly, it led to an almost uniform grey picture for the random and the Perlin noise patterns, as shown in Figures 4.8h still greater than 200 grey levels (see Figures 4.8g and 4.8j). A file is included as supplementary material of [Fouque et al. 2021c]. It contains a picture of each pattern and the corresponding grey level distribution along the coarse graining process described here. This suggests that prospects are good in a multiscale context where cameras with different levels of magnification operate simultaneously [Passieux et al. 2015a], or in the context of large deformations.

Monoscale displacement-field measurements

In order to further evaluate the performance of one of the generated patterns (F µ,0.3σ gl (h)), it is now compared with the random and the Perlin noise patterns for two different kinds of displacement fields. The first kind of displacement field addressed corresponds to subpixel translations and the second one to stretches. These displacement fields were chosen because they are exactly represented by the FE shape functions and do not require any specific care from an algorithmic standpoint (as opposed to large rotations for instance [Neggers et al. 2016,Passieux & Bouclier 2019]). This reduces meshdependent and algorithm-dependent effects as much as possible, the idea being that similar results may be expected for subset-based DIC.

From a reference state image f (i.e. one of the three patterns considered), we generated a deformed state image g. The process for this generation will be described more in detail in each corresponding subsection. Then, in order to account for image noise acquisition, images were assumed to be polluted by a Gaussian noise [START_REF] Blaysat | [END_REF]]. Finally, each picture was quantified and grey level values above 255 and below 0 were set to 255 and 0 respectively (i.e. N gl = 256).

Adding Gaussian noise of variance σ 2 noise to both images f and g is equivalent to corrupting only g with a Gaussian noise of variance 2σ 2 noise [Roux & Hild 2006]. Consequently, we only generated noise for g and set the σ noise value to 3 grey levels, which is a conservative approach regarding the 2 grey levels evaluated in [Roux & Hild 2006]. From a practical point of view, for each measurement, 3 different noises were drawn, and the 3 corresponding displacement fields were measured using a FE-DIC algorithm [Besnard et al. 2006, Passieux 2018, Sun et al. 2005]. For the initialisation, a coarse-graining approach with a decreasing gradient-based Tikhonov regularisation was used (see Section 4.1).

The discrepancy between the measured displacement u m and the imposed one u imp was evaluated at each node of the mesh for the 3 measurements corresponding to each noise. Two quantities of interest were extracted:

1. The systematic error expectancy The second one shows a randomly generated pattern where each pixel has a one-half chance of being equal to 0 or 255. The last column concerns a pattern based on a Perlin noise with a speckle size of approximately 3 pixels. For each column, the first figure is a picture of the pattern considered for n = 8191. The second picture is a zoom on a 100×100-pixel square of the 16384×16384 pattern. The third one is a subsampling of the pattern considered (each grey level value is equal to the integer part of the mean over the corresponding 512×512-pixel square). The last picture is the grey level distribution of the subsampling. Colour scale for pictures goes from 0 (black) to 255 (white). where N n denotes the number of nodes of the mesh, u m αk the measured displacement for node k and direction α ∈ {x, y}, u imp αk the corresponding imposed displacement and • the expectancy, in the sense of the mean over the different noises.

e u = 1 2N n Nn k=1 α∈{x,y} u m αk -u imp αk , ( 4
Finally, picture size was set to 256×256 pixels (i.e. n = 127) and an unstructured triangular mesh, containing 132 nodes and 222 elements, was generated (see Figure 4.9). The mesh was located in the centre of f so that high levels of strain (up to 25%) could be reached without having pixels initially contained in the mesh moved out of the picture.

• Subpixel translations: The generation of g from f for subpixel translations was achieved by a phase shift in the Fourier space [START_REF] Périé | [END_REF]. In our framework, 21 subpixel translations ranging from 0 to 1 pixel in the {x}-direction were imposed by steps of 0.05 pixel. At the same time, the {y}-component of the imposed displacement was kept equal to 0. The multigrid level was set to 1 for the initialisation step.

In order to choose the value for σ, the influence of this parameter was studied as shown in Figure 4.10. We can observe a smooth decrease of the curve representing the systematic error expectancy as σ decreases until the curve reaches a minimum for σ = 0.3σ gl . pattern and the Perlin noise pattern in terms of systematic and random error expectancies. In order to conveniently compare all three patterns, for each expectancy, a scalar value is derived from the 21 measurements as in [Bomarito et al. 2017]. These values are the mean of the systematic error expectancy absolute value and the mean of the random error expectancy over the 21 subpixel translations respectively. Thus, using the proposed pattern generation technique, an 18.9% improvement for random error and a 95.1% improvement for systematic error can be obtained with respect to the random pattern. Compared to the Perlin noise pattern, the random error of the proposed pattern increases by 1.52% while, at the same time, its systematic error decreases by 77.9%.

In what follows, only F µ,0.3σ gl (h), the random and the Perlin noise patterns are considered.

• Stretches: Regarding stretch deformation, we decided to generate image g from a reference state image f such that g represented a symmetric stretch along the {x}-direction. It was performed via an inversion of the transformation φ.

f (x) = g • φ(x) with φ(x) = (x + ε(x -n), y) T had to be inversed. It yielded g(x) = f • φ -1 (x) with φ -1 (x) = ( x-n 1+ε + n, ỹ) T
, where ε denotes the tensile strain and n has been defined in Section 4.2.2. Finally, g was generated by evaluating f at non-integer pixel positions thanks to a bivariate spline approximation. Since expected values for the displacement were much higher than for subpixel translations, we used a two-level multigrid initialisation. Values for ε ranged from 1% to 20% by steps of 1 percentage point.

Figure 4.11 shows both systematic and random errors for F µ,0.3σ gl (h) (orange diamonds), the random pattern (blue circles) and the Perlin noise pattern (green triangles) for the stretch displacement field. For each curve in Figures 4.11a and 4.11b, two different regimes are exhibited. A part where subpixel accuracy for systematic and random error expectancies is achieved (below 9% stretch for the random and the Perlin noise patterns, and 18% stretch for our pattern) and another part where the algorithm obviously converged to a local minimum since errors are about a pixel or above. Table 4.2 is composed in the same way as Table 4.1 except that averages for systematic and random error expectancies were not taken over the whole measurement range but only over the values where all the curves converged (i.e. for the stretch strain value ε ∈ {0.01, 0.02...0.09}). The pattern developed in this work is outperformed by 15.0% and 33.2% in terms of random and systematic error expectancies, respectively, by the Perlin noise pattern for average values of Table 4.2 and by 13.3% in terms of systematic error expectancy absolute value by the random pattern. However, it should be pointed out that a 35.4% improvement compared to a random pattern is obtained for the random error expectancy. But most importantly, the multiscale property and large watershed radius of the generated pattern result in its capability to converge to the global minimum for higher levels of deformation than either the random or Perlin ones; a 100% improvement in these levels is achieved by F µ,0.3σ gl (h) compared to these patterns. This robustness property should be of particular interest for large deformation applications [Chevalier & [START_REF] Chevalier | Tools for multiaxial validation of behavior laws chosen for modeling hyper-elasticity of rubber-like materials[END_REF]] and, most interestingly, in applications like metal forming, where intermediate images between a reference and significantly deformed images cannot be obtained [START_REF] Genovese | [END_REF]. 2.39 × 10 -4 2.03 × 10 -4

Mean of the random error expectancy (pixel)

1.37 × 10 -2 2.12 × 10 -2 1.19 × 10 -2 At this point, for each scale s and each pattern, only five scalar values assessing the corresponding measurement accuracy are derived:

• the mean over the 21 sub-pixel translations of the:

systematic error expectancy absolute value, -random error expectancy,

• the convergence robustness defined as the highest stretch strain value where the pattern managed to converge,

• the mean over the stretch strain values such that all three patterns managed to converge for the:

systematic error expectancy absolute value, -random error expectancy.

Remark For stretches, if a pattern at a given scale did not manage to converge (even for the smallest stretch strain value considered ε = 0.01), we did not include it.

Figure 4.12 shows corresponding results with respect to the scale s for the three patterns. For scales 5 and 6, the Perlin noise pattern could not converge, even for the smallest stretch strain value (ε = 0.01). This can be seen in Figure 4.12e, where the value plotted for the convergence robustness is 0. Strictly speaking, for these scales, the convergence robustness of the Perlin noise pattern is not 0 but is less than 0.01. As a result, and as shown in Figures 4.12c and 4.12d, the Perlin noise pattern was discarded for scales 5 and 6.

Relatively to the random and the Perlin noise patterns, and for the range of scales considered, the DIC errors associated with the pattern proposed in this work do not depend on the scale s. Except for the steep increase of the mean of the systematic error expectancy absolute value between scales 0 and 1 for translations (see Figure 4.12a), our pattern shows stable error values. Regarding the convergence robustness in Figure 4.12e, this output tends to improve as the scale increases for the proposed pattern while the opposite trend can be observed for the other two patterns.

Figures 4.12c Figure 4.12: Top, results for the sub-pixel translations; means are computed over the 21 translations considered. Middle, results for the stretches; means are computed only over the stretch strain values where all three patterns managed to converge (except for scales 5 and 6 where the Perlin noise pattern was discarded): {s = 0, ε ≤ 0.08}, {s = 1, ε ≤ 0.04},{s = 2, ε ≤ 0.05},{s = 3, ε ≤ 0.04},{s = 4, ε ≤ 0.04},{s = 5, ε ≤ 0.06},{ s = 6, ε ≤ 0.04}. Bottom, convergence robustness and legend. Left, mean of the systematic error expectancy absolute value. Right, mean of the random error expectancy. Results are presented for F µ,0.3σ gl (h), the random and the Perlin noise patterns. scale, these means were computed over the range of ε where all three patterns converged. Additionally, Table 4.3 gives the mean of the absolute value of the systematic error expectancy and the mean of the random error expectancy over the whole range of convergence of our pattern for stretches (i.e. {s = 0,

ε ≤ 0.17}, {s = 1, ε ≤ 0.16},{s = 2, ε ≤ 0.20},{s = 3, ε ≤ 0.20}, {s = 4, ε ≤ 0.25},{s = 5, ε ≤ 0.22},{ s = 6, ε ≤ 0.25}
). This shows that the errors for the highest strain rates are similar to those for the smallest ones. If this was not the case, the values of

Conclusion and outlook

In this chapter, the pyramidal approach on which global DIC relies as an initialisation method has been presented. Some of the regularisation functionals encountered in this context, and others, have also been introduced. It allowed us to make a link with multiscale measurements in DIC, as these two topics are closely related. There is indeed an analogy between the pyramid levels and the levels of magnification of the cameras. In both cases, the pattern deposited on test sample surfaces plays an important role in the measurement accuracy and in the range of scales at which measurements may be performed.

A method based on the literature on fractal surface generation has been proposed to generate patterns directly from a desired auto-correlation function. Unlike the methods used so far, this generation process does not require any optimisation loop. This way of generating patterns has offered the opportunity to clarify what a good pattern is in a multiscale DIC context thanks to a priori quality criteria. One of the main unprecedented ideas that has arisen from these discussions is that a correlation-length-free pattern, i.e. a pattern with no typical speckle size, should exhibit interesting multiscale properties. This new way of considering patterns enables a family of truly multiscale patterns (and not only biscale or n-scale patterns) to be defined, based on fractal (self-affine) statistical scale-invariance. Also, it helps overcome some limitations linked to defining speckle sizes when different levels of magnification are involved (usually, a speckle size of 3 to 5 pixels is recommended). At this point, a complementary process could be implemented based on image filtering [START_REF] Pan | Bias error reduction of digital image correlation using Gaussian prefiltering[END_REF], Zhou et al. 2015] with parameters set accordingly to each camera zoom level in order to increase measurement accuracy at each scale considered. In the present work, a way has also been found to use the inverse method [START_REF] Devroye | General principles in random variate generation[END_REF]] to redistribute grey levels for patterns that we generated in order to improve their MIG value. For the distribution function used and these patterns specifically, this process kept their autocorrelation function practically unchanged.

Numerical experiments have been performed for two different displacement fields in order to assess a posteriori the quality of one of the patterns generated compared to a random black and white pattern and a Perlin noise pattern. These experiments were first performed at a single scale, and then for multiple magnification up to a factor 64 (2 6 ). They corroborated the expectations formulated beforehand. The pattern considered showed DIC errors comparable to those of the other two patterns for the first scales, but a much greater robustness. Most importantly, it showed errors and robustness that were stable with respect to the scale whereas these two outputs became significantly degraded for the other two patterns as the scale increased. Making use of the framework developed here, one may think of performing a closedloop optimisation [Lavatelli et al. 2019]. By considering the auto-correlation function shape (with Equation (4.14) parameterised by H or even with a different parameterisation) and the grey level distribution as input parameters, it would be possible to fine-tune them and optimise error expectancies for specific displacement fields.

The practical question of the best way to deposit such a pattern on test sample surfaces is still open. At Institut Clément Ader (ICA), it has been printed on plates for the VIRTUal testing of aerOnautical compoSite structurEs (VIRTUOSE) project [START_REF] Serra | VIRTUOSE Project[END_REF]]. This will allow to assess the use of such a pattern on an experimental setup and in a multiscale SDIC context.

It is now established that multiscale patterns, able to operate at a wide range of scale, may be defined. The next step now consists in developing DIC formulations in a multiscale context, which is the object of the remaining of the manuscript.

Chapter 5

Photometric Digital Image Correlation (PhDIC) formulation

Machine vision should be based on a thorough understanding of image formation. [START_REF] Horn | [END_REF] We present a framework to fill the gap of an inexistent so-called 'forward problem' allowing to generate virtual images to compare actual ones with, in DIC. To do so, we introduce existing approaches from CV, in a condensed and comprehensive manner, and adapt them to bring novel and significant capabilities for full-field measurements in experimental solid mechanics. This results in the development of a general multiview formulation for SDIC, that allows, in particular, to account for surface curvature and lighting effects. More precisely, a test Digital Twin is defined, which makes it possible to compare a model with different observations (real images taken from different viewpoints). We believe that these developments, together with the associated justifications in Chapter 6, stand for the main contribution of the present manuscript. As explained earlier, serious limitations in SDIC arise from strong implicit assumptions in the current frameworks. Among others, this results in the need for different functionals to be minimised in the global frameworks (i.e. one functional for extrinsics calibration and shape measurements and another one for displacement measurements), but most importantly in the lack of physical understanding. Therefore, we believe a unified framework including a rendering model, mapping an 'intrinsic texture' (further details are provided regarding this expression in Section 5.1) to a grey level is needed in SDIC. This would allow to rely on an explicit physical modelling that could be further improved instead of on implicit assumptions. Another benefit from this approach is to take advantage of grey level variations associated with lighting effects and extract information from it, instead of arbitrary correcting images to obtain a better matching. As shown in Figure 5.1, assuming an infinitely distant light source and the heater 'intrinsic texture' to be uniform, we can obtain a strong prior knowledge on the considered shape.

Figure 5.1: White heater. This picture illustrates that some information regarding the shape of a surface may be extracted from its shading. Courtesy of Jean-Charles Passieux.

Interestingly, formulations including rendering models already exist in the CV Community. They are very similar to those used in DIC (probably because DIC emerged from CV in the 80's). Not only do they propose a physical model, but they come up with astonishing results. Among others, these CV formulations encompass multiview picture frameworks accounting for simple lighting effects [Birkbeck et al. 2006], frameworks accounting for spatial sampling [Goldlücke et al. 2014] or even frameworks allowing to retrieve shape, texture and camera poses without any prior knowledge [START_REF] Jancosek | Michal Jancosek and Tomás Pajdla. Multi-view reconstruction preserving weakly-supported surfaces[END_REF][START_REF] Mes | Meshroom open-source 3D Reconstruction Software[END_REF], Moulon et al. 2012]. Surprisingly, such recent works are almost never included as references in related DIC research works. Here, we modestly aim at introducing such formulations to the DIC community while showing the benefits it could yield for performing data assimilation from structural tests [Wang et al. 2015].

The idea is to extend the Global SDIC framework, in which a model of the structure is available (e.g. a FE mesh) while taking advantage of the strong regularisation allowed by the knowledge of the initial object texture, such as the one provided by the master reference state image in Subset-based SDIC. To achieve this goal, we propose in the present work a framework relying on a Textured Digital Twin of the ROI constituted thanks to a physical modelling of the scene (see Figure 7.1) and call this framework Photometric Stereo Digital Image Correlation (PhDIC) as it relies on a photometric error. This allows to elaborate a refined model of the considered test, the Digital Twin, encompassing information regarding the structure shape, surface reflectance properties and scene lighting. Note that, within this framework, data are not compared to one another anymore. Instead, the Digital Twin is used as a common thread to compare pictures with. In this sense, a parallel can be drawn with the work proposed in Digital Volume Correlation by [Leclerc et al. 2015] or in Virtual Image Correlation by [Réthoré & François 2014, Semin et al. 2011]. Note that this would be delicate, if not impossible, to implement when relying on local approaches. However, a slightly different path may be adopted as in [Lava et al. 2020], where virtual images of the model are generated. Then, local SDIC is run on these virtual images and results are compared with actual measurements. Another proposal of [Rohe & Jones 2021] is to generate synthetic images from image rendering tools (namely Blender) to assess subset-based DIC frameworks performance. In these two cases, one should keep in mind that biases are expected from the image generation process associated with the comparison of two different displacement fields.

To build the above-mentioned physical modelling of the scene, we propose to consider the camera as a sensor converting a surface power density into grey levels (in the same way as it is usually considered as a projector mapping a 3D point in the scene to a 2D point in the image). This offers the opportunity to provide a general framework and to show that further assumptions regarding surface response to lighting and light distribution are needed. Especially, simple models relevant in the DIC context can be introduced (e.g. the Lambertian -surface response -and distant point light source -light distribution -ones). This allows to build a general unified photometric functional. Then the Lambertian and distant point light source models are plugged into this functional. A particular care is also taken to review some important and delicate points discussed in CV in a condensed and comprehensible manner for their significance in the experimental mechanics community.

From the light source to the grey level value

Here below, we present the irradiance equation and the associated physical considerations. We first intend to put forward some challenges in computer graphics and rendering (among others). This general overview is then progressively simplified by explicitly making assumptions that are commonly implicitly adopted in the DIC community.

A general overview

As explained earlier, we wish here to describe a general unified framework for SDIC, based on a photometric error. By photometric error we mean an error thought as a distance between a prediction from a model and an observation in a picture. This kind of error can be based on the image irradiance equation [START_REF] Horn | [END_REF]]. Irradiance I (or E) is the amount of light falling on a surface (power per unit area) while radiance L is the amount of light radiated from a surface (power per unit area per unit solid angle). As discussed by [START_REF] Horn | [END_REF]], the slightly more intricate unit for radiance comes from the fact that a surface can emit different amounts of light depending on the emission directions. The image irradiance equation states that the radiance L e coming from a point on an object (what we want to model) is proportional to the irradiance I at the corresponding point in the image (what we observe).

To formalise this equation, we denote x the pixel coordinates in the image of the associated 3D point X in the world reference frame. In this case the image irradiance equation writes:

I(x) = αL e (X, r),
(5.1) r is the unit vector pointing from X to the optical centre of the camera (see Figure 5.2), and α (called the throughput [Cohen et al. 1993]) depends on the f-number and the angle formed by r and the optical axis of the camera. However, this dependency of α on the angle formed by r and the optical axis is usually neglected [START_REF] Horn | [END_REF]] so that α is considered constant. In the following, we will not distinguish the image irradiance from the corresponding grey level value, as we assume the camera sensor to provide a linear relationship between these two quantities. The multiplicative constant is thus included in α. We refer the reader to [Durou 2007, 1.2 Optique photographique] for a more detailed explanation. Note that Equation (5.1) is fundamentally different from the grey level conservation equation [START_REF] Horn | [END_REF],Lucas & Kanade 1981].

Remark This model assumes a homogeneous and transparent medium. Hence the distance from X to the viewpoint was omitted in Equation (5.1). The dependency with respect to light frequency was also omitted, as we usually deal with monochrome cameras in DIC.

Before choosing a simplified model for the radiance emitted by a surface under given lighting conditions (model L e (X, r)), we will first introduce a more general framework. In the general case, [START_REF] Horn | [END_REF]] explains that the radiance dL e emitted from a point depends on the amount of light falling on it, the irradiance dE f , as well as on the irradiance fraction which is reflected per unit solid angle. The radiance dL e also depends on the geometry and light position, as illustrated by specular reflections. Thus, we can locally parameterise the problem thanks to 4 degrees of freedom: 2 for the incident light direction (the incident polar and azimuth angle -respectively θ f and ϕ f ) which allow to define a unit incident vector n f (θ f , ϕ f ) and 2 for the emission direction (respectively θ e and ϕ e ) which allow to define a unit emission vector n e (θ e , ϕ e ). The definition of these angles with respect to the local normal vector n(X) and an arbitrary vector belonging to the tangent plane to the surface can be seen in Figure 5.2. The fraction of incident light coming from the direction (θ f , ϕ f ) reflected in the direction (θ e , ϕ e ) is usually denoted f (θ f , ϕ f , θ e , ϕ e ) (per unit solid angle) and is called the Bidirectional Reflectance-Distribution Function (BRDF). For the sake of simplicity, we omit, in f and in following developments, the space dependency of each quantity (with respect to X), even though we consider a spatially varying BRDF. Remark Since we assume that the only way out for incoming energy is to be reflected, and since the incident light should come from the outside of the surface, both θ e and θ f belong to [0, π/2], as shown in Figure 5.2. Thus, effects such as transmission and subsurface scattering are not accounted for. For an even more general concept than BRDF which is called Bidirectional Scattering-Surface Reflectance-Distribution Function (BSSRDF), we refer the interested reader to [Nicodemus et al. 1977].

θ f φ f (θ f ,φ ) n f f (θ e ,φ e ) n e dΩ f n ( ) X X (θ r ,φ ) n f r I ( ) x r( ) X
Hence, the radiance can be written as a function of the irradiance and the BRDF:

dL e (θ e , ϕ e ) = f (θ f , ϕ f , θ e , ϕ e )dE f (θ f , ϕ f ), (5.2) with dE f (θ f , ϕ f ) = L f (θ f , ϕ f ) n f (θ f , ϕ f ), n
Normal component of the radiance

dΩ f = L f (θ f , ϕ f ) cos θ f sin θ f dϕ f dθ f dΩ f , (5.3)
where L f (θ f , ϕ f ) denotes the incident radiance coming from the direction -n f (θ f , ϕ f ) and dΩ f the solid angle delimited by [θ f , θ f + dθ f ] and [ϕ f , ϕ f + dϕ f ]. Thus, the total radiance emitted by the surface in direction n e is given by the integral over all elementary contributions coming from every single direction:

L e (θ e , ϕ e ) = π/2 θ f =0 2π ϕ f =0 f (θ f , ϕ f , θ e , ϕ e )L f (θ f , ϕ f ) cos θ f sin θ f dϕ f dθ f .
(5.4)

Remark Equation (5.4) is called the reflection equation. In computer graphics, it is called the rendering equation and a second term in the right-hand side may be included. The latter is an outgoing radiance in the case where the surface emits light by itself, in addition to the reflection from incident light. We chose to discard this term as, in structural mechanics, materials generally do not act as light sources.

At this point, we can observe that the emitted radiance coming from (θ e , ϕ e ) depends on the contributions of all incident radiances L f (θ f , ϕ f ). Each of these incident radiances is in turn the solution of the same kind of equation as Equation (5.4) and so on. This is an infinitedimensional problem and further assumptions are needed to be able to model this radiance. The path that we follow is first to model the surface response to light, that is the BRDF f , and then to model the scene lighting, that is L f . Remark Note the difficulty to define an 'intrinsic' texture. The albedo ρ is defined as the ratio of emitted irradiance over incident irradiance and depends in general on the incident radiance distribution:

ρ(Ω f , Ω e , L f ) = Ω f Ωe f (θ f , ϕ f , θ e , ϕ e )L f (θ f , ϕ f ) cos θ f cos θ e dΩ f dΩ e Ω f L f (θ f , ϕ f ) cos θ f dΩ f , 0 ≤ ρ ≤ 1.
(5.5)

Lambertian model

Most of the time in SDIC, the grey level conservation equation is used. In Global SDIC, whether it be for shape measurement, where the grey level associated to a physical point is assumed to be the same in each camera, or for displacement measurement, where the grey level associated to a physical point is assumed to remain constant in time for a given camera, it relies on a Lambertian assumption. That is, the incoming light is assumed to be reflected with equal intensity L e in all directions (θ e , φ e ) (i.e. ∀(X, θ e , φ e ), L e (X, θ e , φ e ) = L e (X)). Obviously, this assumption is not correct if the light field does not meet some strong assumptions or if the motion/rotation of the object is significant for the displacement measurement step (e.g. if the surface orientation with respect to the light changes). Regarding the shape measurement one, the grey level conservation equation Equation (3.28) assumes that the throughput α is the same for all cameras, which may not be the case.

As explained earlier, some works tried to account for surface illumination changes and optical system differences in cameras. The first step toward this goal in DIC was to use a Zeromean Normalised Sum of Squared Differences cost function as a matching criterion between pictures [Tong 2005]. The same idea was introduced earlier by [START_REF] Faugeras | [END_REF]] in CV. Because a correction on the whole ROI was not always sufficient to explain higher residual values in some areas, local corrections were introduced in SDIC on a finite-element basis [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF] for instance, or thanks to low-order polynomials [Dufour et al. 2015b, Charbal et al. 2020] (see Section 3.4.3).

Instead of a Lambertian assumption, or unphysical corrections, we use here a Lambertian model. It follows that f (θ f , ϕ f , X, θ e , ϕ e ) = f (X) = ρ(X)/π [Nicodemus et al. 1977, Appendix C]. From Equation (5.4), we get:

L e (X) = ρ(X) π π/2 θ f =0 2π ϕ f =0 L f (θ f , ϕ f , X) cos θ f sin θ f dϕ f dθ f . (5.6)
Note that ρ is independent of the viewpoint and light direction in the Lambertian model. Thus, we could refer to it as the intrinsic texture, but we prefer a physical designation: the albedo.

From Equation (5.6), we can see that further assumptions are required to model the radiance L e . Naturally, the knowledge of how the scene is lit (L f ) is also an important thing to compute it.

Lighting models

Ambient lighting The grey level conservation equations (3.28) and (3.35) also rely on the assumption that the same amount of light falls from every single direction:

∀(X, θ f , ϕ f ), L f (X, θ f , ϕ f ) = L 0 f .
In CV, it is referred to as ambient lighting but often as a component of a more complex model. Yet, in the case of SDIC, it is assumed to be the only source of light. For this reason, we prefer the denomination diffuse lighting, not to be confused with the Lambertian (or diffuse) model. In this case, together with the Lambertian model, the image irradiance equation (5.1) writes:

I(x) = αL e (X) = αL 0 f ρ(X).
(5.7)

Remark When relying on a ZNSSD correlation criterion, previous equation is equivalent to grey level conservation equations and a (Zero-Mean Normalised) substitute image may be defined such that ∀X, Î(X) = αL 0 f ρ(X) = I(x). In the shape measurement step, if several pictures are at stake (e.g. I 0 1 and I 0 2 ), we obtain Equation (3.28), as

I 0 1 = Î = I 0 2 .
Remark Diffuse lighting is a very strong assumption which may be met only by restricting SDIC to operate on certain types of tests and/or with specific lighting equipments.

The distant point light source Infinitely distant point light sources offer the advantage that neither the surface power, nor the direction of the light depend on the position in the scene. These light sources are parameterised by r ∈ [[1, N r ]]. It is also assumed that indirect illumination (contribution of the radiance emitted by all other points in the scene to the incident radiance to a point) is negligible compared to the direct illumination from the light sources. In this case, L f does not depend on X and

L f (θ f , ϕ f ) = Nr r=1 Φ r δ(θ f -θ r )δ(ϕ f -ϕ r ) sin θ r [Horn 1986],
where Φ r denotes the irradiance and (θ r , ϕ r ) the direction associated to the rays of light source r (see again Figure 5.2). Finally, we get:

L e (X) = ρ(X) π Nr r=1 Φ r cos θ r = ρ(X) π Nr r=1 Φ r n f (θ r , ϕ r ), n(X) .
(5.8)

Hence, the image irradiance equation (5.1) becomes:

I(x) = αL e (X) = ρ(X) Nr r=1 l r n f (θ r , ϕ r ), n(X) , with l r = αΦ r π .
(5.9)

Remark From Equation (5.9), [Horn 1986, Lightness & Color] and [Woodham 1980] it is possible to identify ρ(X)l 1 with one single light (l 1 ) but at least three different non-coplanar lighting conditions (three non-coplanar vectors n f ). An object in the scene of which the albedo is known allows to evaluate l 1 and thus to retrieve ρ(X). It is important to stress that varying lighting conditions is essential to estimate the albedo [START_REF] Mélou | [END_REF], since papers in SDIC pretend that it is possible to recover a so-called intrinsic texture, independent of the experimental setup, with one single light position [Dufour et al. 2015b, Vitse et al. 2021].

Remark Thanks to Equation (5.9), we can make explicit the assumptions on which grey level conservation equation in DIC relies. In the displacement measurement step, it assumes that the scalar product n f (θ r , ϕ r ), n(X) remains constant over time. This results in constraining DIC to operate in setups where displacements and strains are rather small.

Of course, more sophisticated parametric models can be derived in the same way, accounting for an ambient lighting term, together with a Lambertian model and/or specular reflections [Birkbeck et al. 2006, Tianli Yu et al. 2004].

Now that both a substitute to the grey level conservation equation and a way to model the scene radiance have been introduced, the corresponding framework developed for SDIC can be presented.

The photometric functional

Based on the work by [Goldlücke et al. 2014], we consider a set of N c cameras. Each camera takes N t images of the ROI Ω (typically a surface of R 3 ). As previously, the image taken by the camera c ∈ [[1, N c ]] at time t ∈ [[0, N t -1]] is denoted: Finally, we need to introduce theoretically the backprojection operator: which denotes the inverse function of the restriction of P t c to Ωt c vis . As of now, we note that this operator is only introduced for the sake of consistency in our theoretical developments but will not be used in our implementation.

I t c : Π t c → R x → I t c (x) , ( 5
β t c : S t c → Ωt c vis ⊂ Ωt x → X , ( 5
Retrieving deformed states Ωt t on the sole basis of Equation (5.1) constitutes ill-posed problems, similar to the measurement of a displacement field based on the grey level conservation equation between a reference and a deformed state image. Equation (5.1) does not account indeed for spatial averaging (pixelation), that is we have a finite set of equations. On top of that, grey level values, which are the result of the sum of an integral of the irradiance over each photosensor together with noise, are quantised. The existence of a solution is thus not guaranteed. The usual way to deal with these issues is to reformulate the problem (whether it be shape measurement or displacement measurement) as a functional minimisation (see Chapter 3). Thus, a norm of the residual associated to Equation (5.1) is integrated and the functional is built up by adding the integrals of all images together. Finally, the configurations Ωt are sought in a smaller space (finite dimension because of the finite set of equations). Based on CV literature [START_REF] Faugeras | [END_REF], Soatto et al. 2003, Goldlücke et al. 2014], we pretend that the right place to compute these integrals are the silhouettes S t c . The idea is that the relevant quantum of information is the pixel. Let us assume that we can assign, at all times, the same level of confidence to two pixels coming from a same picture taken by camera c and, without loss of generality, to two pixels coming from different pictures. Hence, we may assign a unit weight in every image plane (if our last assumption is not fulfilled, we may assign the weight 1/σ 2 c to the image plane corresponding to camera c by assuming a Gaussian white uniform noise of variance σ c in images shot by c as in [Hild & Roux 2020]). The unit weight associated to the residual norm should thus be assigned in the image plane. We will discuss that later on (see Section 5.3). According to these considerations, the photometric functional, suitable for data assimilation as it relies on every piece of available information [START_REF] Tsiminaki | [END_REF] reads:

F = Nc c=1 Nt-1 t=0 S t c I t c (x) -α t c L e β t c (x), r t c 2 dx.
(5.12)

Here, data assimilation should be understood as a general method to take advantage of all available observations to evaluate quantities of interest. In our case (DIC), these quantities are typically displacements, camera parameters or the albedo for instance.

Remark

As will be shown in Chapter 7, this formulation is easily extended to setups where the number of cameras is not the same at all times (by substituting N t c to N c ) and where cameras extrinsic and/or intrinsic parameters may change.

Since S t c = P t c ( Ωt c vis ), we can then express this functional over the visible parts of the ROI Ωt c vis thanks to integrations by substitutions, in the same way as [Goldlücke et al. 2014]: (5.13) where

F = Nc c=1 Nt-1 t=0 Ωt c vis J t c • P t c (X) I t c • P t c (X) -α t c L e (X, r t c ) 2 dX,
J t c = det ∇P t c • β t c = det ∇β t c -1 = ∂β t c ∂u × ∂β t c ∂v -1 2 and x = (u, v) T . det ∇P t c
and det ∇β t c denote the area elements of the corresponding projection maps. Note that, as in [Goldlücke et al. 2014], we denoted here the differential (Jacobian matrix) in the same way as the gradient operator, to reduce the amount of notation.

In a general framework, computing J t c is complex and costly. For this reason, and by assuming a pinhole camera model without distortions, we give an analytical expression for this area element: (5.14) where X = O w M , O w and O t c,i respectively denote the world reference frame origin and the camera reference frame origin associated to the picture taken by camera i at time t, Z t c,i stands for the Z coordinate of point X in the camera reference frame associated to the picture taken by camera i at time t, and

J t i = f x t i f y t i (Z t c,i ) 2 O t c,i M 2 |Z t c,i | | n, r t i |,
r t i = -O t c,i M / O t c,i M 2
(see Figure 5.3 again). Be careful that, to avoid misunderstandings between indices and reference frames, we indexed cameras with i (instead of c) in previous expressions. Expression (5.14) can be found in the CV literature [START_REF] Soatto | Tales of shape and radiance in multiview stereo[END_REF][START_REF] Delaunoy | Photometric bundle adjustment for dense multi-view 3d modeling[END_REF]. However, no proof is given in the literature reviewed. For this reason, a detailed demonstration and physical interpretation of this equation are given in Appendix B.

We can then integrate over the whole observed region Ωt introducing a visibility function V t c : (5.15) with:

F = Nc c=1 Nt-1 t=0 Ωt J t c • P t c (X) V t c (X) I t c • P t c (X) -α t c L e (X, r t c ) 2 dX,
V t c (X) = 1 if P t c (X) ∈ S t c and β t c • P t c (X) = X 0 elsewhere .
(5.16)

The condition P t c (X) ∈ S t c makes sure that the projection of X lies in the image frame I t c while the condition β t c • P t c (X) = X ensures that the considered point is not hidden due to self-occlusion for instance.

Remark J t

c and V t c naturally appear when the residual is defined with unit weight in the images, no further assumptions are needed for this weighting scheme. This is a direct consequence of the adopted variational formulation [Goldlücke et al. 2014]. Further discussions regarding this matter are presented in Section 5.3.

At this point, we should note that the residual is computed over what is observed in the images. Thus, what is observed is a deformed (or uncalibrated) state Ωt in Equation (5.15). Yet, Ωt is one of the unknowns which should be described. Considering that a model is available, a simple way to do so is to introduce the discrepancy map D t . In the standard SDIC framework [Pierré et al. 2017], D t stands for either a shape correction field S (D 0 ), or a displacement field U defined on the configuration Ω standing for the nominal geometry (which may not be perfectly consistent with the solid mechanics formalism, see next remark). Usually, this discrepancy map belongs to the linear span of a set of chosen shape functions (e.g. FE shape functions [Pierré et al. 2017], splines [Dufour et al. 2015b]), but it should be stressed that no prior assumptions are needed regarding the discrepancy maps which still belong to an infinite-dimensional space at this point. We define:

φ D t : Ω ⊂ W → Ωt ⊂ W X → X + D t (X) ,
(5.17) in this case Ωt = φ D t (Ω), where Ω stands for the ROI nominal geometry (where shape correction S = D 0 is not accounted for).

Remark With such notations, the displacement field associated to the deformation of Ω between the reference state and the deformed state at time t is D t -D 0 . Ω corresponds indeed to an uncalibrated or nominal state before correcting the shape with S = D 0 .

Finally, we can express the functional over the nominal state Ω:

F = t,c Ω det ∇φ D t J t c • P t c V t c • φ D t (X) × I t c • P t c • φ D t (X) -α t c L e (φ D t (X), r t c • φ D t (X)) 2 dX,
(5.18)

Note that with this method, the functional used to identify a shape correction S (D 0 ) and a displacement field U (D t , t > 0) is exactly the same. This offers a consistent, unified formalism throughout the entire framework. Usually in Global SDIC, the functional associated to the extrinsics and shape measurement problem enforces in a weak way that the grey level associated to a physical point should be the same for all cameras (see Figure 3.10b). Thus, it consists of a sum over all camera pairs of the residual norm squared (see Equation (3.29)), while the functional associated to the displacement measurement is built as a sum of another kind of residual norm squared. This other residual is based on the conservation over time of the grey level associated to a given point on a camera by camera basis only [Pierré et al. 2017] (see Equation (3.36)).

Introducing the Lambertian reflectance and the distant point light source models from Equation (5.6) in Equation (5.18) allows finally to write a functional F taking into account a Lambertian model:

F = t,c Ω det ∇φ D t J t c • P t c V t c • φ D t (X) I t c • P t c (X + D t (X)) -ρt X + D t (X) Ns s=1 l t s,c n f (θ t s , ϕ t s ), n(X + D t (X)) 2 dX.
(5.19)

It also offers the possibility to make explicit an often implicit assumption in DIC. If the pattern deposited on the ROI is assumed to exactly follow the deformation of the specimen, and does not depend on the displacement or strain level, we can write:

∀X ∈ Ω, ∀t ∈ [[0, N t -1]], ρt (X + D t (X)) = ρ(X),
(5.20)

where ρ and ρt respectively stand for the albedos in the reference and deformed states. Eventually the PhDIC functional reads, in the case of a Lambertian BRDF, with distant point light sources:

F = t,c Ω ω t c (X) I t c • P t c (X + D t (X)) -ρ (X) Ns s=1 l t s,c n f (θ t s , ϕ t s ), n(X + D t (X)) 2 dX,
(5.21)

where

ω t c = det ∇φ D t J t c • P t c V t c • φ D t .
(5.22)

Discussions

As already evoked, it seems logical to compute the discrepancy between images and the model in the image domain, as the pixel stands for the elementary unit of information. Besides this heuristic justification, the weighting term J t c , that naturally arises when substituting Ωt for S t c between Equation (5.12) and Equation (5.15), is a key driver for defining F with unit weight in the images. It accounts for the foreshortening of the surface in input views (e.g. a surface is well described in a picture when viewed straight on). Hence, this term acts as an automatic regularisation of the variational problem while making the problem intrinsic, i.e. independent of the parameterisation chosen for the ROI. This is clearly established in CV [START_REF] Faugeras | [END_REF], Goldlücke et al. 2014, Soatto et al. 2003]. Also, the weighting term J t i would allow to define a consistent framework with multiple cameras with different resolutions and distances with respect to the specimen since it accounts for the spatial sampling of the surface (as shown in Appendix B with the distance Z t c,i and the focal lengths f x t i and f y t i in the case of a pinhole camera model). No arbitrary relative weights would be needed for more resolved or near-field cameras as this formulation intrinsically defines a weighting scheme through J t i . This is of particular interest in a multiscale context. The weighting term in Equation (5.21) det ∇φ

D J t c • P t c V t c
• φ D t alone explains why SDIC is restricted to a certain class of displacements and strains. Indeed, in most SDIC framework, both terms are assumed to be the same for all pictures and to remain constant over time. Thus, it means det ∇φ

D t ∼ 1 and J t c • P t c V t c • φ D t ∼ J t c • P t c V t c (i.e. X + D t (X) ∼ X or equivalently φ D t ∼ I where I denotes the identity function).
Regarding the last remark of Section 3.3, the soundness of integrating the residual squared of Equation (3.1) over S or to integrate the one of Equation (3.17) over Ω indifferently may now be discussed. In the light of Equation (5.14), we can justify it in a two-dimensional framework, and by assuming a pinhole camera without distortions together with O c M 2 ∼ Z c , as the two functionals Equation (3.2) and Equation (3.18) are the same up to a factor fxfy Z 2 c . As stated above, the use of a model enables us to define a functional based on the sum of actual errors, that is the difference between a model and an observation. Thus the uncertainty associated to the identified discrepancy map D t (standing equivalently for a displacement field U or a shape correction field S) in PhDIC would be reduced compared to the usual SDIC framework (see Section 6.2). This explains why some authors aimed at forming a substitute reference state image, in applications where the level of confidence in this reference is low for instance, by taking a mean over all available pictures [Berny et al. 2018a].

Finally, in the present work, and contrarily to the usual DIC framework, the camera model encompasses not only a projection model, but also a model to define the grey level value depending on the amount of energy received by the camera sensor. This requires the definition of a radiance model for the experimental setup encompassing both a light model and a Digital Twin of the structure, that is a test Digital Twin. The photometric functional is thus based on the comparison between a prediction (or a model) and an observation, instead of arbitrarily correcting observations to make them match, as done in classic DIC for instance. Thus, information is extracted from these discrepancies and it helps regularise the problem. This is the exact opposite from BCC, for instance, which makes the problem even more ill-posed by increasing the number of unknowns instead. In addition, grey scale residuals are considered the most objective way to probe the ability of a model to reproduce an experiment from images in DIC [Neggers et al. 2017] or to determine areas where the geometry should be refined [Kleinendorst et al. 2015]. It is therefore extremely important to analyse and model finely the different sources of grey level variation during an experiment. These variations cannot always be related to displacement alone, especially in stereo where light-geometry interaction effects can be substantial.

Chapter 6

Theoretical developments

The use of the PhDIC framework and of Digital Twins of experimental setups is justified from different viewpoints in this chapter. First, we show the improvements it could yield by establishing the link with usual frameworks. The PhDIC methodology is indeed a general framework as it encompasses usual ones and offers the opportunity to extend the scope of existing schemes. Then, we show by a sensitivity analysis, with weaker assumptions, that the measurement uncertainty is expected to be reduced compared to usual methods. Eventually, the PhDIC approach is justified a priori thanks to arguments stemming from estimation theory.

Link between PhDIC and usual multiview frameworks

The work presented in this section has been published in [Fouque et al. 2021b].

The aim of this section is to clearly establish the link between the PhDIC formulation and usual SDIC ones. We also show that this topic is closely related to the weighting schemes in functionals. In particular, this allows to give the appropriate weight to use in (more standard) SDIC.

State of the Art

We revisit the usual SDIC multiview frameworks with an increasing complexity. That is, we begin this state of the art by recalling the functionals presented in Chapter 3 on simple SDIC setups. For planar or near-planar test samples, every single point of the ROI can be seen at all times by each camera. Yet, when considering more complex geometries, visibility issues arise and the way to deal with these problems is presented. Eventually, the link between the photometric functional (5.12) and usual SDIC frameworks, when relying on the same assumptions is demonstrated. It is shown that the photometric functional encompasses usual frameworks.

Shape measurements

When considering planar (or near-planar) surfaces, it is possible to ensure that every single point of the ROI Ω remains visible by all cameras at all times. In this case, the functional associated to the shape measurement step reads [Pierré et al. 2017, Colantonio et al. 2020] (see Equation (3.29)):

F 1 D 0 , p c c = Nc c=1 c-1 i=1 Ω I 0 i • P i • φ D 0 (X) -I 0 c • P c • φ D 0 (X) 2 dX, (6.1)
where D 0 is the shape correction field and ∀X ∈ Ω, φ D 0 (X) = X + D 0 (X). There are two main drawbacks to this formulation. First, computational costs associated to the problem scale as N 2 c which is not ideal in a multi-camera setup. Second, the problem is extremely ill-posed as shown in Figure 3.9.

In order to cope with the aforementioned problem ill-posedness but also to account for the surface sampling performed by each camera sensor, a shape measurement functional based on a residual thought as the difference between a substitute image and actual images was built [Dufour et al. 2015b]. In this case, Equation (6.1) becomes:

F 1 D 0 , p c c , Î = Nc c=1 Ω I 0 c • P c • φ D 0 (X) -Î(X) 2 dX. (6.2)
The required assumptions to define a substitute image have been defined in Section 5.1. This formulation shows the benefit to tackle the two main issues identified previously. The problem scales as N c (like the displacement measurement one, see Equation (3.36)), and the optimisation procedure relies on an alternating optimisation (fixed-point algorithm) between the shape correction field, extrinsics and the substitute image. This last point (looking for shape and texture, and not only the shape) makes the formulation much less ill-posed than in [Pierré et al. 2017, Colantonio et al. 2020] as it discards the functional kernel directions such as local and global slidings (again see Figure 3.9). Then, authors explicitly dealt with visibility issues by relying on a weighting term based on visibility [START_REF] Chapelier | [END_REF]:

F 1 D 0 , p c c = Nc c=1 c-1 i=1 Ω V c (X)V i (X) I 0 i • P i • φ D 0 (X) -I 0 c • P c • φ D 0 (X) 2 dX, (6.3) where ∀c ∈ [[1, N c ]],
V c is the visibility function associated to camera c, simply defined as:

V c : Ω → {0, 1} X → 1 if X is visible by camera c 0 otherwise.
Note that the previous definition is less precise than Equation (5.16).

Displacement measurements

To engage in displacement measurements, we consider a deformed state of which the same N c cameras shoot the associated pictures I 1 c c . The displacement measurement functional is then given by [Pierré et al. 2017] (see Equation (3.36)):

F 2 (D 1 ) = Nc c=1 Ω I 1 c • P c • φ D 1 (X) -I 0 c • P c • φ D 0 (X) 2 dX. (6.4)
Visibility issues in the displacement measurement step were addressed by resorting on the assumption that a point visible by a camera in the reference state remains visible by this camera at all times, [Hild & Roux 2020]:

F 2 (D 1 ) = Nc c=1 Ω V c (X) I 1 c • P c • φ D 1 (X) -I 0 c • P c • φ D 0 (X) 2 dX. (6.5)
However, no clear justification for the weighting terms associated to the visibility is given, neither for Equation (6.3) [START_REF] Chapelier | [END_REF] nor for Equation (6.5) [Hild & Roux 2020]. Also, the general case, where the displacement field is such that a part of the structure may disappear from view, is not tackled.

Proposed functional

In this section, we wish to thoroughly establish weighting schemes for both shape and displacement measurements in Global SDIC frameworks. To this end, we establish the link between PhDIC and usual SDIC frameworks by explicitly adopting the same assumptions, that is a Lambertian reflectance model and an ambient lighting. In this case, the irradiance equation (5.1) takes the form of Equation (5.7) and a substitute image Î may be defined (see Section 5.1.3). The photometric functional then writes: (6.6) where ω t c is defined in Equation (5.22).

F D t t , p c c , Î = Nt-1 t=0 Nc c=1 Ω ω t c (X) I t c • P c • φ D t (X) -Î(X) 2 dX,
Remark ∀(X, c, t), ω t c (X) ≥ 0.

Remark Let us stress once again that with such notations, the displacement field associated to time t is D t -D 0 . Usually, after the shape measurement step, Ω is updated such that Ω = φ D 0 (Ω). However, if we want to be able to efficiently perform a minimisation with respect to all arguments of F (that is D t t , p c c and Î), constantly updating the integration domain of all integrals may not be the most effective minimisation strategy. Keeping that in mind, defining the displacement on the nominal shape Ω as Dt = D t -D 0 is a small price to pay.

Shape measurements

In this subsection, we develop F from Equation (6.6) so as to establish the link between this formulation and the usual SDIC shape measurement functionals (see Equation ( 6.1) and Equation (6.3)). For that, we consider N t = 1, that is only reference pictures I 0 c are available. Note that with such considerations, F is very close to F 1 in Equation ( 6.2):

F D 0 , p c c , Î = Nc c=1 Ω ω 0 c (X) I 0 c • P c • φ D 0 (X) -Î(X) 2 dX. (6.7)
As in [Dufour et al. 2015b], Î is obtained by minimising F , that is directly (least-squares solution of F minimisation):

∀X ∈ Ω, Nc c=1 ω 0 c (X) = 0, Î(X) = Nc c=1 ω 0 c I 0 c • P c • φ D 0 (X) Nc c =1 ω 0 c (X) . (6.8) Remark If ∃X ∈ Ω, Nc c=1 ω 0 c (X) = 0, then ∀c ∈ [[1, N c ]], ω 0 c (X) = 0 as ∀ (X, c) , ω 0 c (X) ≥ 0.
Practically, it means that the point X cannot be seen by any camera. Hence Î(X) can be set to any arbitrary real number without affecting the value of F . Thus, in what follows, we do not consider this case any longer.

To reduce the amount of notation, we denote f c = I 0 c • P c • φ D 0 and develop Equation (6.7) (Equation (6.6) with N t = 1):

F = Ω Nc c=1 ω 0 c f 2 c -2 Î Nc c=1 ω 0 c f c + Î2 Nc c=1 ω 0 c ,
using the identity (6.8), it follows:

F = Ω Nc c=1 ω 0 c f 2 c -Î2 Nc c=1 ω 0 c .
Making use a second time of (6.8):

F = Ω 1 k ω 0 k i c ω 0 i ω 0 c f 2 c - i c ω 0 i ω 0 c f i f c . As i c ω 0 i ω 0 c f 2 c = i c ω 0 i ω 0 c f 2 i : F = 1 2 Ω 1 k ω 0 k i c ω 0 i ω 0 c f 2 c -2 i c ω 0 i ω 0 c f i f c + i c ω 0 i ω 0 c f 2 i .
We can factor this expression:

F = 1 2 i c Ω ω 0 i ω 0 c k ω 0 k f 2 c -2f i f c + f 2 i ,
which can finally be rewritten:

F = c i<c Ω ω 0 c ω 0 i k ω 0 k (f c -f i ) 2 .
This last equation is very close to the functionals F 1 (6.1) and F 1 (6.3) used in a standard Global SDIC framework for the shape measurement step (see Section 6.1.1). These developments allow to establish a link between a weight assigned to each observation f c , namely ω 0 c , and the associated weight in the usual framework, which should be

ω 0 c ω 0 i k ω 0 k when comparing f c to f i .
In F 1 (6.1) [Pierré et al. 2017], it is (implicitly) assumed that:

ω 0 c = det ∇φ D 0 J c • P c V c • φ D 0 ∼ 1, because (near
) planar surfaces are considered. In this case, considering the correct weighting scheme does not change much the functional expression, as

ω 0 c ω 0 i k ω 0 k ∼ 1 N c
. Note that we just showed that, according to previous assumptions:

F 1 /N c = F 1 .
However, in F 1 (6.3) [START_REF] Chapelier | [END_REF], ω 0 c ∼ V c is assumed. Hence, when considering more complex geometries, the correct weight when comparing f c to f i should be

V c V i k V k , instead of V c V i .

Displacement measurements

Let us now consider the displacement measurement step. Before getting into the general case N t > 2, we establish the link between F (6.6) and usual frameworks considering only two time steps (N t = 2). It is done first by relying on a substitute image Î given by Equation (6.8) (only reference state images are used to build Î), and then by updating Î thanks to data provided by deformed state images.

Incremental displacement measurements

Since N t = 2, the functional F from Equation (6.6) writes as follows:

F D 0 , D 1 , p c c , Î = Nc c=1 Ω ω 0 c (X) I 0 c • P c • φ D 0 (X) -Î(X) 2 +ω 1 c (X) I 1 c • P c • φ D 1 (X) -Î(X) 2 dX = Nc c=1 Ω ω 0 c f c -Î 2 + ω 1 c g c -Î 2 , ( 6.9)
where

g c = I 1 c • P c • φ D 1 .
• Substitute image based on reference state images only: Here, as in [Dufour et al. 2015b] we keep on using the same substitute image based on Equation (6.8). In this case, the functional is minimised with respect to D 1 only and reads:

F = Nc c=1 Ω ω 0 c f c -Î 2 Constant=F 0 + Nc c=1 Ω ω 1 c g c -Î 2 F = F 0 + Nc c=1 Ω ω 1 c g 2 c -2 Îg c + Î2
We can then make use of Equation (6.8):

F = F 0 + Ω 1 Nc k=1 ω 0 k   Nc c=1 Nc i=1 ω 0 i ω 1 c g 2 c -2 Nc c=1 Nc i=1 ω 0 i ω 1 c g c f i + Nc c=1 ω 1 c Nc k=1 ω 0 k Nc i=1 Nc j=1 ω 0 j ω 0 i f i f j   = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k g 2 c -2g c f i + j ω 0 j f i f j k ω 0 k = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 -f 2 i + j ω 0 j f i f j k ω 0 k = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 -f 2 i + f i Î = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 -f i -Î 2 + Î( Î -f i ) = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 -f i -Î 2 + c Ω ω 1 c Î Î -i ω 0 i f i k ω 0 k .
The last sum equals 0 by definition of Î (6.8):

F = F 0 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 -f i -Î 2 = c Ω ω 0 c f c -Î 2 + i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 - c j Ω ω 0 c ω 1 j k ω 0 k f c -Î 2 = i c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 + c Ω ω 0 c (f c -Î) 2 1 - j ω 1 j k ω 0 k = c Ω ω 0 c ω 1 c k ω 0 k (g c -f c ) 2 + c i =c Ω ω 0 i ω 1 c k ω 0 k (g c -f i ) 2 + c Ω ω 0 c (f c -Î) 2 1 - j ω 1 j k ω 0 k . (6.10)
At this point, let us point out that relying on a substitute image for the displacement measurement step exhibits some interesting properties. First, it encompasses usual formulations of the shape measurement step (Equation (6.4) or Equation (6.5)) thanks to the first sum in Equation (6.10) that are, this time, complemented by every spatio-temporal cross-correlations g c -f i , with cameras i = c, which are usually not included. Also, as Equation (6.2) compared to Equation (6.1), Equation (6.6) scales linearly with the number of cameras, unlike Equation (6.10) which scales quadratically. Finally, note that the last term in Equation (6.10) may be neglected if for every point X of Ω,

j ω 1 j (X) k ω 0 k (X)
, that is if every point X is equally well observed in the pictures I 0 c c and in the pictures I 1 c c . Again, these developments establish a link between a weight associated to a given observation and the consistent weighting scheme that should be adopted in the usual framework.

In F 2 (6.4) [Pierré et al. 2017], it is assumed that ω 0 c ∼ ω 1 c ∼ 1 (which makes sense for near-planar surfaces only) and again, considering a consistent weighting scheme

ω 0 c ω 1 c k ω 0 k
when comparing f c to g c only scales F 2 by a constant factor 1 Nc . However, when introducing a visibility function such that ω 0 c ∼ ω 1 c ∼ V c as in F 2 (6.5) [Hild & Roux 2020], the consistent weight when comparing f c to g c should be

ω 0 c ω 1 c k ω 0 k ∼ V c k V k instead of V c .
• Substitute image updating: Here, we perform data assimilation in the sense that F from Equation (6.9) (Equation (6.6) with N t = 2) is minimised with respect to all arguments (that is, D 0 , D 1 , p c c and Î), unlike the previous subsection. For this reason, Î is updated by minimising F :

Î = Nc c=1 ω 0 c f c + ω 1 c g c Nc k=1 ω 0 k + ω 1 k . (6.11)
Hence, we can develop:

F = Ω Nc c=1 ω 0 c f 2 c -2 Î Nc c=1 ω 0 c f c + Î2 Nc c=1 ω 0 c + Nc c=1 ω 1 c g 2 c -2 Î Nc c=1 ω 1 c g c + Î2 Nc c=1 ω 1 c = Ω c ω 0 c f 2 c + ω 1 c g 2 c -2 Î c ω 0 c f c + ω 1 c g c + Î2 c ω 0 c + ω 1 c .
Making use of the expression of Î (6.11) a first time:

F = Ω c ω 0 c f 2 c + ω 1 c g 2 c -Î2 c ω 0 c + ω 1 c ,
and a second time after factoring by 1/ k ω 0 k + ω 1 k in the integral:

F = Ω 1 k ω 0 k + ω 1 k i c ω 0 i + ω 1 i ω 0 c f 2 c + ω 1 c g 2 c - i c ω 0 c f c + ω 1 c g c ω 0 i f i + ω 1 i g i = i c Ω ω 0 i ω 0 c f 2 c -f i f c + ω 0 i ω 1 c g 2 c + ω 0 c ω 1 i f 2 c + ω 1 i ω 1 c g 2 c -g i g c -2ω 0 i ω 1 c f i g c k ω 0 k + ω 1 k . Again as, for instance, i c ω 0 c ω 1 i f 2 c = i c ω 0 i ω 1 c f 2 i : F = i c Ω 1 k ω 0 k + ω 1 k ω 0 i ω 1 c (g c -f i ) 2 + 1 2 ω 0 i ω 0 c (f c -f i ) 2 + 1 2 ω 1 i ω 1 c (g c -g i ) 2 = c Ω ω 0 c ω 1 c k ω 0 k + ω 1 k (g c -f c ) 2 + c i =c Ω ω 1 c ω 0 i k ω 0 k + ω 1 k (g c -f i ) 2 + c i<c Ω ω 0 i ω 0 c k ω 0 k + ω 1 k (f c -f i ) 2 + c i<c Ω ω 1 i ω 1 c k ω 0 k + ω 1 k (g c -g i ) 2 .
(6.12)

Considering only the terms such that i = c in the previous expression of F (first term) allows to retrieve a functional similar to the one used in the usual frameworks for the displacement measurement step (see Equation (6.4) and Equation (6.5)).

Let us stress again that relying on a substitute image Î in a displacement measurement perspective shows the benefit to have a much richer functional than the usual ones. In Equation (6.12), there are indeed terms proportional to (f c -f i ) 2 and (g c -g i ) 2 which are similar to a shape measurement (see Section 6.1.2). The stereo correspondence is thus preserved. There are also terms proportional to (g c -f i ) 2 , i = c (spatio-temporal cross-correlations) which have no counterparts in the usual frameworks.

Data assimilation displacement measurements

In this subsection, we investigate the possibility to minimise the functional F defined in Equation (6.6) with respect to every argument (i.e. D t 0≤t≤Nt-1 , p c 1≤c≤Nc , Î) and show to which extent this functional is suitable for performing data assimilation by, once again, establishing the link with usual frameworks. Here, by data assimilation, we mean benefiting from all available observations to evaluate quantities of interest (that is, for instance in DIC, displacements, shape, camera parameters, substitute image). A key element in such an approach is the level of confidence associated to observations that we have already discussed.

In what follows, to reduce the amount of notation, we will simply write I t c instead of

I t c • P c • φ D t . With such notations, the expression of Î is simply (least-squares solution of F minimisation): Î = Nt-1 t=0 Nc c=1 ω t c I t c Nt-1 s=0 Nc i=1 ω s i .
Note that this expression for Î, stemming from the minimisation of F , is very close to the heuristic approach used in [Berny et al. 2018a]. In the context of heat haze effects [START_REF] Jones | [END_REF], relying on a substitute image based on all available pictures is essential as the confidence associated to the reference picture is low. Then we can develop F from Equation (6.6) (similar treatment as in Section 6.1.3):

F = Nt-1 t=0 Nc c=1 Ω ω t c I t c 2 -2I t c Î + Î2 = Ω t c ω t c I t c 2 -Î2 t c ω t c = Ω 1 r j ω r j t c s i ω s i ω t c I t c 2 - t c s i ω s i ω t c I t c I s i = 1 2 t c s i Ω ω s i ω t c r j ω r j I t c -I s i 2 .
Finally, F can be split in different parts (s = t and s = t):

F = t     s<t   c Ω ω s c ω t c r j ω r j I t c -I s c 2
Similar to a displacement measurement functional (6.13) This final expression for F , obtained when minimising with respect to every single argument of the functional, clearly establishes the link with usual frameworks. We can see that it includes usual shape measurements (see Equation (6.1)) at all times, together with terms similar to displacement measurements (see Equation (6.4) or Equation (6.5)) for all pairs of times, as well as spatio-temporal cross-correlations (comparing I t c to I s i , with cameras i = c, (spatial) and times s = t (temporal)). Once again, it is much richer than the usual functionals.

+ c i =c Ω ω s i ω t c r j ω r j I t c -I s i 2 Spatio-temporal cross-correlations   + c i<c Ω ω t i ω t c r j ω r j I t c -I t i 2 Similar to a shape measurement functional     .

Discussions

In order to establish the links above between F (6.6) and usual frameworks, we had to adopt the same experimental setups. That is, at all times, the number of cameras N c is the same and the cameras are assumed to remain in a fixed position all along the experiment. Let us stress that it does not have to be the case, and that the formulation proposed (6.6) is easily extended to arbitrary number of pictures at each time (N t c ), with moving cameras (p t c ). This would allow to consider experimental setups with cameras supported by robotic arms or even drones, for instance. For these reasons and others that we wish to illustrate in what follows, the functional proposed (6.6) opens up new perspectives in terms of experimental setups. It offers much more flexibility to the experimenter, while providing a much greater robustness, as it increases the amount of data for each problem (camera calibration, shape, displacement).

First, as already evoked, when considering reference state images f c and deformed ones g c that see totally disjoint regions of the ROI, the first term in Equation (6.10), equivalent to usual SDIC frameworks, becomes zero, as the product ω 0 c ω 1 c equals zero. This kind of situation totally incapacitates all DIC software (including SDIC and 2D-DIC). This may arise in the case of large rotations as described in Chapter 7. Yet, the Functional F (6.6) allows to naturally address this issue, thanks to the cross-correlation terms.

Also, when considering large strains, relying on all available pictures with a weight ω t c depending on the displacement field D t , would be particularly helpful to perform a finer sampling of the substitute image [Goldlücke et al. 2014, Tsiminaki et al. 2019]. For large positive strains, det ∇φ D t = det I + ∇D t > 1. This assigns a greater level of confidence to the image I t c , which is consistent with the better sampling achieved by the pixels in I t c of the ROI. In other words, in the case of large (positive) strains, it is unfortunate, in the current frameworks, to identify the substitute image in the reference state images only, as the information in deformed state ones is much more reliable.

On top of that, note that the weighting scheme together with the construction of the Functional (6.6) naturally provide a way to merge results from different times and different viewpoints in order to perform multiscale substitute image identification and, most importantly, multiscale displacement measurements. That is, cameras with different resolutions imaging the ROI [Passieux et al. 2015a]. Currently, the dialogue between measurements performed at two different resolutions is still an open problem.

Then, regarding camera calibration, some research works identify projection parameters on the sole basis of reference state images [Dufour et al. 2015b]. This camera calibration process, while convenient from an experimenter perspective, has the major drawback not to calibrate the whole volume spanned by the object which can result in a stereo correspondence loss. Identifying camera parameters based on the minimisation of Equation (6.6) would allow to calibrate the whole volume spanned by the ROI, precisely because the minimisation would be performed on all positions occupied by the object. Also, this would allow to avoid calibrating the stereo rig at different times (based on targets), as done to prevent temporal drift during long experiments. Regarding this matter, the last terms in Equation (6.13), similar to shape measurement functionals, turns out to be useful.

Finally, this formulation is particularly suitable for spatio-temporal regularisation and one could imagine making use of it to perform SDIC measurements during tests with a single moving camera, or a rotating object (in a tomograph for instance) in front of the fixed camera, relying on similar techniques as in [Jailin et al. 2017].

Conclusion

There are two main results associated to the developments presented herein. First, we established a link between functionals on which global SDIC usually relies and a functional based on the sum of errors between a substitute image and observations from all cameras, at all times: Functional F (6.6). We showed that the latter is much richer than the formers in a displacement measurement context. Based on the consideration of large displacements, camera calibration and stereo correspondence issues, we illustrated that all the terms usually discarded in SDIC frameworks are actually extremely useful. For this reason, Functional F (6.6) appears to (a) be well-suited to perform data assimilation in SDIC, as expected from the construction of this functional based on all available data, (b) stand for an interesting perspective in the formulation of the Stereo Digital Image Correlation problem, as it can be seen as a dense counterpart of bundle adjustment methods [START_REF] Delaunoy | Photometric bundle adjustment for dense multi-view 3d modeling[END_REF], Goldlücke et al. 2014, Sutton et al. 2009].

Sensitivity analysis

We wish to also discuss the use of a model. In Chapter 5, we based our approach on the sole basis of the understanding of the underlying physical principles and involved phenomena, which is an extremely important point to us. But there is another interesting outcome that we can now explain based on our previous developments. To do so, we will rely on the same assumptions and method as [Roux & Hild 2006, Appendix A]. Each image I t c is polluted by a random white noise γ t c , of zero mean, and variance σ 2 (same variance for all cameras at all times but the generalisation to different variances is straightforward and requires to normalise the functional F accordingly). In order to rely on a QGN scheme, the surface gradient [Goldlücke et al. 2014] of the radiance model ∇ Ω L e is required. In usual frameworks,

∇P c ∇I 0 c • P c is indeed substituted to ∇P c ∇I c • P c • φ D based
on the grey level conservation equation (3.35). In the PhDIC one, we may substitute (at least in the case of a good initialisation)

α t c ∇ Ω L e to ∇P c ∇I c • P c • φ D .
The QGN scheme writes: (6.14) where dD collects all degrees of freedom (including with respect to time) and

H PhDIC dD = b PhDIC ,
           H PhDIC = Nc c=1 Nt-1 t=0 Ω ω t c N T α t c ∇ Ω L e α t c ∇ Ω L e T N b PhDIC = Nc c=1 Nt-1 t=0 Ω ω t c N T α t c ∇ Ω L e L e -I t c • P t c • φ D t . (6.15)
∇ Ω L e is considered to be unaffected by noise. Note that, as L e is a model, this assumption is not as strong as the counterpart in [Roux & Hild 2006], where the same assumption is applied to the image gradient ∇I 0 c . In this case, b PhDIC is changed by a quantity:

δb PhDIC = - Nc c=1 Nt-1 t=0 Ω ω t c N T α t c ∇ Ω L e γ t c • P t c • φ D t δR t cPhDIC . (6.16)
In the following, we assume that ω t c does not depend on noise. This may be true in some cases (e.g. 2D-DIC and det ∇D t 1). This results in an expectancy δD = 0 (where • denotes the expectancy, in the sense of the mean over different noises). The measurement covariance matrix is given by δD δD T = (H PhDIC ) -1 δb PhDIC δb T PhDIC (H PhDIC ) -1 . We can then compute δb δb T in our case:

δb PhDIC δb T PhDIC = c t i s Ω Ω ω t c N T α t c ∇ Ω L e δR t cPhDIC (X) × ω s i N T α s i ∇ Ω L e δR s i PhDIC T (X ) dX dX = c t i s Ω ω t c N T α t c ∇ Ω L e (X) × Ω ω s i N T α s i ∇ Ω L e T (X ) δR t c (X)δR s i (X ) dX A t,s c,i (X)
dX.

(6.17)

And δR t cPhDIC (X)δR s i PhDIC (X ) = σ 2 δ ci δ st δ P s i • φ D s (X ) -P t c • φ D t (X) . Thus A t,s c,i (X) = σ 2 δ ci δ st Ω ω t c N T α t c ∇ Ω L e T (X ) δ P t c • φ D t (X ) -P t c • φ D t (X) dX = σ 2 δ ci δ st S t c N T α t c ∇ Ω L e T • φ -1 D t • β t c (x )δ x -P t c • φ D t (X) dx = σ 2 δ ci δ st N T α t c ∇ Ω L e T • φ -1 D t • β t c • P t c • φ D t (X).
Using this expression for A t,s c,i in Equation (6.17) yields:

δb PhDIC δb T PhDIC = σ 2 c t Ω ω t c N T α t c ∇ Ω L e α t c ∇ Ω L e T N •φ -1 D t •β t c •P t c •φ D t . (6.18)
As ω t c includes V t c • φ D t , we have β t c • P t c = I in the previous integral and it may be simplified as:

δb PhDIC δb T PhDIC = σ 2 c t Ω ω t c N T α t c ∇ Ω L e α t c ∇ Ω L e T N . = σ 2 H PhDIC . (6.19)
Finally the displacement covariance matrix is given by:

δD δD T = σ 2 H -1 PhDIC . (6.20)
When considering only one time step (N t = 1), and thus only one displacement measurement for instance, it is worth pointing out that relying on a model exhibits a strong advantage. Besides the possibility not to differentiate pictures, the measurement uncertainty is divided by a factor √ 2 when compared to usual frameworks [Roux et al. 2009, Eq. (12)], or close to such a factor [START_REF] Blaysat | [END_REF].

Let us stress however that noise in the images is not the only source of measurement error. Besides the already discussed issues associated with light, pattern induced bias (or fattening effect) and interpolation schemes also stand for sources of error [Sur et al. 2021].

Estimation theory

In this work, the functionals used to estimate parameters in the different inverse problems encountered (e.g. camera calibration, displacement measurements) have been presented only as a least-squares minimisation of the residual. Also, regularisation terms have been introduced as a way to increase the knowledge on the solution when the ratio available data over unknowns became too low. However, both may be justified from a mathematical viewpoint, provided that some assumptions are met. In the same way as [Szeliski 2010, Appendix B], we consider a general set of measurement equations (e.g. position of target points in pictures, grey level values). Assuming N m different measurements have been performed, the mth measurement vector is modelled:

∀m ∈ [[1, N m ]], y m = f (x, p m ) + γ m , (6.21)
where y m is the measurement vector corrupted by the noise γ m . f is the measurement model associated with a set of parameters x that does not depend on the measurement number m while p m does. Note that, by definition, f (x, p m ) is uncontaminated (or unbiased). Typically in Chapter 2, y m collects the positions of target points detected in picture m xmj j (see Equation (2.28)), the camera model P stands for the model f , the vector x collects the target point positions X, Y and Z together with the camera intrinsic parameters (p int and d), and p m is the extrinsic parameter vector p ext .

The main idea is then to maximise the conditional probability or likelihood L = p x | y (with y = (y T m ) T m ) of x as values for the sought parameters given the observations y.

Remark We implicitly assumed p m was known in the previous paragraph. If it is not the case (like for the photogrammetric calibration) replacing x by x, p (with p = (p T m ) T m ) in what follows allows to circumvent the issue.

In the following, we assume that the noise vector γ m is a zero-mean multi-dimensional normal (Gaussian) random variable with covariance matrix Σ m :

γ m ∼ N 0, Σ m . (6.22)
We also assume that noise vectors are independent. For a more general framework, the interested reader is referred to [Benning & Burger 2018, 3 Variational Modelling].

At this point, the problem is that Equations (6.21) allow, given x and the noise model, to access y. In other words, we can easily model the conditional probability p y | x . Hence we can make use of the Bayes' rule:

p x | y = p y | x p(x) p(y) . (6.23)
Based on aforementioned assumptions on noise vectors (zero-mean normal random variables and independence), it is indeed possible to compute p y | x [Szeliski 2010]: (6.24) where the probability p(γ m ) is simply:

p y | x = m p y m | x = m p y m | f (x, p m ) = m p(γ m ),
p(γ m ) = 1 (2π) N d,m /2 det Σ m 1/2 exp - 1 2 (y m -f (x, p m )) T Σ m -1 (y m -f (x, p m )) , (6.25)
and N d,m is the length of the measurement vector y m .

Maximum likelihood

Without any prior knowledge on x, a possibility is to maximise p y | x directly as we may write:

p y | x = p(x, y) = p x | y = L. (6.26)
Hence we look for:

x = arg max x L = arg max x Nm i=1 1 (2π) N d,m /2 det Σ m 1/2 exp - 1 2 (y m -f (x, p m )) T Σ m -1 (y m -f (x, p m )) = arg max x Nm i=1 1 (2π) N d,m /2 det Σ m 1/2 exp - 1 2 y m -f (x, p m ) 2 Σ -1 m , (6.27) where y m -f (x, p m ) Σ -1 m = (y m -f (x, p m )) T Σ m -1 (y m -f (x, p m ))
is the Mahalanobis distance [Szeliski 2010]. But it is also possible to minimise the negative log-likelihood:

x = arg min x (-log L) = arg min x E = arg min x 1 2 Nm i=1 y m -f (x, p m ) 2 Σ -1 m + 1 2 Nm i=1 log (2π) N d,m det Σ m Constant k = arg min x 1 2 Nm i=1 y m -f (x, p m ) 2 Σ -1 m , ( 6.28)
with E = -log p y | x the cost or energy [Szeliski 2010].

Note that the matrix Σ m -1 plays the role of a weight for each measurement error residual [Szeliski 2010]. For this reason, it is often called the information matrix as it allows to define a level of confidence to each measurement.

Camera calibration When using cameras with the same noise levels and the same algorithms to detect target points in pictures, each covariance matrix Σ m may be assumed equal to σ 2 d 1. This justifies the framework adopted in Chapter 2, where a least-squares minimisation was used.

PhDIC Here comes our final argument to first define residuals in the image planes. In Chapter 5, we justified it indeed first heuristically (pixels are the information quantum) and then mathematically, a posteriori, thanks to the consistent weighting term that appeared when integrating by substitution. Based on the image irradiance equation Equation (5.1), y m corresponds to the set of pixels corresponding to the ROI in a picture I m , the model f is α m L e •β m , x encompasses the BRDF parametrisation and lighting model, and p m includes the camera thoughtput α m together with the camera model parameters (see Chapter 2). Hence, the functional to minimise based on the negative log-likelihood is:

F = Nt-1 t=0 Nc c=1 x∈S t c I t c (x) -α t c L e • β t c (x) 2 (Σ t c ) -1 . (6.29)
Assuming a spatially uncorrelated noise, that is Σ t c = σ t c 2 1, yields:

F = Nt-1 t=0 Nc c=1 x∈S t c I t c (x) -α t c L e • β t c (x) 2 (σ t c ) 2 . (6.30)
If the same cameras are used (which is usually the case) and noise levels do not depend on time (which may correspond to camera temperature equilibrium throughout the experiment) we may assume ∀(c, t), σ t c = σ and thus we may simply use a least-squares minimisation:

F = Nt-1 t=0 Nc c=1 x∈S t c I t c (x) -α t c L e • β t c (x) 2 (6.31)
Eventually, a continuous formulation may be adopted [START_REF] Delaunoy | Photometric bundle adjustment for dense multi-view 3d modeling[END_REF],Goldlücke et al. 2014], as only such a framework allows to rigorously account for surface foreshortening for instance (see Chapter 5):

F = Nt-1 t=0 Nc c=1 S t c I t c (x) -α t c L e • β t c (x) 2 dx. (6.32)
This justifies a priori the approach developed in Chapter 3 as we obtain the exact same expression as Equation (5.12).

Prior models

In some cases, a strong prior knowledge on the solution is available and may be formalised as a probability p(x) in Equation (6.23). We may then write the negative logarithm:

-log L = -log p x | y = -log p y | x -log p(x) + log p(y)
Constant , (6.33) and x is sought as:

x = arg min x E -log p(x). (6.34)
This is called Maximum A Posteriori (MAP) and we may identify the same kind of structure as Tikhonov regularisation [Benning & Burger 2018] (see Section 4.1.2).

Discussions

The same kind of functionals as Equation (6.28) is sometimes called optimal in the usual DIC frameworks (see for instance [Hild & Roux 2020]). However, it is worth remembering that optimality should be understood in the sense of previous developments (maximum likelihood or MAP) which required numerous assumptions. All of them may not be fulfilled. Typically, noise in pictures may not exactly follow a Gaussian distribution and/or be spatially uncorrelated. One of the strongest assumptions is the possibility for the model to exactly explain observations, as the only source of error is assumed to be the Gaussian noise. Practically, it is rarely the case as, for instance, there may be outliers such as specular reflections. But most importantly, in usual DIC frameworks, no unbiased model was available until now, as measurements were performed based on image registrations.

Conclusion

This chapter allowed us to justify the use of the PhDIC framework. First we showed the improvements it could yield by establishing the link with usual frameworks. The PhDIC framework is indeed a general framework as it encompasses usual ones and offers the perspective to extend the scope of existing schemes. Then, we showed by a sensitivity analysis, with weaker assumptions, that the measurement uncertainty is expected to be divided by a factor √ 2 compared to usual methods. Eventually, the PhDIC approach was justified a priori, and we showed that it was optimal, provided that some assumptions are met, thanks to arguments stemming from estimation theory.

After these justifications from a theoretical viewpoint, we propose to apply this framework on a practical test case in the following chapter.

In this practical application of the PhDIC framework, the object is first 'scanned' thanks to numerous multiview pictures (I i ) i allowing to elaborate the test Digital Twin. In the Lambertian, distant point light source assumptions, it includes a shape correction field S, an albedo ρ (or intrinsic texture) together with light intensity and direction. Then, in the displacement measurement step t 1 , the deformed state images (J j ) j are compared to the model and the sought displacement field is identified. Note that both the number of deformed state pictures and the associated viewpoints may be different from the reference ones. These considerations are summed up in an illustration in Figure 7.1. We denoted the time associated to the displacement measurement step t 1 instead of t 0 + ∆t as in Figure 3.10. There are two reasons for that. First, time t 0 may correspond to a scan of the structure without mechanical loading. Second, the displacement measurement step t 1 can also stand for the positioning of the specimen in the test setup (which may apply loads that are usually neglected) and camera removal.

Intrinsics, extrinsics, shape and albedo measurements

In the present section, we showcase the calibration procedure prior to displacement measurement, the latter being, in turn, described in Section 7.2.

Before being able to perform a displacement measurement thanks to SDIC, several prerequisites must first be fulfilled (see Figure 3.8). The cameras should be calibrated (intrinsics and positions relative to one another), and the extrinsics and shape should be measured (position of the model with respect to the camera rig and corrections between nominal and actual shape). The difficulty in this prior phase concerns the shape measurement problem which is extremely ill-posed (see Figure 3.9). Regularisation strategies must therefore be adopted to circumvent this issue. They usually consist in restricting the subspace in which the shape is sought (whether it be in a strong or a weak sense) (see Section 4.1.2).

Another path that could be followed is an increase of the amount of available data [Passieux et al. 2015a], but as detailed by [Goldlücke et al. 2014], obtaining numerical schemes which scale favourably with the number of cameras is not straightforward (in the case of FE-SDIC, it scales as N 2 c as each picture has to be compared with every other one, see Section 6.1.1). This may explain why this possibility has not been fully investigated in SDIC yet.

Setup

In the present work, a single camera (N c = 1) was used to take multiview pictures of a rectangular plate with a circular hole (see Figure 7.2a) in the reference state (t 0 ). For this reason, we use the letter i to index the pictures (I i ) i taken by a same camera in the reference state. N i denotes the number of multiview pictures. The specimen was 20 cm long, 2.35 cm wide and 6 mm thick, while the diameter of the bore was 7 mm. A classic black and white pattern was created by spraying paint on the surface of the sample. A Jai GO-5000C-USB 5 Megapixel camera and a 25 mm macro lens were used. The distance between the sample and the camera was about 1m. The spatial sampling provided by the pictures was about 8 pixels per mm. Also, a single halogen light was placed right behind the camera so that a point visible by the camera was lit as well. The beam was attached on a custom calibration target composed of 8 points printed on an A4 sheet. The target was then fixed on a turntable allowing to take 360-degree pictures of the coupon, as indicated in Figure 7.2b. Let us stress that, with such a setup, the direction of the light with respect to the beam changes for each picture, while remaining the same in the camera reference frames. Also, the turntable was only a convenient way to take multiview pictures of the specimen. It served no metrological purpose. As described below, camera poses were rather identified thanks to the target.

A classic photogrammetric calibration [Garcia 2001] was performed on the pictures containing both the coupon and the target thanks to the in-house calibration software presented in Chapter 2. Usually, at the end of this step, the camera intrinsic parameters are identified and saved but the relative position of the target with respect to the camera, which is also one of the identified quantities, is discarded. Here, we use this knowledge to initialise the extrinsics and shape measurement procedure since it allows to directly estimate the pose of the images with respect to the coupon which is assumed to be fixed in the target reference frame.

We now place ourselves in a FE-SDIC framework [Passieux 2018]. A perfect CAD model of the specimen is first meshed using T3 elements (see Figure 7.2b again). The typical size for elements was 5mm but 2mm elements were used for mesh refinement care around the hole. The hole inner surface was not meshed.

Assumptions

We assume a Lambertian reflectance model for the beam surface. On top of that, we will assume that the only light contribution comes from the light source mentioned in Section 7.1.1 (thus, N r = 1 and we omit the index r in the following). We further consider that for each picture i ∈ [[1, N i ]], the light can be modelled thanks to an infinitely distant point light source with a vector n f (θ i , ϕ i ) given by the z-vector of the camera reference frame associated to picture i: Z c,i (in our convention Z c,i is the unit vector collinear with the optical axis and pointing from the scene towards the camera, see Figure 5.3). Hence, for image i, Equation (5.9) becomes:

I i (x) = ρ(X)l i Z c,i , n(X) . (7.1)
This infinitely distant point light source assumption is valid if the size of the coupon is negligible with respect to the distance between the coupon and the light, which was the case here (ratio of approximately one order of magnitude).

Since here we only focus on shape measurement, we assume that the available FE model allows to consider only slight corrections D 0 (see Equation (5.17)); in other words, the true shape is expected to be close to the nominal shape of the specimen. In this case, quantities in Equation (5.21) can be computed on the reference state geometry:

n(X + D 0 (X)) ∼ n(X), det ∇φ D 0 ∼ 1 and ((J i • P i ) V i ) • φ D 0 ∼ (J i • P i ) V i .
This simplifies the formulation and allows to compute once and for all the normal field and the weighting term on the reference geometry.

We further assume that the camera can be well described by a pinhole camera model, without distortions. This assumption, which is practically true with the optical system used in this work, comes with the benefit to have an analytical expression for J i (see Equation (5.14)) which can be computed exactly (see Appendix B) and

J i = f x f y Z 2 c,i O c,i M 2 |Z c,i | | n, r i |.
Remark Another benefit associated with the pinhole camera model is the smaller number of parameters and thus the need for a smaller number of calibration pictures.

As the light and camera are close to each other, the size of the coupon is also negligible with respect to the distance between the coupon and the camera. We can thus further simplify the

weight J i . Indeed ∀(X, X 0 ) ∈ Ωi vis 2 , O c,i M (X) 2 ∼ |Z c,i (X)| ∼ |Z c,i (X 0 )|
and since there is only one camera we do not need to consider the factor

f x f y Z 2 c,i
and may take J i = | n, r i |, as in [Birkbeck et al. 2006].

Probably the most challenging concept in this framework is the visibility function since its value at a point X depends on the camera position and orientation but also on the model geometry. To deal with this issue, we take advantage of the fact that the coupon, notwithstanding the hole, is a convex shape. Hence we can assess the value of V i (X) based on the sign of n(X), r i (X) whose absolute value is equal to J i • P i (X). Thus we use in the following [(J i •P i )V i ](X) = ( n(X), r i (X) ) + , where (•) + denotes the positive part function [Horn 1986, 16 Extended Gaussian Images].

Since the camera extrinsics (with respect to the coupon) have already been calibrated, we know the positions of the images relative to one another. Thus we can consider that we do not need to estimate each image pose with respect to the coupon but rather the position of the coupon (only 6 parameters) with respect to the virtual camera rig formed by the pictures (see Figure 7.2b). Ultimately the functional to be minimised writes:

F p ext 0 , S, ρ, (l i ) i = N i i=1 Ω ( n, r i ) + (X) I i •P i (X+D 0 (X))-ρ(X)l i Z c,i , n(X) 2 dX, (7.2)
which is very close to the functional used by [Birkbeck et al. 2006]. p ext 0 denotes here the extrinsic parameters of image 0 with respect to the coupon (all other images are then positioned thanks to the photogrammetric calibration). Here D 0 is the sum of the shape correction field S and the rigid-body displacement U ext associated to p ext 0 (D 0 = S + U ext ). Now that all required assumptions have been clearly stated and expressed in mathematical terms, it is possible to define the process used to minimise this functional.

Discretisation and interpolation

In order to compute the integrals in Equation (7.2) for instance, the specimen surface is discretised thanks to integration points and the integral over the physical domain is rewritten as a sum over all these integration points. The idea is to define integration points once and for all in the physical domain following the same strategy as in Figure 3.4, where as many integration points as the number of pixels in a finite-element are used. There is however a slight difference because of our multiview setup. Depending on the image used to observe it, the number of pixels in an element can be very different. For this reason, we decided to define an integration point density d in points/mm. The number of points along each direction of our T3 elements was then defined as the product of the density with the associated edge length. To choose the value for the density, we used the pictures where the greatest number of pixels was reached for a given physical area. That is, we made the opposite choice to [Dufour et al. 2015b] where the coarsest mapping is used. In our case, this led us to choose d = 8 points/mm.

To avoid undesirable oscillations at the free edges of the specimen (edges that belong to only one element), we removed integration points closer than ε to these edges [Baconnais et al. 2020]. In this work the chosen distance was ε = 0.5mm.

Regarding both picture subpixel interpolation and gradient computation, a regular bi-cubic spline interpolation was used.

Remark An important speed-up was obtained by computing the integrals only over Ω ivis (making use of the visibility function to exclude the integration points with zero weight).

Minimisation strategy

The process that we used to minimise (7.2) was a fixed-point algorithm consisting in an alternating optimisation algorithm. The reason for that is the problem ill-posedness. In addition to the usual sliding modes [Pierré et al. 2017], one should also be aware of the bas-relief ambiguities. We refer the interested reader to [Belhumeur et al. 1999] explaining (for orthographic projection models though) that a surface object is indistinguishable from a generalised bas-relief transformation of the geometry and an appropriate scaling of the albedo.

Before describing in detail the way it was implemented, we describe the search directions we used and how we managed to minimise with respect to each of these directions.

Extrinsics and shape

To minimise F with respect to extrinsics and shape, the extrinsic displacement field U ext and the shape correction field S were sought in subspaces of lower dimension than the linear span of the shape functions associated to the FE mesh. A Gauss-Newton iterative minimisation scheme together with a Ritz-Galerkin reduced order method, similar to those presented in Chapter 3, were used. Both S and U ext were written as a linear combination of elementary displacement fields (similar treatment as in [Colantonio et al. 2020, Equations (4) and(5)]).

• Extrinsics: One difficulty with this approach is that it is not straightforward to write the coupon rigid-body displacement as a linear combination, as rotations involve sine and cosine functions of rotation angles. To circumvent this issue, we used infinitesimal rotations around the centre of the coupon (the position is assumed to be well initialised).

The exact position of the hole was not known prior to experiment and is hard to measure experimentally, thus 2 additional degrees of freedom (translations) were added to be able to precisely measure its position. This step was needed to avoid nodes at the boundary of the hole to 'fall' into the void.

According to these considerations, we can denote U ext = N R ext p ext 0 where p ext 0 collects now both the extrinsics and the 2 additional degrees of freedom associated to the hole position.

• Shape: It is common to measure the shape correction along the normal at the nodes of the mesh [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF], Pierré et al. 2017, Chapelier et al. 2021]. Defining the normal at a node is not straightforward, and it is usually done by computing the mean over the normals of neighbouring elements. This definition is satisfactory for nodes located in the bulk of the surface, it is not when considering nodes located on an edge or a corner. In this work, we used a k-means clustering algorithm [START_REF] Jiawei | [END_REF]] to be able to detect nodes where 2 (edge) or even 3 (corner) different degrees of freedom were needed to consistently measure the shape. We defined a maximum value for the half-angle of the cone circumscribed to all the normals of a cluster. For each node, a k-means clustering algorithm was called with only one cluster over the set of normals of neighbouring elements. Then the number of clusters was increased until either there were 3 clusters or in each cluster the angle formed by all normals and the cluster centre was less than the maximum defined half angle. Finally, the node was affected the cluster centres as degrees of freedom. For instance, for a node located on an edge of our beam, this allows to measure a shape correction along both relevant directions, namely the normals of the planes intersecting in the edge. We can write S = N R shape s.

• Gauss-Newton algorithm: To minimise F with respect to extrinsics and shape, a Gauss-Newton algorithm was used [Pierré et al. 2017]. The Gauss-Newton update at each iteration results from the following linear system:

H D 0 dD 0 = b D 0 ; (7.3) with:              H D 0 = N i i=1 Ω ( n, r i ) + N T ∇P i ∇I i • P i ∇P i ∇I i • P i T • φ D 0 N b D 0 = N i i=1 Ω ( n, r i ) + N T ∇P i ∇I i • P i • φ D 0 ρ l i Z c,i , n -I i • P i • φ D 0 , ( 7
.4) where ∇I i denotes the image gradient. The Ritz-Galerkin method then writes: (7.5) in the case of extrinsics, R = R ext and dq = dp ext 0 , in the case of shape, R = R shape and dq = ds.

R T H D 0 R dq = R T b D 0 ,
Remark Compared to the usual framework [Pierré et al. 2017] (see Section 3.5.2), considering only p ext 0 as parameter and not every single p ext i allows to use the exact same algorithm as for the shape measurement to calibrate the extrinsics.

Remark Practically, we used a slightly different visibility function Vi , as in [Birkbeck et al. 2006]. We considered that a point was visible not only when n(X), r i (X) > 0 but when n(X), r i (X) > V isionT hre > 0. Because we found that results tend to be more accurate when increasing V isionT hre, the value V isionT hre = 0.4 was determined as the greatest possible value leading to a non-singular matrix H D (the top face of the beam was not 'seen' for higher values of V isionT hre). Thus, in the following, ( n(X), r i (X) ) + will rather stand for the product of (J i • P i ) and Vi . An interesting outlook here could be to compare the optimal value of V isionT hre (when there are enough observations to ensure that the matrix H D does not become singular) with the value corresponding to the optimal Stereo-angle found in [Balcaen et al. 2017].

Light intensity

We chose to calibrate the light once and for all. To do so, we arbitrarily set the value for the albedo of the white sheet standing for the target to 1. We then considered 4 points located on the sheet (see the blue dots in Figure 7.2a). These points are denoted by X light p p . An overdetermined system was solved in the least-squares sense for each picture i to retrieve l i thanks to Equation (7.1) by taking n(X light p ) = -Z w :

∀p ∈ [[1, 4]], I i • P i (X light p ) = -l i Z c,i , Z w . (7.6)
Remark We tried to use the same framework as [Birkbeck et al. 2006] where spheres are used to calibrate both the light intensity and light direction. This has the benefit to allow the identification of an ambient term, which is not the case here because the normals n(X light p ) are all the same (the matrix associated to the overdetermined system would not have full rank). The specular reflection can easily be detected and allows to obtain the direction of the source while the Lambertian part of the surface allows to get the other parameters. However, in our case this yielded poorer results than the method described above. We believe that this is because the infinitely distant point light source assumption was not completely valid. There was only one sphere and thus the light intensity information was only valid around the sphere. The 4 points used above allowed to obtain a less accurate but more general value for l i .

Albedo estimation This part is probably the easiest one since a closed-form solution for the albedo ρ minimising (7.2) can be derived (standard linear least-squares problem):

∀X ∈ Ω, ρ(X) = N i i=1 ( n, r i ) + l i Z c,i , n (X)I i • P i (X + D 0 (X)) N i j=1 n, r j + l 2 j Z c,j , n 2 (X) . (7.7)
Remark This definition for ρ(X) is in some sort a weighted average of all available observations of the physical point X. In this sense, it makes it similar to the definition of Î in [Dufour et al. 2015b] (actually denoted f in the cited work). However, considering both foreshortening and lighting effects shows the benefit to obtain a much sharper albedo, see Figure 7.3.

Remark

In order to obtain a speed-up in computation time, no interpolation scheme was used to evaluate the numerator in Equation (7.7), that is the nearest neighbour pixel was used to evaluate I i • P i (X + D 0 (X)). No significant changes in the identified shape nor albedo were observed regardless of the interpolation scheme.

(a) Texture identified at initialisation without accounting for lighting effects (similar to f in [Dufour et al. 2015b]). Results are presented in grey levels.

(b) Albedo (dimensionless) identified at initialisation accounting for lighting effects.

Figure 7.3: Texture and albedo identified at the initialisation step. Note that the units are different: we make a distinction between the texture (in grey levels) and the albedo (dimensionless). To avoid bias, the scale chosen is the amplitude of each data set. Accounting for lighting effects clearly results in a much sharper gradient for the albedo.

Alternating optimisations

The structure of the iterative algorithm used herein to minimise the functional is presented in Figure 7.4. On top of the procedure detailed above consisting in minimising with respect to different variables, we also made use of a multiscale (or coarse-graining) initialisation process (see Section 4.1.1). At the beginning, the discrepancy map D 0 was initialised to 0 and pictures were considered at a scale s = N s = 3. At the scale s, pixels in the initial pictures were aggregated by groups of 2 s × 2 s resulting in coarser images. The density of points d introduced in Section 7.1.3 was set accordingly (d/2 s ). A Tikhonov regularisation term [Pierré et al. 2017] was added to the functional (gradient L 2 -norm squared). This term had a decreasing amplitude with the scale, until no regularisation was used for s = 0.

Two different values were used as stopping criteria, namely stagnation with respect to the discrepancy map and with respect to the functional, defined respectively by Res D = 10 -5 and Res F = LoopRes F /4 s with LoopRes F = 10 -3 in Figure 7.4. Dividing LoopRes F by 4 s for each scale allows to demand a better precision at the fast-to-compute coarsest scales which are known to be less accurate.

Results

In Figure 7.5 are shown the initialisation and convergence states with the mesh superimposed on pictures. We can see that our method allows to recover the specimen shape even though the object size was overestimated. Note that, in this framework, there is neither a need for selecting points 'by hand' or automatically [Passieux et al. 2015a, Pierré et al. 2017] nor for fiducial marks in the pattern on the object. The regularisation of the extrinsics and shape measurement problem was rather obtained thanks to much more images than usual and the multiview setup associated to the 3D mesh which allow to measure the specimen edges.

The total number of available pictures of the specimen was 72 (approximately 5 • between each pose). To evaluate the methodology described herein, we decided to form two different disjoint sets of pictures. Each of them contained 36 pictures with approximately 10 • between each pose (see again Figure 7.2b). This allowed us to apply independently our method on these The mesh is superimposed on pictures.

two sets to compare the identified shapes and albedos. The results presented in Figure 7.5 were obtained thanks to one of these two sets.

Camera calibration In Figure 7.6 is shown the reprojection error standard deviation associated to the calibration step described in Section 7.1.1, for one of the 36-picture sets. These results are satisfactory since the total reprojection error standard deviation is equal to 0.18 pixel. Albedo In order to compare the retrieved albedos between the two different 36-picture image sets, we decided to compute the normalised albedo difference, defined as follows: 2(ρ 1 -ρ 2 )

ρ 1 +ρ 2
, where ρ 1 and ρ 2 stand for the albedo associated to each image set, and ρ 1 and ρ 2 denote their mean values. The distribution of this quantity is plotted in Figure 7.7. The difference mean value (0.004) is small compared to the difference standard deviation (2.7%). We chose to compare this last quantity to the normalised camera noise level, since we believe it is the relevant quantity to compare the normalised albedo to, in the same way as we will compare the shape measurement error to the calibration reprojection error (see Section 7.1.5). We estimated the camera noise thanks to 9 pictures taken for 6 different poses (54 pictures in total) and obtained a normalised camera noise mean and a normalised camera noise standard deviation respectively equal to 1.9 × 10 -16 and 1.3%. We can see that the standard deviations of the normalised camera noise and normalised albedo differences are of the same order of magnitude. However, we will see later on, that the camera noise is not the only error source that we identified. Shape In Figure 7.8 are shown the projections of the integration points coordinates difference along each direction of the world reference frame (see Figure 7.2a). For the X and Z direction, the mean value is small compared to the standard deviation associated to the projected coordinates difference. Distinguishing the 3 directions allow to see a particularity in the Y direction since .7: Distribution of the normalised albedo difference (defined at integration points) for the two different sets of pictures at convergence state. The mean value is small compared to the standard deviation. It should be noted that camera noise is of the same order of magnitude as the albedo difference standard deviation. in this case the mean is much greater than the standard deviation. We interpreted these results based on the solution non-uniqueness and an illustration is shown in Figure 7.9. Basically, the idea is that the framework described herein does not prevent from converging to different parameterisations of a same geometry. Thus, integration points can describe the same geometry and lie at different places on the surface of the object, even though they were defined at the exact same place on the initialisation mesh. This can also partly explain the slightly larger standard deviation associated to normalised albedo difference than the one associated to the normalised camera noise. Each integration points stands for an albedo at a slightly different place in each one of the considered image set. The reason we do not end up with odd results in Figure 7.7 is thanks to the pattern which smoothly varies in space.

Displacement measurement

In this section, we aim to show, through an easy-to-setup test case, the potential of the developed PhDIC methodology to measure displacement fields that would be extremely delicate, if not impossible, to capture with the usual SDIC framework. Before really entering into the details of the specific test case, let us present the general way to use this framework to perform arbitrary displacement measurements.

Ground-breaking methodology for displacement measurement

To underline the novelty of the approach, we will denote the deformed state images by (J j ) j . Note that since we are almost able to generate virtual images (with the digital twin) to compare our observations to, they do not need to be associated to so-called reference state images I i (see Figure 7.1 as a reminder). Among others, this implies that the number of deformed state images N j does not have to be equal to the number of reference state images N i and can be greater or Figure 7.8: Distribution of the coordinates difference (defined at the integration points) for the two different sets of pictures at convergence state. Each direction is plotted on a different graph and the reference frame chosen to define X, Y and Z is the beam (or world) reference frame. The unit chosen for measuring a distance is the pixel in order to compare it to the calibration reprojection error. One should keep in mind that it makes only little sense to measure distances in pixels, as in a general multiview setup it is not straightforward to convert a distance from mm to pixels. Here, because the camera stood always about the same distance from the beam, an 8 pixels/mm constant of proportionality was used to convert the measures. For X and Z directions, the mean value is smaller than the standard deviation. This last quantity is of the same order of magnitude as the calibration reprojection error (0.18 pixel). Regarding the Y direction, we can clearly see a bias (∼ -0.15 pixel), much larger than the standard deviation, corresponding to a shift along the Y direction. We suggest an interpretation of this bias in the measurements in Figure 7.9 based on the solution non-uniqueness. Figure 7.9: From a same mesh (black one), one can obtain two different parameterisations of the surface. We can interpret that by the coupling between shape and extrinsics measurement. The black mesh, corresponding to the initialisation step, is slightly larger than the real object. During the extrinsics calibration step, the mesh can slide indifferently in one direction or another, which is represented respectively by the yellow and purple mesh. Finally, the same shape is measured, but integration points lie at different positions. It is due to the problem ill-posedness (solution non-uniqueness). The crosses stand for the integration points. smaller, and also can correspond to different camera poses.

Remark This is what was suggested in Section 6.1.4: N 0 c = N i and N 1 c = N j .

For the displacement measurement, the associated part of Functional (5.18) (generalised with different number of cameras per time step N t c ) writes (Lambertian model and distant point light sources):

F = N j j=1 Ω det ∇φ D 1 J j • P j V j • φ D 1 (X) × J j • P j (X + D 1 (X)) -ρ(X) Nr r=1 l r n f (θ r , ϕ r ), n(X + D 1 (X)) 2 dX. (7.8)
Remark We chose to update all variables identified in the previous section. More precisely, in the equation above, Ω stands for an update of Ω by φ D 0 , and ρ together with n for the associated albedo and normal vector field. This allows to obtain the usual definition for the displacement D 1 (that is on the reference geometry, not on the nominal one).

Remark As we consider D 0 to be identified in the previous step (see Section 7.1), the above equation discards t = 0 compared to the whole functional Equation (5.18).

With the same assumptions as in Section 7.1 regarding the weights J i and the light model, we can simplify Functional (7.8). We also make explicit the variables with respect to which F will be minimised:

F D 1 , l j j = N j j=1 Ω det ∇φ D 1 n, r j + • φ D 1 J j • P j • φ D 1 -ρl j Z c,j , n • φ D 1 2 .
(7.9)

Remark In this proposed implementation, ρ is no longer an unknown since we consider it has been identified in previous Section 7.1.

At this point, we suggest two ways to proceed. Either the integral is rewritten on an initial Ω, which is very close to the actual deformed geometry, and the same assumptions as in Section 7.1

apply (n(X + D 1 (X)) ∼ n(X), det ∇φ D 1 ∼ 1 and n, r j + • φ D 1 ∼ n, r j +
) or we keep integrating on the reference geometry Ω as in Equation (5.21). Choosing the first approach may lead to two main benefits. First, the very same Gauss-Newton minimisation scheme as in Section 7.1.4 can be employed because a satisfactory initialisation does require to be built as it is a prerequisite for any gradient-based iterative minimisation scheme. Second, choosing the second approach is extremely costly computationally speaking, since, for instance, the normals and the visibility function should be constantly updated. However, this allows to obtain discrepancy maps D t t with consistent definitions as their definition domain remains the same (that is Ω). Also, assuming D 1 0 is a first good guess, we may approximate n(X + D 1 (X)) ∼ n(X + D 1 0 (X)), for instance, and thus compute the normal vector field once and for all, instead of constantly updating it.

Remark Note that in both cases, the functional is unchanged. What changes is the integration domain.

Remark

The two proposed resolution procedures can be viewed as the counterparts of the two main variants in computational solid mechanics to perform geometrically non-linear analysis, i.e. the updated Lagrangian (first method) and the total Lagrangian (second method) strategies [Bouclier et al. 2015, ten Thije et al. 2007[START_REF] Oliver | [END_REF].

However, for the sake of pedagogy, we chose to follow the second approach, that is integrating over the reference domain Ω, as it allows to showcase a slightly different minimisation algorithm.

Example of a large rigid-body displacement measurement

The simple but illustrative test case that we considered was a large rotation (90 • ). In a usual SDIC framework, this would yield to two whole faces from the reference image replaced by two others as shown in Figure 7.10a and Figure 7.10b, and thus would make the usual SDIC framework fail. We can further simplify Equation (7.9). As we want to measure a rigid-body rotation, we have det ∇φ D 1 = 1. Also, the light intensity l j is obtained in the same way as in Section 7.1.4. Finally, we used only one picture (N j = 1), since we theoretically have only 6 degrees of freedom, a single image should yield enough information.

Remark This is because we used only one picture for the deformed state that we performed the minimisation on parts of the whole Functional Equation (5.18). One more picture does not yield a lot more information than already available. Note that this process (Sections 7.1 and 7.2) could be used as an initialisation for the minimisation of the whole Functional (5.18).

Once again, a Gauss-Newton minimisation algorithm is used:

H D 1 dD 1 = b D 1 , (7.10) with:                    H D 1 = N j j=1 Ω N T n, r j + ∇P j ∇J j • P j ∇P j ∇J j • P j T • φ D 1 N b D 1 = N j j=1 Ω N T n, r j + ∇P j ∇J j • P j • φ D 1 × ρ l j Z c,j , n • φ D 1 -J j • P j • φ D 1 .
(7.11) Since an acceptable initialisation is needed for the Gauss-Newton minimisation algorithm to work, we kept on using the same linearisation for the extrinsics measurement (R ext ) as in Section 7.1.4. For initialisation purposes, we computed the position of the camera relatively to the target thanks to the same kind of algorithm as the calibration one. At each step (each time a new dD 1 was computed), both n and r j were updated. The algorithm successfully converged with a single picture as long as we initialised with a displacement field D 1 0 corresponding to a position of the mesh not further away than 4 • from the actual specimen (see Figure 7.10c).

Remark This rigid-body rotation measurement can also be seen as a repositioning of the camera with respect to the object. Thus, it opens the possibility of experimental setups where (some) cameras move around the object during tests. Figure 7.10: This Figure allows to explain the approach followed and how difficult it would be to measure this 90 • rotation in a usual SDIC framework. As the displacement measurement step in classic SDIC would be based on the comparison between reference state images I 0 i 7.10a and deformed state ones I i 7.10b that observe different faces of the specimen.

Conclusion

A downgraded version of the PhDIC framework (in the sense that a minimisation is performed for each time step, independently as in usual frameworks, and not for the whole functional), in the case of a Lambertian model and a single distant point light source, was described and applied to a parallelepiped beam with an open hole. A turn table was used to scan the object thanks to 360 • multiview pictures. This allowed to recover both the shape and albedo of the beam. For this step, the in-house calibration software that we implemented, and of which details are provided in Chapter 2, turned out to be particularly useful. The results between two completely different sets of input images were compared, demonstrating the accuracy and robustness of the approach.

Of course, building a textured digital twin comes at a price, but this price is related to the increase in the amount of available data. In return, associated data assimilation allows a strong regularisation of the shape measurement step. This step, which may be very delicate in usual Global SDIC frameworks, is critical for performing accurate displacement measurements.

Conclusion and Perspectives

In this last chapter, we recall the main points and results of the present manuscript as well as the perspectives that have been brought. The future works that will be engaged during the Photometric stereo Digital Image Correlation and full-field measurement Techniques (PreDICT) project are also introduced.

In a few words, I wrote the PreDICT project during the final year of this PhD thesis for several reasons. First, I realised the vast potential of the PhDIC framework, while at the same time being perfectly aware that this final year was far too short to engage in all the paths that I wanted to explore. Second, the price for all these possibilities is that the present research work is still a proof of concept, quite far from a demonstrator able to operate in an industrial context. That is why I took my chance and applied to the Price of the 'Fondation Jean-Jacques et Felicia Lopez-Loreta pour l'excellence académique', in the hope that the time invested could generate even more time in the future to work, with a project team, on these fascinating topics. I am particularly grateful to the 'Conseil de la Fondation', which decided to give this Price to, and thus funded, the PreDICT project.

Summary

In this manuscript, we presented DIC as a measurement method allowing to extract more data from large-scale tests than usual strain gauges for instance (see General introduction). After having introduced some prerequisites (Chapters 1 to 3), we have underlined that making DIC operate on such cases with complex geometries is a challenging task (see Chapter 3).

There are two main identified issues. First, complex geometries are associated with inhomogeneous strain states and thus may require multiscale measurements, that is cameras observing the ROI with different resolutions. This allows to obtain both local data in high-strain zones and global data regarding the boundary conditions. However, we showed that available patterns could compromise the use of a multiscale approach, given the fact that they are often adapted to a single scale (or, at best, two scales). Also, using multiple cameras and making images from different resolution levels contribute appropriately to measurements was still an open problem. Usually, measurements were performed on a scale by scale basis (that is constituting different functionals, one for each scale), and then, for all points, only the most resolved available scale was considered. In such a process, a lot of information is basically thrown away, which is unfortunate regarding the low number of tests performed at the largest scales.

Second, much greater displacements are expected at a given strain on large structures. It results in assumptions of the usual frameworks becoming harder and harder to meet as the scale increases. For instance, the assumption that a point may remain visible all the time by the same camera strongly restricts DIC operation domain. Another aspect of large displacements is associated with severe changes in the surface normal orientation as time goes on. Depending on the light distribution, this may result in grey level variations over time of a same point, but most importantly in changes in the surface sampling performed by cameras and thus of the level of confidence that can be assigned at the different time steps considered.

Working on DIC formulation in a multiscale context was thus a necessary long-term objective. But before following it, we first ensured that these possible theoretical developments were not pointless by defining multiscale patterns able to operate at a wide range of scales (see Chapter 4). For one of them, we showed that its intrinsic statistical scale invariance results in a measurement error independent of the considered scale as well, which was not the case for usual patterns.

Then, by a thorough physical modelling and explicit assumptions, a general DIC formulation has been proposed (see Chapter 5). It is based on CV and the definition of a Digital Twin of tests. This Digital Twin includes a refined geometric model of the structure together with its albedo (or a parametrisation of the BRDF, in a more general framework), and a light model. It shows the benefit to naturally tackle the issues of visibility, weighting of observations (scale, foreshortening, surface pixel sampling) and thus of data assimilation. Another benefit is the possibility to develop new formulations based on this work by questioning the soundness of some assumptions, that is to be able to model physical phenomena.

In Chapter 6, we justified the wording 'general' of the proposed formulation as we showed that it not only encompasses usual frameworks but also allows, in particular, to perform measurements when the reference and deformed state images are totally disjoint regions, thanks to spatio-temporal cross-correlation terms, and prevents temporal drift, or stereo correspondence loss. Also, as a dense counterpart of Bundle Adjustment, it would allow to reduce the uncertainty associated to the framework by considering one single problem instead of several nested problems (see Figure 3.8) where uncertainties accumulate. Thanks to a sensitivity analysis based only on displacement measurement, it was demonstrated that a reduction in the measurement uncertainty is expected when relying on a Digital Twin (or model). The possibility to differentiate this model instead of pictures was also discussed, as usual sensitivity analysis of QGN schemes based on the image gradients do not address this issue. Eventually, the choice for the functional definition to integrate first in the images was justified from the estimation theory viewpoint. It allowed to discuss the conditions of validity regarding the optimality of this definition (e.g. model that perfectly fits data up to noise level, Gaussian noise distribution).

In the last chapter (Chapter 7) a practical test case following the shape measurement and displacement measurement steps was presented. Results regarding shape, albedo and displacement measurement have been presented and discussed. This allowed to show practically that light and visibility issues are naturally addressed with the proposed formulation. Typically, it was demonstrated that this methodology allows to measure a (rigid-body) motion that totally incapacitates existing SDIC frameworks.

Perspectives Test cases

Regarding the choice of the large 90 • rotation as a displacement to measure: it was convenient in the context of the COVID-19 lockdown as it did not require any testing machine, and as I could set the experiment in my living room. However, a more industrial test case would allow to better showcase our work. First, one could think of tackling large deformations with large strains that incapacitate usual frameworks or require special treatments because reference state images do not sample the initial surface in the same way as the most deformed ones (e.g. tire sidewall correlation in the context of impact on panels, elastomer traction). Another possibility is to tackle large deformations in a test case where visibility becomes the main issue (for instance Carpentier Joint based hinges [START_REF] Pellerin | [END_REF], or buckling modes [START_REF] Jacquet | [END_REF]).

Clearly, PhDIC is based on radiometric principles like IRT. Coupling both measurement methods in one framework could be interesting [START_REF] Jones | [END_REF]. Usually, IRT assumes a test surface of uniform emissivity, however DIC requires a pattern to provide reliable measurements. Hence, depositing a pattern in the infrared spectrum and then characterising the emissivity at each surface measurement point in the same way that we measured the albedo in Chapter 7 may allow to both measure a displacement field and a temperature one, using only infrared cameras.

DIC, as an optical measurement method, is subject to physical phenomena which may lead to large measurement biases or even may prohibit its use in challenging experiments [START_REF] Jones | [END_REF]] (e.g. fluid-structure interaction applications because of refraction, or hightemperature mechanical tests because of heat haze effects [Berny et al. 2018a[START_REF] Jones | [END_REF]), taking into account these effects in the DIC framework by identifying them could allow to overcome the issues faced. For instance, following the steps of [START_REF] Cassidy | [END_REF] would enable the experimenter to identify the interface and follow an underwater test. Following a method similar to [START_REF][END_REF] would make it possible to identify refractive index gradients and be able to cope with heat waves in DIC.

Of course, as the PhDIC formulation is well-suited to address camera repositioning (see Chapter 7), having drones or robotic arms equipped with cameras [Ben Abdallah et al. 2019, Kalaitzakis et al. 2021] stands for the next step in DIC. This would allow, coupled with multiscale processes, to sweep the test surface, prior to tests, to scan it and measure its shape and reflectance properties. During tests, some cameras may carry on skimming the surface, looking for hot spots that could be focused on. Yet, this scan prior to tests might not be very practical for the experimenter and using different light positions while keeping cameras at the same position instead, as done in Shape-from-Shading [START_REF] Mélou | [END_REF], could be more relevant depending on the experimental setup.

In some cases, the assumption that the pattern is expected to exactly follow the deformation of the specimen and that its texture does not depend on the displacement or strain level (5.20) is not valid. In large strain contexts, the texture of a point may indeed depend on the strain level. For instance, measuring displacements on lumbar belts [START_REF] Molimard | [END_REF]] is complex because of such effects. Developing a law for ρt depending on the strain level is an interesting outlook that could help regularising the displacement measurement problem.

Digital Twin of tests

Concerning the Digital Twin, a lot of work remains to be done. As explained earlier, in the context of large-scale tests, issues associated with light may become important (see Figure 7.11) thus it could be interesting to integrate a BRDF parametrisation including not only a diffuse (or Lambertian) component but also a specular one as in [Birkbeck et al. 2006] for instance. Other light models may also be used. The ultimate goal being to find reflectance and light models suitable to a wide variety of large-scale tests. However, this may not be the only way to deal with light problems. By keeping a Lambertian model for instance, in the context of light issues, the assumptions of Section 6.3 fall apart, as the model does not perfectly explain the observations and it may be relevant to consider other residual norms (e.g. L 1 ) less sensitive to outliers. Another benefit associated with changing the residual norm used from L 2 to L 1 would be the better performance of a super-resolution framework [Mitzel et al. 2009, Figure 7]. Besides the fact that it is a beautiful application of inverse problems, it opens up perspectives to extract even more information from pictures [Goldlücke et al. 2014]. Also, it could be interesting to evaluate whether a super-resolved texture allows to obtain a super-resolved displacement measurement as well.

In the present work, no multiscale Digital Twin of tests has been defined, as pictures were all taken at the same resolution. It could thus be interesting to use this framework in a test involving a multiscale aspect. A lot of questions would require to be answered such as the efficient dialogue between data and the Digital Twin. We believe this would find an answer by being able not only to compare integration points grey level value with a prediction, but to generate virtual images (see [Balcaen et al. 2017, Lava et al. 2020, Rohe & Jones 2021] for the generation of virtual images in DIC) to directly compare, pixelwise, actual ones with. Carefully addressing the integration performed by a photosensor (and more generally by the optical system) to give a pixel is probably a key element regarding these questions. For instance, issues like unfocusing could be addressed and, again, more information could be extracted from pictures, as it has been shown that accounting for these effects allowed robust displacement measurements for low depth of field systems [Réthoré et al. 2014]. Typically, the transition in Section 6.3 between the discrete functional, optimal in the sense of the estimation theory, and the continuous one, used herein, neglected these aspects. The generation of virtual pictures of a scene would also allow to define a multiscale pattern in the physical space instead of in the image one, to be able to print it easily. In Chapter 4, we discarded indeed the process associated with image acquisition.

Algorithms

Obviously, the proposed test case is still too close to usual frameworks and does not take full advantage of the proposed formulation. A lot could still be done and the perspective to see Chapter 7 as an initialisation for the whole PhDIC functional minimisation is promising. As already discussed, it should be possible to calibrate cameras on the object (including intrinsics, and possibly distortions) [Dufour et al. 2016] as well as light model parameters and identify the shape and displacements by resorting on only one functional and thus one final minimisation. This would allow to reduce measurement uncertainties compared to those associated to nested optimisation problems. Investigating the possibility to derive the model instead of pictures in the minimisation schemes is also an interesting associated perspective to reduce uncertainties.

There is an important aspect that we did not discuss yet, the computational cost. As explained in the General introduction, we believe that DIC, and especially local approaches, focused on developing fast running algorithms. Here, we based the developed approach on physical considerations which may drastically increase computational costs (bear in mind that, ultimately, quantities stemming from Chapter 7 would only be an initialisation for the whole Functional (5.18) minimisation). Hence, we adopt here an orthogonal viewpoint where opposing local and global approaches makes little sense for they should rather be seen as complementary tools. For instance, local approaches speed and robustness [Le Besnerais et al. 2016] could be used to initialise the PhDIC framework on small patches of a large structure that could not be instrumented by standard DIC procedures (because of light and/or visibility issues). A parallel may be drawn here with [Yang & Bhattacharya 2019]. Another possibility is to couple the PhDIC framework with High Performance Computing (HPC) tools in order to obtain results faster [START_REF] Bouclier | [END_REF],Bouclier 2020,Merta et al. 2014,Passieux et al. 2015b[START_REF] Ronovskỳ | Elastic image registration based on domain decomposition with mesh adaptation[END_REF], Tournier et al. 2019, Wang & Kemao 2018].

PreDICT

Eventually, the PreDICT project is probably the most exciting perspective of this work, as it will strive for developing a large number of the aforementioned points. During this five-year project three different PhD thesis topics will be proposed. They break down as follows:

• visibility function, large-deformation/large-strain test cases and whole PhDIC functional minimisation,

• light model and BRDF parameterisations, IRT-DIC coupling,

• super-resolution and investigations regarding residual norms.

Also two research engineers will complement the project team all along the five years. The first one, will support the PhD students research work while at the same time develop algorithmic strategies to allow moving cameras during tests. The second one will provide support regarding numerical efficiency and software development, as a final demonstrator able to operate on largescale applications should be implemented by the end of the project. Personally, I will lead this project towards its final objective, that is the demonstrator, while obviously relying on the experience and expertise of those who helped me go further than I could imagine.

Appendices

Appendix A

Plane-stress state and elastic energy minimisation

In the plane-stress state assumption, the displacement u = (u, v) T is a function of only two spatial coordinates, say x and y. Also, using the Voigt notation:

σ =    σ x σ y √ 2σ xy    = C ε = C    ε x ε y √ 2ε xy    with C = E 1 -ν 2    1 ν 0 ν 1 0 0 0 1 -ν    .
(A.1)

We also assume no dynamic effect and the absence of volume force, thus:

∇.σ = 0. (A.2)
Using Equation (A.1) in previous equation yields: , we obtain:

           ∂ 2 u ∂x 2 + ν ∂ 2 v ∂x∂y + 1 -ν 2 ∂ 2 v ∂x∂y + ∂ 2 u ∂y 2 = 0 1 -ν 2 ∂ 2 v ∂x 2 + ∂ 2 u ∂x∂y + ∂ 2 v ∂y 2 + ν ∂ 2 u ∂x∂y = 0 . (A.4)
We can rewrite it so as to reveal the displacement Laplacian:

           ∆u + 1 + ν 2 ∂ 2 v ∂x∂y - ∂ 2 u ∂y 2 = 0 ∆v + 1 + ν 2 ∂ 2 u ∂x∂y - ∂ 2 v ∂x 2 = 0 (A.5)
In the case where ν = -1, we thus have ∆u = 0 over S which is the strong formulation associated with minimising 1 2 S ∇u 2 Following the exact same steps for the derivatives along the v-direction we obtain: At this point, we can compute: In each Sub-figure, the dependency of the physical length corresponding to a given pixel length in the image plane on one of the contributions is illustrated. For this, the orange dash-dotted line and the green dashed line have the exact same values for all parameters except those from the considered contribution, also both correspond to the same length in pixel in the image plane. The orange dash-dotted line stands for a unit contribution. Thus, the ratio of the orange dash-dotted line to the green dashed line lengths allows to retrieve the contribution.

F .        X c,u = - Z c + Z c,u (u -u 0 ) f x = - Z c f x + Z c,u X c Z c Y c,u = - Z c,u (v -v 0 ) f y = Z c,u Y c Z c (B.
                         X c,v = Z c f y n, e c
β ,u × β ,v = O w M ,u × O w M ,v = O c M ,u × O c M ,v = Z 2 c f x f y   n,
J = β ,u × β ,v -1 2 = f x f y Z 2 c O c M 2 |Z c | | n,
In Figure B.1a, we have similar triangles and thus the ratio of the orange dash-dotted line to the green dashed line lengths is f Zc . In Figure B.1b and Figure B.1c, one must recall that J is defined locally. Hence, both the orange dash-dotted line and the green dashed line should be considered as infinitesimal and the solid blue rays going from the optical centre to the edges of these lines can be considered parallel, although it is obviously not the case in the illustrations. For this reason, the ratio of the orange dash-dotted line to the green dashed line lengths is | n, r | in

Résumé long

Nous avons choisi d'inscrire ce travail dans le contexte aéronautique et celui de la pyramide des essais. Les cycles de conception aéronautiques sont en effet basés sur de nombreux essais aux plus basses échelles (par exemple, caractérisation d'éprouvettes élémentaires), et leur nombre diminue à mesure que l'échelle considérée augmente, jusqu'à n'avoir quasiment plus qu'un seul essai à l'échelle d'un avion complet (d'où la forme de pyramide). En raison de ce nombre d'essais qui diminue et d'une densité de mesure faible (du fait des coûts et du temps d'instrumentation), les modèles deviennent de moins en moins représentatifs des structures réelles au fur et à mesure de la montée dans les échelles. Des coefficients de sécurité sont donc régulièrement introduits, ce qui peut être vu à la fois comme un frein à l'innovation mais aussi comme un risque potentiel de sécurité. Nous présentons la Corrélation d'Images Numériques (CIN) comme une méthode de mesure permettant d'extraire davantage d'informations des essais aux plus grandes échelles par rapport à des mesures ponctuelles comme des jauges de déformations par exemple. Il s'agit d'une méthode de mesure optique, peu coûteuse. Elle repose sur la mise en correspondance d'images de la région d'intérêt dans un état de référence avec des images déformées dans le but d'extraire le champ de déplacement. Cette opération est facilitée par un mouchetis déposé sur la surface de mesure. Lorsque plusieurs caméras sont disponibles, des surfaces éventuellement non planes peuvent être considérées et un champ de déplacement tridimensionnel peut être mesuré. On parle alors de Stéréo Corrélation d'Images Numériques (SCIN). L'intérêt d'utiliser ces mesures surfaciques denses est de pouvoir mieux alimenter les modèles numériques et d'améliorer leur capacité à prédire le comportement des structures aux grandes échelles. Toutefois, un certain nombre de limitations de ces méthodes et notamment de la SCIN entravent l'extension de leur utilisation sur des essais aux géométries complexes.

Dans un premier temps, nous présentons dans ce manuscrit l'état de l'art dans les domaines de la CIN et de la SCIN, en particulier par le prisme des méthodes globales, et plus spécifiquement celui de la SCIN éléments finis. Cela nous permet de détailler les modèles de caméras utilisés, et qui sont présents dans bon nombre d'implémentations de ces méthodes car ils permettent de relier la position d'un point dans l'espace physique à la position qu'il occupe dans les images considérées. Cela offre ensuite la possibilité de retrouver le déplacement d'un point par la mise en correspondance des images. Nous détaillons l'implémentation, réalisée au cours de cette thèse, du procédé de calibration photogrammétrique et utilisée dans l'application pratique finale de scan multi-vues. Ensuite, les méthodes CIN et SCIN classiques sont progressivement exposées. Cela permet de faire apparaître dans le détail les deux limitations principales identifiées.

La première d'entre elles est liée au fait que les structures complexes sont souvent associées à des états de déformations hétérogènes. Il en résulte que des instrumentations multi-échelles, c'est-à-dire des caméras qui observent la région d'intérêt à différents niveaux de résolution, deviennent pertinentes. Cela permet en effet d'obtenir à la fois des informations locales précises dans les zones de fortes déformations et des informations globales sur la géométrie et les conditions aux limites. Cependant, nous montrons que les mouchetis usuels peuvent compromettre l'utilisation de méthodes multi-échelles car ils ne sont bien souvent adaptés qu'à une seule échelle (ou, dans le meilleur des cas, à deux échelles). En outre, parvenir à fusionner les données issues d'images prises à différentes échelles en une unique mesure de manière adéquate constitue une question ouverte dans le formalisme actuel. En règle générale, les mesures sont réalisées à chaque échelle indépendamment, puis, pour chaque point, seule la mesure la plus résolue est conservée. En pareil cas, une grande quantité d'information est tout simplement ignorée, ce qui est fâcheux étant donné le peu de tests réalisés aux grandes échelles.

La seconde limitation principale identifiée est liée aux plus grands déplacements attendus, à déformation donnée, sur de grandes structures. Cela peut en effet se traduire en des hypothèses (des approches usuelles) devenant de plus en plus dures à satisfaire à mesure que l'échelle augmente. Par exemple, garantir qu'un point de mesure reste visible à tout moment par une même caméra restreint grandement le domaine d'utilisation de la CIN. Un autre aspect des grands déplacements est lié à d'importantes variations de la normale à la région d'intérêt au cours du temps. Dépendant de la distribution de lumière incidente, cela peut se traduire en des variations de niveaux de gris dans les images d'un même point au cours du temps, et, plus important encore, en une modification de l'échantillonnage de la surface réalisé par les caméras et donc du niveau de confiance que l'on peut attribuer à cette observation au cours du temps.

Travailler sur la formulation du problème de CIN dans un contexte multi-échelle était donc un objectif à long terme de cette thèse. Mais avant de s'y employer, nous nous sommes d'abord assurés que ces éventuels développements théoriques n'étaient pas vains en définissant des mouchetis multi-échelles capable d'opérer sur une large gamme d'échelles. Pour l'un d'entre eux, nous montrons que ses propriétés statistiques intrinsèques d'invariance par changement d'échelle se traduisent en une erreur de mesure indépendante de l'échelle considérée. Cela n'est pas le cas des mouchetis usuels.

Ensuite, par une modélisation physique rigoureuse et grâce à des hypothèses explicites, une formulation CIN générale est proposée. Elle est basée sur des outils issus de la vision par ordinateur, une approche photométrique, raison pour laquelle nous l'appelons CIN photométrique (PhDIC, de l'anglais Photometric Digital Image Correlation), et la définition d'un jumeau numérique des essais. Ce jumeau numérique comprend un modèle géométrique raffiné de la structure incluant son albedo (ou une paramétrisation de la fonction de réflectivité bidirectionnelle, dans une approche plus générale), et un modèle d'éclairage. Cela offre l'avantage de naturellement prendre en compte les problèmes de visibilité, pondération des observations entre elles (échelles, échantillonnage en pixel des surfaces) et donc de l'assimilation de données. Un autre bénéfice est la possibilité de développer de nouvelles formulations basées sur ce travail par la remise en question du bien-fondé de certaines hypothèses, ou autrement dit la capacité à modéliser des phénomènes physiques.

L'emploi du terme "général" pour qualifier la formulation proposée est justifié par la démonstration du fait que, non seulement, elle englobe les approches usuelles mais qu'elle permet aussi, en particulier, de réaliser des mesures quand les images de référence et déformée représentent des zones complétement disjointes, grâce aux termes d'intercorrélation spatio-temporels, et prévient la dérive temporelle, ou perte de stéréo-correspondance. En plus, en tant que pendant dense de l'ajustement de faisceaux ('Bundle Adjustment'), elle permettrait de réduire les incertitudes de mesure associées aux approches SCIN en ne considérant qu'un seul problème de minimisation au lieu de plusieurs problèmes les uns à la suite des autres où les incertitudes s'accumulent tout au long de la chaîne. Grâce à une analyse de sensibilité, basée uniquement sur la mesure en déplacement, nous montrons qu'une réduction importante de l'incertitude de mesure est attendue du simple fait de se baser sur un jumeau numérique (ou modèle). La possibilité de différencier ce modèle au lieu des images est également examinée, car les analyses habituelles évacuent ce problème. Enfin, le choix, dans la définition de la fonctionnelle, d'intégrer en premier lieu dans les images est justifié du point de vue de la théorie de l'estimation statistique. Cela permet de discuter des conditions de validité du caractère optimal de cette définition (par exemple, un modèle qui explique parfaitement les données au bruit près ou encore une distribution de bruit Gaussienne).

Un cas test suivant les étapes de mesure de forme et de déplacement sur une éprouvette trouée est proposé. Les résultats concernant la forme, l'albedo, la mesure de déplacement sont présentés et critiqués. Cela permet de montrer, en pratique, que les problèmes de lumière et visibilité sont naturellement traités par la formulation proposée. Il est notamment démontré que cette méthodologie permet de mesurer un déplacement (de corps rigide, en l'occurrence) mettant complètement en défaut les approches SCIN conventionnelles.

De nombreuses perspectives sont suggérées à ce travail. Il est question de cas d'essais qui permettraient de mettre davantage en valeur la formulation, par un intérêt industriel et une difficulté à instrumenter avec les méthodes classiques de SCIN. On peut penser, par exemple, à de grandes transformations avec de grandes déformations nécessitant des traitements particuliers dans les approches conventionnelles du fait que l'image de référence n'a plus rien à voir avec l'image déformée (essais d'impact de pneus sur un panneau aéronautique, traction d'élastomère). Une autre possibilité est celle des grandes transformations qui induisent une perte de visibilité (modes de flambement).

De toute évidence, la PhDIC est basée sur des principes radiométriques, de la même manière que la Thermographie Infra-Rouge (TIR). Coupler les deux méthodes de mesure en une seule et même approche pourrait être intéressant. Habituellement, la TIR suppose une surface d'émissivité uniforme, alors que la CIN requiert un mouchetis pour produire des mesures fiables. Ainsi, déposer un mouchetis dans le spectre infra-rouge, puis caractériser l'émissivité en chaque point de mesure, de la même manière que nous avons mesuré l'albedo dans notre dernière application pratique, pourrait permettre de mesurer à la fois un champ de déplacements et un champ thermique, à l'aide d'une seule caméra infra-rouge.

Aussi, nous avons montré dans ce travail que la PhDIC permettait le repositionnement de caméras en cours d'essais. Dès lors, des essais où tout ou partie des caméras sont montées sur des drones et/ou des bras robotiques semble être la prochaine étape à franchir pour la SCIN. Cela permettrait, couplé à des processus multi-échelles, de balayer les surfaces d'essais, avant les essais, pour les scanner et mesurer leurs formes et propriétés de réflectance. Au cours des essais, certaines caméras pourraient continuer à parcourir les surfaces à la recherche de points chauds sur lesquels se focaliser ensuite. Cependant, la phase de scan préalable ne semble pas forcément être des plus pratiques pour l'expérimentateur, et utiliser différentes positions de l'éclairage à la place, tout en conservant les caméras aux mêmes positions, à la manière du Shape-from-Shading, pourrait s'avérer plus pertinent dépendant de la configuration d'essai.

Concernant le jumeau numérique, beaucoup reste à faire. Comme expliqué précédemment, dans le contexte des grandes transformations, les problèmes associés à la lumière peuvent devenir difficiles à traiter, il pourrait donc être intéressant d'intégrer une paramétrisation de la fonction de réflectivité bidirectionnelle incluant non seulement une composante diffuse (ou Lambertienne), mais également une composante spéculaire. D'autres modèles d'éclairage, que celui de la lumière située à l'infini, pourraient également être utilisés. L'objectif ultime serait de pouvoir trouver des modèles de réflectance et de lumière adaptés à une large gamme d'essais aux grandes échelles. Cependant, cela ne constitue pas l'unique manière de procéder. En conservant un modèle Lambertien malgré des problèmes liés à l'éclairage (réflexions spéculaires par exemple), les hypothèses liées au caractère optimal de la formulation deviennent caduques car le modèle n'explique plus les données au bruit près. Il peut donc devenir pertinent d'étudier l'utilisation d'autres normes de résidus moins sensibles aux données aberrantes.

General Multiview Stereo Digital Image Correlation Formulation: towards multiscale digital twins of complex experimental setups Keywords: Data assimilation, Photometric Digital Image Correlation, Fractal pattern, Structure scan, Large Deformations, Physically based model Abstract: To meet the future aircraft energy efficiency requirements, breakthrough innovations will be needed, including in structural mechanics. Currently, aircraft design cycles rely on a large number of tests performed at the smallest scales. At the same time, measurement point densities tend to be smaller at the largest scales because of costs and time associated with setting up sensors. This results in a lack of predictability of the associated models. Based on these observations, extending the operability of full-field measurement techniques to larger scales stands for a relevant way to increase the amount of available data and thus address the aforementioned predictability issue.

In this work, we present the main challenges associated with this upscaling to complex structures when considering the Stereo Digital Image Correlation (SDIC) measurement technique. SDIC is a popular, contactfree, optical measurement method. It relies on the matching in pictures taken by cameras of features contained in patterns deposited on possibly non-planar test sample surfaces. This allows then to measure the three-dimensional displacement field on these surfaces.

The first identified issue when considering larger scales is the need for cameras observing tests at different levels of resolution. For instance, it is clear that, for efficiency purposes, regions associated with high strain gradients should be more resolved than those with lower ones. However, patterns are usually designed so as to reach an optimum value of 3 to 5 pixels per speckle. This guideline is not relevant in a multiscale context, where cameras with different pixel sizes are used. The use of fractal (self-affine) patterns with intrinsic multiscale properties is thus investigated. It is shown that for one of the proposed patterns and two different displacement fields, the associated scale invariances result in measurement uncertainties independent of the scale at which the pattern is considered, which is not the case for usual patterns.

Then, when considering complex structures, questions are raised such as lighting or visibility which are often addressed on a case-by-case basis in current SDIC frameworks. Here, we draw inspiration from Computer Vision developments and propose a general SDIC formulation based on a physical modelling of the image formation process and photometry. For this reason we call it Photometric DIC (PhDIC). It addresses all of the aforementioned issues at once, by relying on a digital twin of the experimental setup. In addition, this digital twin allows the generation of virtual pictures to compare directly actual data with. Thus, PhDIC is not based on the comparison between reference pictures and deformed ones. This offers the possibility to intrinsically account for complex deformations that totally incapacitates usual DIC frameworks.

A theoretical justification of the PhDIC formulation is proposed. Its general features are exhibited by showing that this formulation encompasses the usual SDIC frameworks, while being able to operate over a wider range of cases. Also, a sensitivity analysis demonstrates that measurement uncertainties are expected to be reduced compared to usual SDIC. Finally, based on estimation theory, we make explicit the assumptions under which the formulation is optimal.

A practical test case on 360-degree multiview pictures is investigated and a large rotation which incapacitates usual approaches is measured over the considered surface.

Finally, a broader vision is adopted, and unprecedented perspectives made possible by PhDIC are suggested.

Formulation Générale de la Stéréo Corrélation d'Images Numériques Multi-vues:

vers des jumeaux numériques multi-échelles d'essais complexes.

Mots-clés: Assimilation de données, Corrélation Photométrique d'Images Numériques, Mouchetis fractal, Scan de structure, Grandes transformations, Modèle basé sur la physique Résumé: Afin d'atteindre les futurs objectifs d'efficacité énergétique du secteur aéronautique, des innovations de rupture seront nécessaires, y compris en mécanique des structures. Actuellement, les cycles de conception aéronautiques reposent sur un grand nombre d'essais aux plus basses échelles. Aussi, les densités de mesure ont tendance à être plus faibles aux grandes échelles à cause des coûts et du temps associés à la mise en place des capteurs. Cela résulte en une faible prédictibilité des modèles associés. Partant de ce constat, il paraît judicieux d'étendre le domaine d'application des méthodes de mesure de champs à des échelles plus importantes afin d'y accroître la quantité de données disponibles et de s'attaquer ainsi à ces problèmes de prédictibilité. Dans ce travail, nous présentons les principaux verrous liés à la mesure aux grandes échelles sur structure complexe par des méthodes de Stéréo Corrélation d'Images Numériques (SCIN). La SCIN est une méthode de mesure optique, peu coûteuse. Elle repose sur la mise en correspondance d'images de la surface de mesure, éventuellement non plane, prises par différentes caméras, à différents instants. Cette mise en correspondance, facilitée par le dépôt d'un mouchetis, permet la mesure d'un champ de déplacement tridimensionnel.

Aux grandes échelles, le premier verrou identifié est lié au besoin d'avoir des caméras dont les résolutions sont différentes. Par exemple, pour des raisons d'efficacité, les régions aux grands gradients de déformations doivent être mieux résolues que celles où ils sont plus faibles. Cependant, les mouchetis sont généralement déposés de sorte que la valeur optimale de 3 à 5 pixels par tache soit atteinte. Cette règle perd de son sens dans un contexte multi-échelle, puisque la taille des taches dépend de la caméra. Dans ce contexte, l'utilisation de mouchetis fractals (auto-affines) aux propriétés intrinsèquement multi-échelles est étudiée. Nous montrons, pour un tel mouchetis et deux champs de déplacement différents, que les invariances par changement d'échelles associées se traduisent en des incertitudes de mesure indépendantes de l'échelle à laquelle le mouchetis est considéré, ce qui n'est pas le cas des mouchetis usuels.

Ensuite, lorsque l'on s'intéresse à des structures complexes, des problématiques telles que l'éclairage ou la visibilité émergent. Elles sont habituellement traitées au cas par cas. Dans ce travail, nous nous inspirons de travaux issus de la Vision par Ordinateur et proposons une formulation SCIN générale basée sur une approche photométrique. Nous l'appelons donc CIN Photométrique (PhDIC, de l'anglais Photometric Digital Image Correlation). Elle répond à toutes les problématiques précédentes à la fois, en s'appuyant sur un jumeau numérique de l'essai. La PhDIC n'est donc pas basée sur la comparaison d'images de référence à des images déformées mais d'images réelles à un jumeau numérique. Cela offre la possibilité d'intrinsèquement prendre en compte des transformations complexes mettant en échec les méthodes conventionnelles.

Une justification théorique de la formulation PhDIC est proposée. Son caractère général est justifié en montrant que cette formulation englobe les approches SCIN classiques, tout en étant capable d'opérer sur une plus large gamme d'essais. Aussi, une analyse de sensibilité établit qu'une réduction des incertitudes de mesure par rapport à la SCIN classique est attendue. Enfin, à partir de la théorie de l'estimation statistique, nous justifions l'optimalité de cette formulation, tout en explicitant les hypothèses nécessaires.

Un cas d'application sur des images multi-vues à 360 degrés est étudié. Une grande rotation mettant en défaut les méthodes classiques y est mesurée.

Enfin, une vision plus large est adoptée, et des perspectives inédites rendues possibles par la PhDIC sont suggérées.

  tions (θ f , ϕ f ) (falling) and (θ e , ϕ e ) (emitted) ρ Albedo Î Substitute image (= αL 0 f ρ in the case of an ambient lighting of uniform radiance L Jacobian matrix of P (absolute value) V Visibility function D Discrepancy map (indifferently shape correction or displacement) ρ Albedo in the deformed state

  Annual Report of the Council Air traffic has been vulnerable to external factors including disease outbreaks

Figure 1 :

 1 Figure 1: Revenue Passenger-Kilometres (thick line) and Freight Tonne-Kilometres (thin line). Reproduced from [OACI 2020].

Figure 3 :

 3 Figure 3: Dynamic strain sensors on the ground vibration testing of an F-16 aircraft, illustrating the coarse measurement point density on a large-scale test. Extracted from [dos Santos et al. 2015].

Figure 4 :

 4 Figure 4: Falcon 7X fuselage local buckling during static tests. Figure reproduced from [Barrière 2014].

Figure 5 :

 5 Figure 5: Illustration of the test-simulation dialogue allowed by DIC in the context of an impact on a panel. Images are courtesy of DGA Aeronautical Systems.

Figure 1 . 1 :

 11 Figure 1.1: Direct problem.

  Figure 1.2: Inverse problems classification illustration.

Figure 2 . 1 :

 21 Figure 2.1: Schematic drawing of the different coordinate systems introduced in the pinhole model.

Figure 2

 2 Figure 2.3: "Block diagram" of the pinhole model.

Figure 2

 2 Figure 2.4: Definition of the skew angle, accounting for the non-orthogonality of pixel lines.

  Figure 2.6 sums up previous notations and steps.

  Figure 2.5: Pictures of a chessboard without distortions (a), with δθ = 0 (b), r 1 = 0 (c), r 2 = 0 (d), r 3 = 0 (e), d 1 = 0 (f), d 2 = 0 (g), p 1 = 0 (h), p 2 = 0 (i).

Figure 2 . 7 :

 27 Figure 2.7: Schematic drawing and notations for a stereo rig

Figure 2 . 8 :

 28 Figure 2.8: Target picture and detected target reference frame axis.

Figure 2 . 9 :

 29 Figure 2.9: Initialisation for camera 0. On the Figure are shown the 16 pictures that have been taken by camera 0 of the target. We superimposed the evaluated positions of the target points on the picture as the yellow dots. σ 0 = 1.22 pixels (error standard deviation).

Figure 2 .

 2 Figure 2.10: Initialisation for camera 1. On the Figure are shown the 16 pictures that have been taken by camera 1 of the target. We superimposed the evaluated positions of the target points on the picture as the yellow dots. σ 1 = 4.07 pixels (error standard deviation).

  Figure 2.11: σ ci (error standard deviation in pixels for picture i, camera c) at two different stage of the calibration process, before and after optimisation. Error standard deviations are also presented for two independent optimisation processes (one for each camera) as described in Section 2.2.2 and for a single optimisation process as described in Section 2.2.3.

Figure 3

 3 Figure 3.1: Left, reference state image f and considered point x. Right, deformed state image g and identified displacement u(x) associated to the point x. Pattern from [Fouque et al. 2021c].

Figure 3 . 2 :

 32 Figure 3.2: Displacement tracking strategy of local DIC (top) and global DIC (bottom). Extracted from [Wang & Pan 2016].

Figure 3 . 5 :

 35 Figure 3.5: Two points M 1 and M 2 lying on a line passing through camera centre O c have the same projection m 1 . The use of another camera allows to distinguish M 1 from M 2 thanks to the two different projections m 1 and m 2 . Inspired from [Garcia 2001, Figure 1.10].

Figure 3

 3 Figure3.6: Subset SDIC framework: pixel subsets are defined in the master (left) reference state (t 0 ) image f L and their coordinates are sought in the slave (right) reference state image (f R ) based on standard two-dimensional DIC processes. The knowledge of the camera relative positions allows then to triangulate the position of each subset centre in the world reference frame R w and obtain the shape S of the specimen. In the deformed state (t 0 + ∆t), the procedure is the same except that subsets remain defined in the master reference state image (f L ). A temporal matching is thus performed with the left deformed state image (g L ) and with the right one (g R ). Again, a triangulation allows to obtain the position of subset centres in the deformed state (S ). The displacement field is retrieved by subtraction: U = S -S. Modified from[Pierré et al. 2017].

Figure 3 . 7 :

 37 Figure 3.7: Three major strategies for the stereo-correlations in stereo-DIC: (a) stereo matching for initial left and right images, and temporal matching for subsequent left and right image series; (b) temporal matching for left image series and stereo matching for all the right images using the original left image as the reference image (same strategy as Figure 3.6), and; (c) temporal matching for the left images and stereo matching for all the right images, using the current left image as the reference image. Extracted from [Pan 2018].

Figure 3 . 9 :

 39 Figure 3.9: Illustration of the shape measurement ill-posedness. (a) Theoretical position of the mesh after extrinsics calibration and shape measurement. (b) Global or (c) Local sliding during the shape measurement. Extracted from [Pierré et al. 2017].

  Figure 3.10: Diagrams introducing in a graphical manner the current Subset-based and Finite-Element Multiview Stereo Digital Image Correlation frameworks.

Figure 4

 4 Figure 4.1: N s -level multigrid or pyramidal scheme with N s = 3.

Figure 4

 4 Figure 4.2: Principle of the geometric regularisation: the multilevel spline optimisation process is coupled with the multiscale initialisation of shape measurement. Extracted from [Colantonio et al. 2020].

  (i) S,f from Equation (3.14) are respectively replaced by M S + w R F R Γ and b

  (a) Biscale optimised pattern. Extracted from [Bomarito et al. 2018]. (b) Biscale multispectral pattern. Extracted from [Mathew et al. 2018].

Figure 4

 4 Figure 4.4: Biscale patterns defined for multiscale DIC.

Figure 4

 4 Figure 4.5: Deformation test example. Lines of each color represent a cross section of the correlation surface at a given stretch. Extracted from [Bomarito et al. 2017].

Figure 4 .

 4 Figure 4.6 shows slices of different auto-correlation functions. The solid blue line represents the desired auto-correlation function A of Equation (4.14) for n = 127, H = 1/22 and τ x = 0.

Figure 4 . 6 :

 46 Figure 4.6: Comparison between desired (solid blue line, A) and actual (dashed orange line, A h ) auto-correlation functions for n = 127, H = 1/22 and τ x = 0. The green dash-dotted line is the auto-correlation function of F µ,0.3σ gl (h), a MIG-improved version of h.

  , 4.8i, 4.8k and 4.8l. However, for our correlation-length-free pattern, the subsampled picture dynamic range was F µ,σ (h), σ = σ gl , MIG= 58.0. F µ,σ (h), σ = 0.3σ gl , MIG= 84.6. Randomly generated pattern (each pixel can be equal to either 0 or 255 with a 50% chance), MIG= 109.9.

Figure 4 . 7 :

 47 Figure 4.7: Grey level distributions of h and F µ,σ (h) for different values of σ and grey level distribution of a randomly generated pattern. MIG is also given for each pattern. Histograms and MIG values are given for n = 127.

  Nn k=1 α∈{x,y}(u m αk -u imp αk -e u ) 2 , (4.17) F µ,0.3σ gl (h), n = 8191.

Figure 4 . 8 :

 48 Figure4.8: The first column is for a pattern generated using the proposed technique, F µ,0.3σ gl (h). The second one shows a randomly generated pattern where each pixel has a one-half chance of being equal to 0 or 255. The last column concerns a pattern based on a Perlin noise with a speckle size of approximately 3 pixels. For each column, the first figure is a picture of the pattern considered for n = 8191. The second picture is a zoom on a 100×100-pixel square of the 16384×16384 pattern. The third one is a subsampling of the pattern considered (each grey level value is equal to the integer part of the mean over the corresponding 512×512-pixel square). The last picture is the grey level distribution of the subsampling. Colour scale for pictures goes from 0 (black) to 255 (white).

Figure 4

 4 Figure 4.9: Mesh used for the FE-DIC measurements superimposed on h (n = 127).

  in pixel (b) Random error expectancy (in pixel) for subpixel translations.

Figure 4 .

 4 Figure 4.10: Systematic and random error expectancies for sub-pixel translations. Results are presented for the random pattern, the initial pattern h, F µ,σ (h) for five different values of σ ranging from 0.2σ gl to σ gl , and the Perlin noise pattern.

Figure 4 .

 4 Figure 4.11: Systematic error expectancy absolute value and random error expectancy for stretches. Results are presented for F µ,0.3σ gl (h), the random and the Perlin noise patterns.

  Figure4.12 shows corresponding results with respect to the scale s for the three patterns. For scales 5 and 6, the Perlin noise pattern could not converge, even for the smallest stretch strain value (ε = 0.01). This can be seen in Figure4.12e, where the value plotted for the convergence robustness is 0. Strictly speaking, for these scales, the convergence robustness of the Perlin noise pattern is not 0 but is less than 0.01. As a result, and as shown in Figures 4.12c and 4.12d, the Perlin noise pattern was discarded for scales 5 and 6.Relatively to the random and the Perlin noise patterns, and for the range of scales considered, the DIC errors associated with the pattern proposed in this work do not depend on the scale s. Except for the steep increase of the mean of the systematic error expectancy absolute value between scales 0 and 1 for translations (see Figure4.12a), our pattern shows stable error values. Regarding the convergence robustness in Figure4.12e, this output tends to improve as the scale increases for the proposed pattern while the opposite trend can be observed for the other two patterns.Figures 4.12c and 4.12d show the mean of the systematic error expectancy absolute value and of the random error expectancy for stretches for the three patterns considered. For each

Contents 5 . 1

 51 From the light source to the grey level value . . . . . . . . . . . . 89 5.1.1 A general overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.1.2 Lambertian model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.1.3 Lighting models . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5 . 2 :

 52 Figure 5.2: Bidirectional Reflectance Distribution Function (BRDF) parameterisation.

  Figure 5.3: Diagram introducing different applications and notations used herein.

Contents 6 . 1

 61 Link between PhDIC and usual multiview frameworks . . . . . . 102 6.1.1 State of the Art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.1.2 Shape measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.1.3 Displacement measurements . . . . . . . . . . . . . . . . . . . . . . . 106 6.1.4 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 6.3.1 Maximum likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.3.2 Prior models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.3.3 Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 7 . 1 :

 71 Figure 7.1: Illustration of the practical use made of the Photometric Stereo Digital Image Correlation (PhDIC) in this chapter.

  (a) Example of an input image used for the shape and albedo scan. Superimposed blue dots stand for the light calibration points. The world reference frame is shown in orange. Both were not present in the original picture. (b) Image positions relatively to the coupon. The mesh has been enlarged for reasons of visibility and thus sizes and distances are not to scale.

Figure 7

 7 Figure 7.2: Calibration, shape and albedo measurement setup.

Figure 7

 7 Figure 7.4: Structure of the algorithm used for the functional alternating optimisations.

  Figure 7.5: Comparison between initialisation and convergence states of the specimen shape. The mesh is superimposed on pictures.

Figure 7 . 6 :

 76 Figure 7.6: Target point reprojection error standard deviation for each picture of one of the considered image sets. The total reprojection error standard deviation across all pictures is equal to 0.18 pixel.

  Figure7.7: Distribution of the normalised albedo difference (defined at integration points) for the two different sets of pictures at convergence state. The mean value is small compared to the standard deviation. It should be noted that camera noise is of the same order of magnitude as the albedo difference standard deviation.

  Distribution of the integration points coordinates difference projected along X-direction at convergence state. Distribution of the integration points coordinates difference projected along Y-direction at convergence state. Distribution of the integration points coordinates difference projected along Z-direction at convergence state.

  (a) One of the multiview pictures used for the extrinsics, shape and albedo calibration. Above is shown the position of all these multiview pictures with respect to the beam. In a usual SDIC framework, this picture would stand for the reference state image I 0 i for displacement measurement.(b) One of the multiview pictures taken for the displacement measurement of a 90 • rotation. We chose this picture (instead of the actual picture used) for illustration purposes: it is really close to Figure7.10a. The actual position of the only picture used to perform the displacement measurement is shown above. The mesh superimposed on the picture is the initialisation used (94 • rotation). In a usual SDIC framework, this picture would stand for the deformed state image I i corresponding to I 0 i in Figure7.10a.(c) Position of the mesh superimposed on a deformed state image at convergence state.

Figure 7 .

 7 Figure 7.11: Specular reflection appearing on SDIC images of an airplane tyre sidewall during experiments mimicking the tyre/runway interaction. Courtesy of DGA Aeronautical Systems.

  x + νε y,x + (1 -ν)ε xy,y (1 -ν)ε xy,x + ε y,y + νε x,y 0

  5)We have two equations but three unknowns (X c,u , Y c,u and Z c,u ). To obtain a third equation, we can derive Equation (B.4) with respect to u: be rewritten n, O w M ,u = 0 with n = ∇g/ ∇g 2 . Also we know thatO w M = O w O c + O c M . Making use of both yields n, O c M ,u = 0 which is the third equation. By decomposing O c M ,u as X c,u e c x + Y c,u e c y + Z c,u e c z with (e cx , e c y , e c z ) the orthonormal basis associated with the camera reference frame, we obtain:n, O c M ,u = n, e c x X c,u + n, e c y Y c,u + n,e c z Z c,u = 0. (B.7) By using Equation (B.5) in Equation (B.7), we can get an expression of Z c,u . And Z c,u allows then to retrieve X c,u and Y c,u thanks to Equation (B.5):

  y n, O c M -n, e c x O c M × e c y -n, e c y e c x × O c M + n, O c M e c M -X c n, e c x -Y c n, e c y y n, O c M n (B.10) Since, by construction, n 2 = 1, we finally obtain:

  identified contribution of the weighting term J ( f Zc , | n, r | and OcM 2 Zc ) in Figure B.1.

  Figure B.1b. Similarly, in Figure B.1c this ratio equals (sin α) -1 = OcM 2 Zc . [Zhou et al. 2015] Yihao Zhou, Chen Sun, Yuntao Song and Jubing Chen. Image pre-filtering for measurement error reduction in digital image correlation. Optics and Lasers in Engineering, vol. 65, pages 46 -56, 2015. Special Issue on Digital Image Correlation doi:10.1016/j.optlaseng.2014.04.018. (Cited on pages 47 and 85.) Abstracts Formulation Générale de la Stéréo Corrélation d'Images Numériques Multi-vues: vers des jumeaux numériques multi-échelles d'essais complexes.

field measurements in Digital Image Correlation

  ext Extrinsics (t x t y t z α β γ) T p int Intrinsics (f x f y u 0 v 0 ) T d Distortions (δθ r 1 r 2 r 3 d 1 d 2 p 1 p 2 ) T

	x		Notations
			
	p	Camera model parameters	 
			 
	W	Physical space
	P	Camera parameter space
	Π	Image plane
	I	Domain of Π corresponding to the picture
	P	Camera projection model
	1	Identity matrix
	N p	Number of points (calibration target)
	N i	Number of images
	N c	Number of cameras
	xij	Detected pixel position of point j in picture i
	θ	Unknown vector associated with a single camera calibration
	Θ	Unknown vector associated with a multiple camera calibration
	Full-f	Reference state image
	g	Deformed state image
	u	Two-dimensional displacement vector
	U	Three-dimensional displacement vector
	N (i) S	2D Digital Image Correlation Gauss-Newton algorithm Hessian matrix
		at iteration i, for the functional defined over S
	b	Right-hand side in Gauss-Newton minimisation schemes
	f	Mean value of f
	∆f	Standard deviation of f
	I t c	Image taken at time t by camera c
	N t	Number of time steps
	S	Three-dimensional shape correction field
	V h	Finite-dimensional subspace for the Ritz-Galerkin method
	H	Stereo Digital Image Correlation Hessian matrix
	R shape	Reduction matrix for the shape measurement step
		su sv s) T
	R	Rotation matrix
	T	Homogeneous rigid body matrix

p x

Pixel coordinates of a point M (x = m pix ) X World coordinates of a point M (X = M w ) j Shape function associated with the jth degree of freedom N Matrix collecting the shape functions, such that, for instance, u(x) = N (x) u M

Notations xi Multiscale approaches in Digital Image Correlation

	s	Scale in multigrid or pyramidal schemes
	N s	Number of levels of the multigrid or pyramidal scheme
	h	Discrete quantised pattern
	n	2n + 2 is the number of pixels in each direction of h
	N gl	Number of quantisation levels (grey levels)
	τ	Shift (or translation) vector (τ x τ y )
	A h	h auto-correlation
	F T	Fourier Transform
	F T -1	Inverse Fourier Transform
	H	Hurst exponent
	ε	Stretch strain value

Contents 2.1 Camera modelling

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.1 The pinhole model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.2 Distortion models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2 Photogrammetric calibration process . . . . . . . . . . . . . . . . 21 2.2.1 Parameters initialisation . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.2 Process for a single camera . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.3 Multiple-camera setup . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.3 Example on an actual stereo rig . . . . . . . . . . . . . . . . . . . . 30 2.3.1 Initialisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Contents 3.1 Two-dimensional Digital Image Correlation . . . . . . . . . . . . . 38

	3.1.1 Grey level conservation or strong formulation . . . . . . . . . . . . . 38
	3.1.2 Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . 39
	3

.2 Ritz-Galerkin method . . . . . . . . . . . . . . . . . . . . . . . . . . 41

	3.2.1 General formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
	3.2.2 Local approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
	3.2.3 Global approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.3 Including camera models . . . . . . . . . . . . . . . . . . . . . . . . 43 3.4 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . 46

	3.4.1 Integration schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
	3.4.2 Subpixel image evaluation and image gradient . . . . . . . . . . . . . 47
	3.4.3 Correlation criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
	3.5

Extension to Stereo Digital Image Correlation . . . . . . . . . . . 49

	3.5.1 Subset-based Stereo Digital Image Correlation . . . . . . . . . . . . 49
	3.5.2 Finite-Element-based Stereo Digital Image Correlation . . . . . . . . 50
	3.6

Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

	3.6.1 Stereo Digital Image Correlation and large-scale applications . . . . 58
	3.6.2 Deposited patterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.1 Pyramidal approach for initialisation . . . . . . . . . . . . . . . . . 62 4

.1.1 Coarse-graining process . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.1.2 Regularisation methods in Digital Image Correlation . . . . . . . . . 63 4.

2 Fractal pattern for Multiscale Digital Image Correlation . . . . . 66 4

.2.1 Pattern quality criteria . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.2.2 Multiscale pattern generation from auto-correlation function . . . . 70 4.2.3 Pattern performance evaluation . . . . . . . . . . . . . . . . . . . . . 74 4.3 Conclusion and outlook . . . . . . . . . . . . . . . . . . . . . . . . . 84

Table 4

 4

	0.03		Systematic error expectancy in pixel
	0.02
	0.01
	0.00
	-0.01
	-0.02
	-0.03	0.0	0.2	0.4 Translation value 0.6	0.8	1.0
	(a) Systematic error expectancy (in pixel) for sub-
	pixel translations.

.1 sums up the performance achieved by F µ,0.3σ gl (h) with respect to the random

  Table 4.1: Comparison of the pattern F µ,0.3σ gl (h) with the random and the Perlin noise patterns. Each line represents the mean over the 21 translations of Figure 4.10, for systematic error expectancy absolute value and for random error expectancy.

	Pattern	F µ,0.3σ gl (h) Random pattern	Perlin noise (3 pix/speckle)
	Mean of the systematic
	error expectancy	4.19 × 10 -4	8.60 × 10 -3	1.89 × 10 -3
	absolute value (pixel)
	Mean of the random error expectancy (pixel)	1.29 × 10 -2	1.59 × 10 -2	1.27 × 10 -2
			More precisely, the relation

Table 4 .

 4 2: Comparison of the pattern F µ,0.3σ gl (h) with the random and the Perlin noise patterns. Each line reports the mean over the first 9 measurements of Figure4.11 (ε ∈ {0.01, 0.02...0.09}), for either systematic (upper line) or random (lower line) error expectancies.Multiscale displacement-field measurements This section aims to demonstrate the suitability of one of the proposed patterns in a multiscale DIC setup in terms of expected errors on displacement field measurements. The process used for this purpose was similar to that described in Section 4.2.3 and thus generated a 16384 × 16384-pixel picture for each of the three patterns considered (i.e. F µ,0.3σ gl (h) for n = 8191, the random pattern and the Perlin noise pattern with 3-pixel-wide speckles, see Figures 4.8a to 4.8c). For each pattern and for every scale s ∈ [[0; 6]], we generated a 16384 2 s × 16384 2 s -pixel image via s steps of the coarse graining process used in Section 4.2.3. Then from this aggregated picture only a 256 × 256-pixel portion located in the centre was considered and extracted. Finally, exactly the same displacement measurements as in Section 4.2.3 were performed, the only difference being that, here, ε could reach 25% (instead of 20% in Section 4.2.3).

  Table 4.3: F µ,0.3σ gl (h) results for stretches and the scales considered. For each scale, means are computed over the whole range of convergence of this pattern.

	Scale s	Mean of the systematic error expectancy absolute value (pixel)	Mean of the random error expectancy (pixel)
	0	4.98 × 10 -4	1.67 × 10 -2
	1	2.26 × 10 -4	1.32 × 10 -2
	2	3.05 × 10 -4	1.48 × 10 -2
	3	4.80 × 10 -4	1.64 × 10 -2
	4	3.44 × 10 -4	1.49 × 10 -2
	5	2.74 × 10 -4	1.60 × 10 -2
	6	3.48 × 10 -4	1.40 × 10 -2

Table 4

 4

.3 would not be so close to the values shown in Figures 4.12c and 4.12d for F µ,0.3σ gl (h).

Remark Error values for F µ,0.3σ gl (h) at the scale s = 0 differ quite a bit from the values of Section 4.2.3 (see Tables 4.2 and 4.3 for instance). This can be explained by the fact that the patterns considered are not exactly the same. In Section 4.2.3, the pattern was generated with n = 127 to produce a 256 × 256-pixel picture. In this section it was generated with n = 8191 and then a 256 × 256-pixel picture was extracted.

  Ωt c vis denotes the visible part of Ωt in the picture taken by camera c at time t such that there is a one-to-one relation between Ωt c vis and S t c thanks to the projection map P t c .

	.10)
	where Π t c stands for the image plane. Note that Π t c depends on t, as we adopt a general
	formulation where camera positions may change. Some of the introduced notations are presented
	in Figure 5.3 for better understanding. P t c is the camera model associated to Π t c and depends on camera parameters p t c . The notation Ωt is used to refer to the deformed state of the ROI Ω at time t. We also define the silhouettes S t c = P t c ( Ωt ) ∩ I t c = P t c ( Ωt c vis ) where I t c (⊂ Π t c ) stands
	for the bounded domain of Π t c corresponding to the image. Remark Strictly speaking, in Equation (5.10), I t c is defined over I t c . Yet, as images are usually
	interpolated, this makes it possible to define I t c over the whole image plane Π t c .

Acknowledgements

Chapter 7

Application on a real test case

We now propose a practical application of PhDIC on real images of an open-hole specimen. Only two time steps are considered, the reference and the deformed ones. As a first proof of concept, a framework slightly downgraded compared to the one of Chapter 5 is proposed, as, for reasons that will be explained later on, the whole Functional (5.18) is not minimised. Instead, it is minimised in the reference state to calibrate the extrinsics, the shape and the reflectance model (in the same way as the shape measurement step). Then the part of the functional associated with the deformed state is minimised considering only the displacement as a variable (in the same way as the displacement measurement step in usual frameworks). The associated results are presented successively.

Pinhole model functional weighting term

This appendix originates from [Fouque et al. 2021a].

We have the identity

. However, the backprojection operator is very costly to evaluate for a camera model accounting for distortions, and its gradient even more. Instead, we propose an approximation of J assuming a pinhole camera model. In this case we can write:

where u 0 and v 0 denote the optical centre pixel coordinates, f x and f y the product of the focal length and the camera sampling parameter along each direction, (X c , Y c , Z c ) the coordinates of a point M in the camera reference frame (denoted R c ), and u and v the pixel coordinates of M in the image.

We also introduce a world reference frame (denoted R w ). A rigid-body transformation makes the link between the coordinates of M in the different reference frames:

which will alternatively be written:

(with O c and O w the origin of R c and R w respectively), R is a rotation matrix. We also consider that the surface to which M belongs can be parameterised as: