Frédéric Sur

Benoît Blaysat

Michel Grédiac

On biases in displacement estimation for image registration, with a focus on photomechanics

Keywords: Image registration, digital image correlation, fattening effect, Savitzky-Golay filtering, experimental mechanics, stereo-imaging, optical flow

Image registration under small displacements is the keystone of several image analysis tasks such as optical flow estimation, stereoscopic imaging, or full-field displacement estimation in photomechanics. A popular approach consists in locally modeling the displacement field between two images by a parametric transformation and performing least-squares estimation afterwards. This procedure is known as "digital image correlation" (DIC) in several domains as in photomechanics. The present article is part of this approach. First, the estimated displacement is shown to be impaired by biases related to the interpolation scheme needed to reach subpixel accuracy, the image gradient distribution, as well as the difference between the hypothesized parametric transformation and the true displacement. A quantitative estimation of the difference between the estimated value and the actual one is of importance in application domains such as stereoscopy or photomechanics, which have metrological concerns. Second, we question the extent to which these biases could be eliminated or reduced. We also present numerical assessments of our predictive formula in the context of photomechanics. Software codes are freely available to reproduce our results. Although this paper is focused on a particular application field, namely photomechanics, it is relevant to various scientific areas concerned by image registration.

Introduction

Two images I and I of a scene being given, image registration consists in estimating a bidimensional displacement field φ which permits to map pixels from one image to the corresponding pixels in the other image. Assuming that corresponding pixels have the same intensity, for any pixel x, the following relation thus holds: I(x) = I (x + φ(x)). When small displacements are sought (typically below one pixel), the displacement field is often estimated locally by imposing the preceding relation over corresponding subsets of the domains of I and I . While image registration is the very first step of many image analysis tasks, the case of small amplitude displacements is of particular interest in optical flow estimation [START_REF] Fleet | Optical flow estimation[END_REF][START_REF] Fortun | Optical flow modeling and computation: A survey[END_REF], disparity estimation in stereoscopic imaging [START_REF] Delon | Small baseline stereovision[END_REF][START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF][START_REF] Scharstein | High-resolution stereo datasets with subpixel-accurate ground truth[END_REF][START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF], or extensometry in experimental mechanics [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties: a review[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF], to mention some representative examples. Although we will focus our presentation on displacement estimation in photomechanics, which is a branch of experimental solid mechanics dedicated to full-field measurements from images, many aspects of the problem discussed in this paper are of general interest. However, in this lattest application domain, no occlusions are present, in contrast to stereoscopic imaging or optical flow estimation.

Problem statement

In photomechanics, digital image correlation (DIC) [START_REF] Hild | Digital image correlation: from displacement measurement to identification of elastic properties: a review[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF] is of prime interest to measure displacement fields on the surface or in the bulk of materials subjected to thermomechanical loads. The term DIC equally refers to methods which are based on crosscorrelation (CC), on sum of squared differences (SSD), or on normalized CC or SSD as well [START_REF] Pan | Equivalence of digital image correlation criteria for pattern matching[END_REF]. DIC methods are based on two images I and I of the surface of the specimen taken before and after deformation, respectively. The specimen shall be marked beforehand with a contrasted random pattern, called a speckle pattern, as shown in Section 3. The aim is to retrieve the displacement field u such that for any x, I(x) = I (x + u(x)). Strain fields are then deduced by differentiation. These fields are used to observe various phenomena, which occur on the surface of deformed specimens, and which are revealed by displacement or strain heterogeneities. They can eventually be used to identify parameters governing constitutive equations, which are then used to design structural components [START_REF]Full-field measurements and identification in solid mechanics[END_REF]. It is therefore of prime importance to retrieve displacement and strain maps affected by the lowest possible measurement errors, so that the resulting identified parameters are themselves as reliable as possible. The displacement field usually features values well below one pixel. The so-called local subset-based DIC consists in registering I on I by optimizing some criterion as the above-mentioned CC or SSD over subsets of the image domains. These criteria being equivalent under mild assumptions [START_REF] Pan | Equivalence of digital image correlation criteria for pattern matching[END_REF], we focus on the sum of squared differences (SSD), which is also used in block matching for stereoscopic imaging [START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF]. Since I and I are known only at integer pixel coordinates and since subpixel accuracy of u is sought, interpolation is required. In order to estimate the displacement u at a pixel x, the following SSD, defined over a subset Ω x of the image of the specimen surface centered at x, is minimized with respect to a displacement φ x defined over Ω x :

x i ∈Ωx I(x i ) -I (x i + φ x (x i )) 2 (1)
Here I denotes a continuous interpolation of I (I being sampled at integer coordinates), and Ω x is a set of M pixels (x 1 , . . . , x M ). Minimization is performed with respect to a displacement field φ x expected to approximate the actual unknown displacement u over Ω x . We do not consider a weighted SSD criterion, but the formulas in the remainder of the paper would easily adapt.

An estimation φ(x) of the displacement field over the whole specimen surface is eventually obtained by taking, at any pixel x, the value of φ x at the center of the subset Ω x .

Displacement estimation is an ill-posed problem because of under-determination. Indeed, at each pixel, a bidimensional displacement must be retrieved. Moreover, information is lost from the component of the displacement orthogonal to the image gradient. This is the so-called aperture problem [START_REF] Fleet | Optical flow estimation[END_REF][START_REF] Fortun | Optical flow modeling and computation: A survey[END_REF]. Consequently, in experimental mechanics (and in general applications as well [START_REF] Bouthemy | Robust model selection in 2D parametric motion estimation[END_REF]), the displacement φ x is usually sought as the linear combination of N shape (or basis) functions (φ j ) 1≤j≤N , such that the parameters (λ j ) 1≤j≤N minimize SSD(λ 1 , . . . , λ N ) =

x i ∈Ωx I(x i ) -I x i + N j=1 λ j φ j (x i ) 2 (2)
It may be noted that classic ways to deal with the aperture problem are either to consider a constant displacement φ x over Ω x and a first-order Taylor expansion of the preceding equation, giving fast approaches à la Lucas-Kanade [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], or to consider global smoothing with regularization constraints, giving approaches à la Horn-Schunck [START_REF] Horn | Determining optical flow[END_REF]. In experimental mechanics, it is often preferred to numerically minimize SSD because an accurate estimation of φ is required and setting the hyperparameter involved by regularization is not an easy task, which rules out such approaches.

From now on, we assume that Λ = (λ 1 , . . . , λ N ) T minimizes the expression of SSD given by Equation 2, and φ x = N j=1 λ j φ j denotes the corresponding local displacement field. For example, as recalled in [START_REF] Bornert | Assessment of digital image correlation measurement errors: Methodology and results[END_REF], zero-order shape functions are such that N = 2 and φ 1 (x) = (1 0) T , φ 2 (x) = (0 1) T , giving a constant φ x over Ω x . First-order shape functions are such that N = 6 and functions φ 3 (x) = (x 0) T , φ 4 (x) = (0 x) T , φ 5 (x) = (y 0) T , φ 6 (x) = (0 y) T are added to φ 1 and φ 2 . Second-order shape functions are such that N = 12 and embed the following additional functions: φ 7 (x) = (x 2 , 0) T , φ 8 (x) = (xy, 0) T , φ 9 (x) = (y 2 , 0) T , φ 10 (x) = (0, x 2 ) T , φ 11 (x) = (0, xy) T , φ 12 (x) = (0, y 2 ) T .

The origin of the axis is often the center of the subset Ω x , so that the displacement estimated at x by minimizing the SSD criterion over Ω x is given by φ(x) = φ x (0, 0) = (λ 1 , λ 2 ) T for any order. Moreover, one can see that first-and second-order shape functions satisfy ∂φ x /∂x(0, 0) = (λ 3 , λ 4 ) T and ∂φ x /∂y(0, 0) = (λ 5 , λ 6 ) T .

The displacement field is supposed to be smooth enough so that approximating u by φ makes sense. In particular, no occlusions (an object hide the part of another object in stereoscopy) or missing parts (as a mechanical fracture) are allowed. Displacement is locally invertible if the mapping from x to x + φ(x) has a non-singular Jacobian matrix, thanks to the inverse function theorem. In the case of polynomial shape functions, the determinant of this matrix at the center of the subset is (1 + λ 3 )(1 + λ 6 ) -λ 4 λ 5 = 1 + λ 3 + λ 6 + λ 3 λ 6 -λ 4 λ 5 . Of course, zero-order shape functions (corresponding to simple translations) give invertible displacement since, in this case, λ i = 0 for i ∈ {3, 4, 5, 6} and the determinant is nonzero. Higher-order shape functions also yield locally invertible displacement fields if the λ i 's are small with respect to 1, which holds true under the classic small strain hypothesis in photomechanics.

Motivation

The motivation of our work is the observation that, in photomechanics, displacement fields are often impaired by spurious small fluctuations. These fluctuations are often believed, in this community, to be caused by sensor noise, interpolation bias, or the numerical scheme minimizing the DIC criterion. Nevertheless, recent works consider the random marking on the surface of the specimen [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF][START_REF] Lehoucq | The effect of the ill-posed problem on quantitative error assessment in digital image correlation[END_REF][START_REF] Wang | Theoretical analysis on the measurement errors of local 2D DIC: Part I temporal and spatial uncertainty quantification of displacement measurements[END_REF], in addition to the aforementioned causes, the authors of [START_REF] Fayad | Spatial DIC errors due to pattern-induced bias and grey level discretization[END_REF] coining the term pattern-induced bias. From a different perspective, it is known that the local parametric estimation of a displacement field is intrinsically biased, the retrieved displacement being merely the convolution of the true displacement by a Savitzky-Golay filter [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]. In this latter approach, specimen marking does not play any role, which is contradictory with the presence of a pattern-induced bias. The main motivation of this study is thus to see if this contradiction can be removed.

Contributions

The goal of the present paper is to provide a theoretical basis for these different viewpoints and thus to show that they can coexist. More precisely, we shall give the relation between the retrieved Λ (thus, the retrieved displacement φ) and the actual unknown displacement u over any subset Ω x . By revisiting papers dealing with disparity estimation in stereo-imaging [START_REF] Blanchet | Fattening free block matching[END_REF][START_REF] Delon | Small baseline stereovision[END_REF][START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF][START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF], we show in Section 2 that the estimation of Λ is mainly the sum of two terms, namely a term depending on the unknown displacement u and the approximation induced by the shape functions, and a term caused by the interpolation error. We propose a characterization of these terms, which actually depend on the gradient distribution over the imaged specimen. This unifies the points of view of stereoscopy (fattening effect) and photomechanics (pattern-induced bias). In Section 3, numerical experiments validate the theoretical developments and assess their limitations. To limit the size of the paper, we discuss how sensor noise propagates from the images to the estimated displacement field in a separate companion research report [START_REF] Sur | On biases in displacement estimation for image registration, with a focus on photomechanics -Extended version[END_REF].

Related work

Stereoscopy and photomechanics share the same objective, namely estimating displacement fields very accurately. The proposed contribution echoes different papers from these research areas.

First, the proposed approach relies on the calculation presented in [START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF] (and, to some extent, in [START_REF] Delon | Small baseline stereovision[END_REF]) in the context of disparity estimation in stereo-imaging. In this research field, disparity plays the role of the displacement field considered in the present paper. On the one hand, disparity may be estimated at any pixel as a constant monodimensional displacement between small patches extracted from the stereo image pair. Indeed, disparity is collinear with a given direction, since images are rectified so that the epipolar lines are parallel to each other. In this context, the authors of [START_REF] Delon | Small baseline stereovision[END_REF][START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF] give predictive formulas for quantifying the fattening effect, i.e., the bias in disparity estimation caused by image gradient distribution. On the other hand, in experimental mechanics, local non-constant displacements are estimated as linear combinations of shape functions and they are not constrained to be collinear with a given direction. The main difference with disparity is that displacement fields are usually smooth and have tiny fluctuations, the strain components (defined from the partial derivatives of the displacement) being below 10 -2 for many materials and load intensity. Gradient distribution is also likely to affect displacement estimation, giving pattern-induced bias. A first attempt at characterizing this bias in photomechanics is available in [START_REF] Lehoucq | The effect of the ill-posed problem on quantitative error assessment in digital image correlation[END_REF].

In addition to pattern-induced bias or fattening effect, the interpolation scheme required for registration is certainly a source of error. It is a common assumption in stereo-imaging that the input images satisfy Shannon-Nyquist sampling conditions, cf. [START_REF] Delon | Small baseline stereovision[END_REF][START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF][START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. This involves that continuous images can be perfectly interpolated from the Fourier coefficients without any interpolation bias. However, it is mentioned in [START_REF] Reu | All about speckles: Aliasing[END_REF][START_REF] Sur | Determining displacement and strain maps immune from aliasing effect with the grid method[END_REF] that aliasing, although hardly noticeable to the naked eye on the raw input images, may strongly affect the estimated displacement and strain fields in photomechanics. This motivates the use of bilinear or bicubic interpolation schemes in this field [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF] instead of Fourier interpolation. The drawback is that I matches I only at integer pixel coordinates: non-integer pixel coordinates are thus affected by interpolation error, as numerically illustrated in [START_REF] Bornert | Shortcut in DIC error assessment induced by image inerpolation used for subpixel shifting[END_REF]. Seeing interpolation as a convolution filter, the authors of [START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF] have proposed a characterization of the interpolation-induced bias giving rise to the famous "S-shape function". The authors of [START_REF] Su | Elimination of systematic error in digital image correlation caused by intensity interpolation by introducing position randomness to subset points[END_REF] reduce the interpolation bias by sampling the subset Ω at non-integer pixels.

As mentioned in the introduction, the most popular approach in experimental mechanics is certainly to locally approximate the displacement field by a linear combination of shape functions, which may in turn undermatch the true displacement. This is another source of error discussed in [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF], where it is shown that the retrieved displacement field is the convolution of the true displacement with a Savitzky-Golay filter characterized by the order of the shape functions and the size of the subset. This characterization is used in [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF] and [START_REF] Xu | Theoretical estimation of systematic errors in local deformation measurements using digital image correlation[END_REF] to define some metrological parameters, and in [START_REF] Grédiac | A robust-to-noise deconvolution algorithm to enhance displacement and strain maps obtained with local DIC and LSA[END_REF] to restore displacement fields through a dedicated deconvolution procedure.

Notation and reminder

In what follows, •, • and | • | denote the bidimensional Euclidean product and norm, respectively. The gradient of any 2D function is denoted by ∇. We identify any vector with the corresponding column matrix, in boldface letters. The transpose of any matrix A is denoted by A T .

We denote by O Landau's "big-O" for a variable tending to 0. We recall that if f and g are 2D-valued functions, f = O(|g|) means that for some K > 0, any small enough x

satisfies |f (x)| ≤ K|g(x)|. If f = O(|g|) and A is some 2-column matrix, then Af = O(|g|)
and in particular h, f = O(|g|) for any constant 2-D vector h.

In the remainder of this paper, we will make use of partial derivatives of interpolated 2-D functions. While some interpolation schemes, such as bilinear interpolation, do not provide us with derivatives at integer pixel coordinates, it should be noted that derivatives of the interpolated functions will be calculated at points like x i + φ(x i ), x i being an integer pixel coordinate. In most situations, φ(x i ) has a subpixel value, thus x i + φ(x i ) has a noninteger value and the derivatives are well defined.

Estimating displacements by minimizing SSD

This section gives a closed-form expression of the displacement φ estimated by minimizing the SSD criterion (defined in Equation 2) as a function of the actual unknown displacement u. The aim is to emphasize the role of systematic errors caused by image texture and interpolation scheme. Non-noisy images are considered in this section. Images affected by signal-dependent noise are investigated in a separate research report [START_REF] Sur | On biases in displacement estimation for image registration, with a focus on photomechanics -Extended version[END_REF].

Since the present section deals with the minimization problem at a given pixel x and in order to simplify notations, we do not write the index x to remind that the subset is centered at a pixel x, and we simply write Ω and φ instead of Ω x and φ x .

Relation between retrieved and actual displacement fields

The goal here is to express Λ as a function of the images I and I , of the interpolated continuous image I , and of the unknown displacement field u.

By definition, Λ = (λ 1 , . . . , λ N ) T is a stationary point of the SSD given by Equation 2. For any j ∈ {1, . . . , N }, taking the derivative with respect to λ j thus gives:

x i ∈Ω I(x i ) -I (x i + φ(x i )) • ∇ I (x i + φ(x i )) , φ j (x i ) = 0 (3)
Let δ be the difference between the unknown displacement field and the retrieved one, such that for any x i ∈ Ω, δ(x i ) = u(x i ) -φ(x i ). It quantifies the undermatching of shape functions mentioned in [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]. First,

I(x i ) -I (x i + φ(x i )) = I(x i ) -I (x i + φ(x i )) + I (x i + φ(x i )) -I (x i + φ(x i )) (4)
By definition of u,

I(x i ) -I (x i + φ(x i )) = I (x i + u(x i )) -I (x i + φ(x i )) (5)
A first-order Taylor series expansion allows writing:

I(x i ) -I (x i + φ(x i )) = ∇I (x i + u(x i )), δ(x i ) + O(|δ(x i )| 2 ) (6)
Denoting by DI = I -I the interpolation error (equal to 0 at integer coordinates), we conclude from Equation 4that

I(x i ) -I (x i + φ(x i )) = ∇I (x i + u(x i )) , δ(x i ) + DI (x i + φ(x i )) + O(|δ(x i )| 2 ) (7)
Second, another Taylor series expansion gives:

∇I (x i + φ(x i )) -∇I (x i + u(x i )) = -δ(x i ) T HI (x i + φ(x i )) δ(x i ) + O(|δ(x i )| 3 ) = O(|δ(x i )| 2 ) (8)
where HI (x i + φ(x i )) is the Hessian matrix of I at x i + φ(x i ). Consequently,

∇ I (x i + φ(x i )) = ∇I (x i + u(x i )) + ∇ I (x i + φ(x i )) -∇I (x i + φ(x i )) + ∇I (x i + φ(x i )) -∇I (x i + u(x i )) = ∇I (x i + u(x i )) -D∇I (x i + φ(x i )) + O(|δ(x i )| 2 ) (9)
where D∇I = ∇I -∇ I is the gradient interpolation error. Plugging Equations 7 and 9 into Equation 3 gives:

x i ∈Ω ∇I (x i + u(x i )), δ(x i ) + DI (x i + φ(x i )) + O(|δ(x i )| 2 ) × ∇I (x i + u(x i )) -D∇I (x i + φ(x i )) + O(|δ(x i )| 2 ), φ j (x i ) = 0 (10)
For any i ∈ {1, . . . , M } and j ∈ {1, . . . , N }, let L u i,j = ∇I (x i + u(x i )) , φ j (x i ) , such that L u is a M × N matrix. We eventually obtain:

x i ∈Ω L u i,j ∇I (x i + u(x i )), δ(x i ) + L u i,j DI (x i + φ(x i )) + O(|δ(x i )| 2 ) + O |D∇I (x i + φ(x i )) | • |δ(x i )| + O |DI (x i + φ(x i )) | • |D∇I (x i + φ(x i )|) + O |DI (x i + φ(x i )) | • |δ(x i )| = 0 (11)
or in a simpler way:

x i ∈Ω L u i,j ∇I (x i + u(x i )), δ(x i ) + x i ∈Ω L u i,j DI (x i + φ(x i )) + O(|δ| 2 ) + O |DI | • |D∇I | + O |DI | • |δ| + O |D∇I | • |δ| = 0 (12)
where

O(|f |) denotes x i ∈Ω (O|f (x i )|). Since δ(x i ) = u(x i ) -φ(x i ) = u(x i ) -N k=1 λ k φ k (x i
), the preceding equation gives:

x i ∈Ω L u i,j ∇I (x i + u(x i )), u(x i ) - x i ∈Ω L u i,j N k=1 λ k L u i,k + x i ∈Ω L u i,j DI (x i + φ(x i )) + O(|δ| 2 ) + O |DI | • |D∇I | + O |DI | • |δ| + O |D∇I | • |δ| = 0 (13)
If we denote by G the vector of components G i = ∇I (x i + u(x i )), u(x i ) for any i ∈ {1, . . . , M }, we obtain with Equation 13 the following matrix relation:

(L u ) T L u Λ = (L u ) T G + (L u ) T DI + O(|δ| 2 ) + O |DI | • |D∇I | + O |DI | • |δ| (14)
Matrix (L u ) T L u is invertible as soon as the columns of L u are linearly independent (it is a Gramian matrix), which holds if the φ j form a valid basis and if the gradient is not equal to zero over the whole subset Ω. We assume in the following that these assumptions hold.

We have finally demonstrated the following theorem.

Theorem 1 If Λ minimizes the SSD criterion of Equation 2 registering image I over I, then the following equality holds:

Λ = ((L u ) T L u ) -1 (L u ) T G + ((L u ) T L u ) -1 (L u ) T DI + O(|δ| 2 ) + O |DI | • |D∇I | + O |DI | • |δ| ( 15
)
where u is the actual displacement field, and for any i ∈ {1, . . . , M } and j ∈ {1, . . . , N }, the following notations hold: With this presentation, it is not straightforward to see how the estimation of Λ is affected by a systematic error.

L u i,j = ∇I (x i + u(x i )) , φ j (x i ) , G i = ∇I (x i + u(x i )), u(x i
Let us now decompose (for instance, in the least-squares sense) the unknown, true displacement u over the basis functions (φ k ). We introduce the N -dimensional vector Λ u = (λ u 1 , . . . , λ u N ) T such that u = N k=1 λ u k φ k + where is the residual in the decomposition of u over the basis (φ k ) k∈{1,...,N } . Consequently, for any i ∈ {1, . . . , M },

G i = N k=1 λ u k ∇I (x i + u(x i )), φ k (x i ) + ∇I (x i + u(x i )), (x i ) (16)
If we denote by E the vector of components E i = ∇I (x i + u(x i )), (x i ) , we obtain G = L u Λ + E. Consequently, Theorem 1 gives the following corollary.

Corollary 1 If Λ minimizes the SSD criterion of Equation 2 registering image I over I, then the following equality holds: [START_REF] Grédiac | Comparing several spectral methods used to extract displacement fields from checkerboard images[END_REF] with, in addition to the notations of Theorem 1, E i = ∇I (x i + u(x i )), (x i ) for any i ∈ {1, . . . , M }, where is the component of u out of the basis (φ k ).

Λ = Λ u + ((L u ) T L u ) -1 (L u ) T E + ((L u ) T L u ) -1 (L u ) T DI + O(|δ| 2 ) + O |DI | • |D∇I | + O |DI | • |δ|

Discussion

Skipping second-order terms from Theorem 1 gives:

Λ = ((L u ) T L u ) -1 (L u ) T G + ((L u ) T L u ) -1 (L u ) T DI (18)
and from Corollary 1:

Λ = Λ u + ((L u ) T L u ) -1 (L u ) T E + ((L u ) T L u ) -1 (L u ) T DI (19)
It is easy to see that if interpolation is perfect (DI = 0), and if the function basis is expressive enough to perfectly represent the true displacement u (E = 0) or at least if the degree of the polynomial basis is large enough (E 0) , then Λ is equal to the sought Λ u .

These assumptions are, however, optimistic. They are discussed in the remainder of this section.

Interpolation bias

In Equation 18 the term ((L u ) T L u ) -1 (L u ) T DI quantifies the effect of subpixel interpolation error. Since the underlying image I is unknown, it is not possible to bound a priori this error, except if some additional information is available. For instance, well-sampled images satisfying the Shannon-Nyquist condition give no interpolation error when using Fourier interpolation as in [START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF][START_REF] Szeliski | Sampling the disparity space image[END_REF] and in this case, DI = 0. As discussed in the introduction, real images from photomechanical experiments probably do not satisfy these hypotheses.

It should be noted that most image interpolation schemes (such as Fourier, bilinear, bicubic, or Lanczos interpolations) are actually linear, in the sense that any interpolated value is a weighted mean of image values at integer pixel coordinates [START_REF] Getreuer | Linear methods for image interpolation[END_REF][START_REF] Lehmann | Survey: interpolation methods in medical image processing[END_REF]. For any x ∈ R 2 , there exists a row matrix P (x) of size M such that the interpolation of I at x is given by I (x) = P (x)I , where I momentarily denotes the matrix of image values reshaped as a column vector. For instance, row-vector P (x) has four non-zero entries in bilinear interpolation, sixteen in bicubic or bicubic spline interpolation [START_REF] Getreuer | Linear methods for image interpolation[END_REF][START_REF] Lehmann | Survey: interpolation methods in medical image processing[END_REF], which are the most popular schemes in photomechanics.

The interpolation bias thus writes:

((L u ) T L u ) -1 (L u ) T (I -P (x)) I (20)
It has been derived by another approach in [START_REF] Blaysat | Effect of interpolation on noise propagation from images to DIC displacement maps[END_REF] (refining [START_REF] Réthoré | Experimental investigation of localized phenomena using digital image correlation[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements[END_REF]), and experimentally assessed in [START_REF] Blaysat | On the propagation of camera sensor noise to displacement maps obtained by DIC -An experimental study[END_REF]. It can be noted that the image gradient is involved in L u and affects the interpolation bias. While interpolation error may be neglected (either by assuming the Shannon-Nyquist condition to be satisfied or by using high-order interpolation schemes), the next section deals with the so-called "undermatched subset shape functions" which potentially gives pattern-induced bias.

Undermatched shape functions and pattern-induced bias

In Equation 18, the term ((L u ) T L u ) -1 (L u ) T G links Λ with the gradient of I and the actual displacement u. Let us recall that, for any i ∈ {1, . . . , M } and j ∈ {1, . . . , N },

L u i,j = ∇I (x i + u(x i )) , φ j (x i ) (21)
G i = ∇I (x i + u(x i )), u(x i ) . (22)
and

E i = ∇I (x i + u(x i )), (x i ) . (23)
We can see that we find again the well-known aperture problem: only the displacement component collinear with the image gradient plays a role in G i . Moreover, Theorem 1 (Corollary 1, respectively) shows that each component (the error on each component, respectively) of Λ is a weighted mean of G i (of the undermatching error E, respectively), the weights being proportional to the squared components of the gradient along the shape functions (here, ((L u ) T L u ) -1 acts as a normalization).

For didactic purposes, we explicit the relation in a simple, yet realistic case, where displacements is sought as a pure translation along the x-direction. This situation corresponds to disparity estimation in stereo-imaging with rectified images. We have u(x i ) = (u(x i ) 0) T , N = 1 and φ 1 (x) = (1 0) T . In this case, L u = (I x (x i + u(x i ))) 1≤i≤M is a column-vector, and:

((L u ) T L u ) -1 (L u ) T G = M i=1 (I x (x i + u(x i ))) 2 u(x i ) M i=1 (I x (x i + u(x i ))) 2 (24)
where I x denotes the partial derivative of I along direction x.

In this case, the term related to the interpolation error simplifies into:

((L u ) T L u ) -1 (L u ) T DI = M i=1 I x (x i + u(x i )) DI (x i ) M i=1 (I x (x i + u(x i ))) 2 (25)
When u = (u, 0) is a constant displacement, Equation 24 simplifies to u. This observation is consistent with Corollary 1 since in this case, the displacement is perfectly represented by the basis function φ 1 as, for any x, u(x) = u φ 1 (x).

It should be noted that Equation 24is exactly the relation given in [START_REF] Blanchet | Fattening free block matching[END_REF][START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF]. When the sought displacement is not constant over Ω, it is a weighted mean of the actual displacement, the weights being the squared image derivatives. This justifies the well-known fattening effect in stereo-imaging (also called adhesion effect in [START_REF] Delon | Small baseline stereovision[END_REF]): points lying in the neighborhood of edges have a disparity essentially governed by the edge points, which have large gradient values. Fattening effects can be seen for example in [START_REF] Sabater | How accurate can block matches be in stereo vision[END_REF][START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]: foreground objects (which have a larger disparity than the background) appear fatter than they are because background pixels near their edges inherit their disparity.

In experimental mechanics the relevant formula is given by Theorem 1. The estimated displacement field given by Λ is biased since it is a weighted sum of the actual displacement at the pixels belonging to the subset, the weights increasing with the squared gradient at the pixels (more precisely, with the squared component of the gradient collinear with the shape functions). In particular, the displacement returned at the center of the considered subset is essentially the average of the actual displacement taken at high-gradient points, even if the actual displacement is different at the considered point. Although this corresponds to the intuition, we are not aware of earlier papers giving a rigorous description of the phenomenon in the context of experimental mechanics. As a consequence, we keep on calling this phenomenon "pattern-induced bias" after [START_REF] Fayad | Spatial DIC errors due to pattern-induced bias and grey level discretization[END_REF], since the expression "fattening effect" does not seem to be adequate in experimental mechanics, nothing really becoming "fatter".

Note that the relation between Λ and G in Theorem 1 can be seen as a generalized convolution, with a spatially varying kernel.

From pattern-induced bias to Savitzky-Golay filtering

The authors of [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF] claim that the estimated displacement is simply the convolution product between the true displacement and the Savitzky-Golay (SG) kernel (a low-pass filter [START_REF] Savitzky | Smoothing and differentiation of data by simplified leastsquares procedures[END_REF][START_REF] Luo | Properties of Savitzky-Golay digital differentiators[END_REF]), causing the so-called matching bias [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF]. The SG filter only depends on the degree of the polynomial shape functions and on the size of the considered subset Ω [START_REF] Savitzky | Smoothing and differentiation of data by simplified leastsquares procedures[END_REF]. In particular, it does not depend on the gradient of the underlying image. Moreover, the discussion of Section 2.2.2 concludes that the relation between the estimated displacement and the true one can be seen as a convolution with a spatially varying kernel. The claim of [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF], backed by results from [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF] or [START_REF] Grédiac | A robust-to-noise deconvolution algorithm to enhance displacement and strain maps obtained with local DIC and LSA[END_REF], therefore seems to contradict Theorem 1 and pattern-induced bias discussed in the preceding section. In the remainder of this section, we explain how these two viewpoints can be accomodated.

The case of stationary random patterns

To establish the relation with the SG filter, the ground hypothesis of [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF] is that, with our notations, φ minimizes

x i ∈Ωx |u(x i ) -φ(x i )| 2 (26)
over a subset Ω x centered at pixel x, for any x. Let us write for a while the x-and y-components of φ as:

φ x (x) = 0≤i+j≤d α i,j x i y j (27)
φ y (x) = 0≤i+j≤d β i,j x i y j ( 28
)
where d is the maximum degree of the shape functions. Thus, λ 1 = α 0,0 , λ 2 = β 0,0 ,

λ 3 = α 1,0 , λ 4 = β 1,0 , λ 5 = α 0,1 , λ 6 = β 0,1 , etc.
With J x = (1, x, . . . , x d , y, . . . , yx d-1 , y 2 , . . . , y 2 x d-2 , . . . , y d ) T , α = (α 0,0 , α 1,0 , . . . , α d,0 , α 1,0 , . . . , α 1,d-1 , α 2,0 , . . . , α 2,d-2 , . . . , α d,0 ) T , and β = (β 0,0 , β 1,0 , . . . , β d,0 , β 1,0 , . . . , β 1,d-1 , β 2,0 , . . . , β 2,d-2 , . . . , β d,0 ) T , Equation 26 writes:

x i u 1 (x i ) -J T x i α 2 + u 2 (x i ) -J T x i β 2 (29)
Consequently, minimizing Equation 26 with respect to α and β amounts to solving the normal equations J T α = u 1 and J T β = u 2 where matrix J collects all column vectors J x . In other words, α = (JJ T ) -1 Ju 1 and β = (JJ T ) -1 Ju 2 . Since the displacement estimated over the subset Ω is simply (λ 1 , λ 2 ) = (α 0,0 , β 0,0 ), each component of the retrieved displacement is the convolution of the components of u by a Savitzky-Golay [START_REF] Savitzky | Smoothing and differentiation of data by simplified leastsquares procedures[END_REF] filter of order equal to the degree of the shape functions and of support given by the dimensions of Ω. With first-and second-degree shape functions, the partial derivatives of the displacement field being given by (λ 3 , λ 4 , λ 5 , λ 6 ), they are also given by SG filters, as explained in [START_REF] Luo | Properties of Savitzky-Golay digital differentiators[END_REF][START_REF] Savitzky | Smoothing and differentiation of data by simplified leastsquares procedures[END_REF]. The justification that the sought φ minimizes Equation 26 is only based on heuristic arguments in [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]. We shall see that the very specific nature of the speckle patterns used in experimental mechanics permits to justify this claim.

Since I(x i ) = I (x i + u(x i )), if we neglect interpolation error and identify I with I , we obtain the following first-order approximation as in Equation 6:

I(x i ) -I (x i + φ(x i )) = ∇I (x i + u(x i )), u(x i ) -φ(x i ) (30)
Provided this first-order approximation holds (for instance because an initial guess of the solution is available), minimizing the SSD criterion thus amounts to minimizing:

x i ∈Ω g 2 i cos 2 (θ i ) |u(x i ) -φ(x i )| 2 (31)
with g i the norm of ∇I (x i + u(x i )) and θ i the angle between ∇I (x i + u(x i )) and the error u(x i ) -φ(x i ). As can be seen, this is not Equation 26. It can be noted that the component of uφ orthogonal to the gradient of I does not play any role, which is consistent with the preceding discussion about the aperture problem. This was also mentioned in [START_REF] Lehoucq | The effect of the ill-posed problem on quantitative error assessment in digital image correlation[END_REF].

As explained in the introduction, specimens tested in experimental mechanics are marked with random speckle patterns which can be modeled as stationary textures. We can thus safely assume that the gradient norms (g i ) x i ∈Ω are indentically distributed random variables, as well as the angles (θ i ) x i ∈Ω , and that at each pixel x i , g i and θ i are independent. Nevertheless, these random variables are spatially correlated. A common assumption is that spatial correlation vanishes with the distance between pixels as in natural images [START_REF] Su | Theoretical analysis on performance of digital speckle pattern: uniqueness, accuracy, precision, and spatial resolution[END_REF].

Let

X i = g 2 i cos 2 (θ i ) |u(x i ) -φ(x i )| 2 for any 1 ≤ i ≤ M .
The classic law of large numbers does not hold here because of spatial correlations. However, generalizations such as Bernstein's weak law of large numbers [8, Ex. 254 p. 67] still hold. Assuming that the variance of the X i is bounded, i.e., there exists c > 0 such that Var(X i ) ≤ c, and assuming also that spatial correlations vanish with distance, i.e., Cov(X i , X j ) → 0 when

|x i -x j | → +∞, we obtain 1 : Var M i=1 X i = M i=1 Var(X i ) + 2 M i=1 j∈[1,M ]\V i Cov(X i , X j ) + 2 M i=1 j∈[1,M ]∩V i Cov(X i , X j ) (32)
1 The proposed calculation is adapted from https://math.stackexchange.com/questions/245327/ where, for any ε > 0, V i is the set of indices j such that for any i and j ∈ V i , Cov(X i , X j ) ≤ ε. Here, B\A denotes the relative complement of a set A with respect to B. Let N i be the cardinality of N\V i which is a finite set since Cov(X i , X j ) → 0 when |x i -x j | → +∞, and

N = max i N i .
Since for any i, j, Cov(X i , X j ) ≤ c by Cauchy-Schwartz inequality, the following upper bound holds:

Var M i=1 X i ≤ M c + 2M N c + 2M 2 ε (33) Thus, Var 1 M M i=1 X i ≤ c/M + 2N c/M + 2ε (34)
Chebyshev's inequality implies that the random variable

1 M M i=1 X i -1 M M i=1 E(X i ) tends to 0 in probability as M → +∞.
In other words, this justifies that minimizing Equation 31 amounts to minimizing

M i=1 E(X i ) = E(g 2 )E(cos 2 (θ)) x i ∈Ω |u(x i ) -φ(x i )| 2 (35)
as soon as the size M of the domain Ω is "large enough". Equation 34 states that, for a given M , the approximation is as tight as the bound c on the variance is small, or as the spatial correlation of the image gradients quickly vanishes, giving a small N .

Towards an optimum pattern with respect to pattern-induced bias?

The result of the preceding section can be interpreted as follows: minimizing Equation 31to estimate φ is similar to minimizing Equation 26, provided that the size M of the subset Ω is large enough. From Equation 34, this is all the more valid as the speckle pattern is fine (giving quickly vanishing spatial correlations, thus a small N for a given ε) and as the variance of the X i is small (giving a small c). Since Var(X i ) is proportional to Var(

g 2 i cos 2 (θ i )) = E(g 4 i )E(cos 4 (θ i )) -E 2 (g 2 i )E 2 (cos 2 (θ i )) = 3E(g 4 i )/8 -E 2 (g 2 i )
/4 (assuming that the θ i are uniformly distributed in the interval (0, 2π), which is sound with an isotropic speckle pattern, we indeed obtain E(cos 2 (θ)) = 1/2 and E(cos

4 (θ)) = 3/8). As a consequence, Var(g 2 i cos 2 (θ i )) = 3Var(g 2 i )/8 + E 2 (g 2 i )/8 (36)
A fine pattern minimizing this quantity should have lower pattern-induced bias. One can see that a concentrated gradient distribution with a low gradient average value is of interest. However, such a speckle pattern is still to be designed.

Link with a fattening-free criterion in stereoscopic imaging

Disparity in stereo-imaging being a 1-D displacement, θ i = 0 in Equation 31 and g i = I x (x i + u(x i )) where I x denotes the partial derivative of I along the epipolar line. This motivates the authors of [START_REF] Blanchet | Fattening free block matching[END_REF] to estimate φ by minimizing the following weighted SSD criterion, denoted SSD in what follows:

x i ∈Ωx I(x i ) -I (x i + φ(x)) 2 max(| I x (x i + u(x i ))| 2 , κ) (37)
where κ > 0 avoids divisions by zero.

In stereo-imaging, φ is sought as a constant displacement. The solution of SSD minimizes i |u(x i ) -φ| 2 and does not depend on the image gradient. It is shown in [START_REF] Blanchet | Fattening free block matching[END_REF] that estimating φ by minimizing SSD at any pixel gives an estimated displacement at x which is the mean of the u(x i ) over Ω x . Interestingly, this is consistent with the preceding section, since a zero-order Savitzky-Golay filter is a simple moving average with a kernel constant over its domain.

Numerical assessment

The goal of this section is to provide the reader with illustrative didactic experiments, and to assess the validity of the predictive formulas given by Theorem 1 and Corollary 1. We also discuss to what extent PIB can be eliminated or decreased. In the proposed numerical assessments, the real displacement fields are known, and we are able to compare the estimated displacement to this ground truth. It should be noted that the numerical scheme actually used to minimize the SSD criterion over the subsets (Equation 2) (see, e.g., [START_REF] Passieux | Classic and inverse compositional Gauss-Newton in global DIC[END_REF]) is not important here. Nevertheless, the stopping criterion must be set carefully so that the stationarity assumption, which is the ground of Section 2, is valid. In practice, we use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method with a cubic line search implemented in Matlab's fminunc function.

Since the predictive formulas are based on the ideal, continuous images I and I whose derivatives are required, we shall first consider in Section 3.1 images and deformation fields expressed as simple closed-form expressions. We also assess the effect of gray-level quantization. Nevertheless, the true derivatives of quantized images are, of course, unknown. Section 3.2 deals with synthetic speckle images which mimic real images used in experimental solid mechanics.

The numerical experiments proposed in this section can be reproduced with datasets and Matlab codes available at the following URL: https://members.loria.fr/FSur/software/PIB/ The interested reader can also easily modify the parameters and datasets for further investigations.

Image pairs given by a closed-form expression

In this section, we make use of images defined by closed-form expressions in order to assess the proposed formulas in a controlled experimental setting.

Displacement estimation

We define the image I of the reference state and the image I of the deformed state at any pixel of coordinates (x, y) as sine waves by the following equations:

I(x, y) = 2 b-1 (1 + γ sin(2π(x + u(x)) 1.5 /p 1.5 )) (38)
I (x, y) = 2 b-1 (1 + γ sin(2πx 1.5 /p 1.5 )) (39)
where (x, y) spans a 200 × 15 pixel domain, b = 8 (so that the gray level of both images spans an 8-bit range), γ = 0.9 is the contrast, p = 50 pixels governs the varying period of the sine wave, and u(x) is the ground-truth displacement field, supposed to be restricted along the x-axis. See Figure 1. Note that in this section, images are not quantized over b bits.

We seek for a constant displacement field (φ, 0) on each subset Ω x of size 15 × 15 pixels 2 distributed along the x-axis. The size of the subset is chosen in accordance with the characteristic scale of the modulated sine wave giving images I and I . Such a constant displacement corresponds to zero-order shape functions. We therefore calculate at any abscissa x the quantity φ(x) = (φ(x), 0) minimizing the SSD criterion:

x i ∈Ωx I(x i ) -I (x i + φ(x)) 2 . ( 40
)
The numerical assessment in the present section deals with a 1-D displacement along the x-axis in an image which varies only along the x-axis. The gradient is thus always collinear with the displacement. Consequently, the aperture problem manifests itself only because of vanishing image gradients: there is no loss of information orthogonally to the gradient as in the general 2-D case.

The 1-D displacement considered here is the case of interest of stereoscopy. It is possible to also implement the SSD criterion of [START_REF] Blanchet | Fattening free block matching[END_REF], recalled in Section 2.3.3 above. We set the value In each of these three cases, the first row depicts, superimposed on the same graph: the ground truth u, the displacement φ(x) estimated from SSD, the displacement φ(x) estimated from SSD, the displacement predicted by Theorem 1, and Savitzky-Golay filtering of u. In each of the three cases, the second row depicts, superimposed on the same graph : the differences (bias estimations) between u on the one hand, and displacement retrieved with SSD, predicted displacement, displacement retrieved with SSD, output of the SG filter on the other hand.

of κ to achieve the best trade-off between reducing the fattening effect and numerical stability, see Section 3.2.5.

In the SSD and SSD criteria, I denotes a continuous image. Since we use images given by closed-form expressions, image values at non-integer pixels are available. Real experiments require, however, to interpolate images. In this illustrative experiment, we consider three possibilities: 1) using Equation 39 which allows avoiding any interpolation scheme (in this case we also use the closed-form expression of the gradient in the quasi-Newton scheme), 2) using bilinear interpolation, or 3) using bicubic interpolation.

We also consider two displacement fields: a) a constant u(x) = 0.2 pixel, b) a lowfrequency sine wave u(x) = 0.2 sin(2πx/q) where q = 90 pixels. We also discuss in [START_REF] Sur | On biases in displacement estimation for image registration, with a focus on photomechanics -Extended version[END_REF] the case of a high-frequency sine wave (not shown here). In all cases, the largest displacement 1) closed-form expression (first column), 2) bilinear interpolation (second column), and 3) bicubic interpolation (third column). In each of these three cases, the first row depicts, superimposed on the same graph: the ground truth u, the displacement φ(x) estimated from SSD, the displacement φ(x) estimated from SSD, the displacement predicted by Theorem 1, and Savitzky-Golay filtering of u. In each of the three cases, the second row depicts, superimposed on the same graph: the differences (bias estimations) between u on the one hand, and displacement retrieved with SSD, predicted displacement, displacement retrieved with SSD, output of the SG filter on the other hand. value is 0.2 pixel. Figures 2 and3 show various plots. Each of these figures permits discussing displacement fields on the top and biases (systematic errors) on the bottom, as a function of the interpolation scheme (from left to right: closed-form expression, bilinear interpolation, and bicubic interpolation). Concerning the displacement, we plot the ground truth u (thin green line), the displacement φ(x) retrieved by minimizing the SSD criterion (blue), the displacement φ(x) retrieved from the modified criterion SSD defined in Equation 37 (red), the predicted displacement (cf. Theorem 1, yellow) and the filtering of the ground-truth displacement by the Savitzky-Golay filter of order 0 and frame length 15 (purple). Concerning the biases, we plot the differences between the ground truth displacement on the one hand, and the retrieved displacement with SSD, the predicted displacement, the retrieved displacement with SSD, and the output of the SG filter on the other hand.

Concerning the constant displacement a) in Figure 2, since it can be represented with zero-order shape functions, the error (thus E) in Corollary 1 is null: no marking bias should be noticed. With the closed-form expression of image I (case 1), the interpolation error DI is also null. It can be seen, indeed, that all curves are superimposed in this case. With bilinear interpolation (case 2), Theorem 1 predicts an error term caused by the interpolation error. We can see that the retrieved displacement fits well the prediction, the blue and yellow curves being superimposed. We can also see that correcting the marking bias (which is, here, non-existent) significantly amplifies the interpolation error, giving the erratic red curve. When looking closely at the curves showing the biases, we can notice a small difference between the retrieved bias and the predicted bias, which could probably be explained by higher-order error terms than the first-order terms of the present calculation. With bicubic interpolation (case 3), the interpolation error is very small, all the curves being close to each other. With the constant displacement considered here, we can see that interpolation still causes a very small drift in the estimation, giving the increasing bias. Correcting the marking bias with SSD does not give an erratic red curve but a few spurious estimations can be seen. Interpolation errors seem to be amplified.

Concerning the low-frequency sine displacement b) in Figure 3, it cannot be represented by zero-order shape functions, thus marking bias should affect the retrieved displacement, the error term E being non-null in Theorem 1. We can see that with closed-form expression (case 1), thus no interpolation error, the retrieved displacement fits perfectly the prediction, the blue and yellow curves being superimposed. The marking bias causes departures from the ground truth displacement, whose amplitude is governed by the gradient of the underlying image, giving a rather chaotic bias curve. The marking bias is perfectly removed with the SSD criterion: the retrieved displacement indeed fits the output of the SG filter, as predicted by the theory. With bilinear interpolation (case 2), an additional interpolation error is predicted by Theorem 1, although it is difficult to see a difference between the blue and yellow curves. The retrieved displacement globally fits the prediction. With bicubic interpolation (case 3), the curves are close to case 1. When using the SSD criterion, we can see that cases 1 and 3 fit well the predicted Savitzky-Golay filtering (green and red curves are superimposed), but as for displacement a), bilinear interpolation error (case 2) is amplified, giving an erratic red curve.

Effect of quantization

In the preceding section, image intensity is not quantized. We now perform experiments with quantized images in order to illustrate its impact on the predictive formulas. Bicubic interpolation is used in order to minimize interpolation bias, as illustrated earlier. We perform the same experiment as the one described in Figure 3.

Figure 4 shows the result of displacement estimation by minimizing SSD and SSD. The former estimation should match the predicted one, and the latter should match the output of SG filter. We can see that the proposed predictive formulas are quite accurate as soon as quantization is performed over 10 bits (the yellow and blue curves are superimposed), and that retrieving the output of the SG filter with SSD requires to quantize image intensity over 12 bits (so that red and purple curves are superimposed).

Speckle patterns: assessing the prediction for pattern-induced bias

The previous section shows numerical assessments with smooth images and simple groundtruth displacement fields. The present section presents numerical assessments based on speckle pattern images corresponding to typical usage cases in experimental mechanics. BSpeckleRender software2 [START_REF] Sur | Rendering deformed speckle images with a Boolean model[END_REF] is used to render synthetic speckle images of size 500 × 100 pixels: one image corresponds to a reference state and another one corresponds to a deformed state. Speckles are designed to mimic real patterns such as the ones that can be seen in [START_REF] Lavatelli | Closed-loop optimization of DIC speckle patterns based on simulated experiments[END_REF][START_REF] Su | Elimination of systematic error in digital image correlation caused by intensity interpolation by introducing position randomness to subset points[END_REF][START_REF] Su | Quality assessment of speckle patterns for DIC by consideration of both systematic errors and random errors[END_REF] for instance. The displacement field u = (u x , u y ) is given by the following closed formula:

u x = 0 u y = 0.5 cos 2π y-50 5+45x/500 (41)
The speckle pattern is deformed along the y-direction, the prescribed displacement being a sine wave along the x-direction of maximum amplitude of 0.5 pixel, whose period ranges from 5 pixels (x = 0) to 50 pixels (x = 500). Such a deformation field is relevant in order to highlight the frequency response of the SG filter, see [START_REF] Fayad | Spatial DIC errors due to pattern-induced bias and grey level discretization[END_REF][START_REF] Grédiac | A robust-to-noise deconvolution algorithm to enhance displacement and strain maps obtained with local DIC and LSA[END_REF][START_REF] Grédiac | Comparing several spectral methods used to extract displacement fields from checkerboard images[END_REF][START_REF] Grédiac | On the optimal pattern for displacement field measurement: random speckle and DIC, or checkerboard and LSA?[END_REF] for instance. A 12-bit quantization is used, following the prescription of Section 3.1.2. Note that these images are smaller than in assessment datasets used in some recent works [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF][START_REF] Grédiac | Extracting displacement and strain fields from checkerboard images with the Localized Spectrum Analysis[END_REF][START_REF] Grédiac | Comparing several spectral methods used to extract displacement fields from checkerboard images[END_REF][START_REF] Grédiac | On the optimal pattern for displacement field measurement: random speckle and DIC, or checkerboard and LSA?[END_REF] because of the time needed to render 12-bit images, and in order to facilitate the reading of the graphs. Because of the discrete nature of the rendered speckle images, interpolation is also required. Bicubic interpolation is used, following Section 3.1.1. Consequently, estimation error is likely to be only caused by the ((L u ) T L u ) -1 (L u ) T E term in Corollary 1.

Figure 5 shows two speckle images and the imposed (ground truth) displacement field. Images (before and after deformation) with the large speckle pattern are discussed in Section 3.2.2. Section 3.2.3 concerns images with the fine speckle pattern. The goal here is to compare the displacement field retrieved with SSD minimization to its counterpart given by the predictive formula for the pattern-induced bias (PIB). For didactic purpose, several values are tested for both the subset size and the order of the shape functions. In the remainder of this section, we show the displacement maps (in each direction) retrieved by SSD minimization, the maps predicted by Theorem 1, the ground-truth (GT) displace- ment, and the output of the Savitzky-Golay filtering of the GT displacement. We also show a cross-section plot of the displacement field along its middle-line y = 50 where GT displacement is constant, either null along the x-direction, or equal to 0.5 pixel along the y-direction.

Influence of the subset size

Figure 6 shows the evolution of the difference between the retrieved displacement and the output of the Savistzky-Golay filter as a function of the subset size, for zero-order shape functions used in sterescopy. The results (not shown) for first-order shape functions are similar. The fine speckle pattern images are used here. As discussed in the following sections, the predictive formulas for the retrieved displacement are quite accurate. These figures illustrate the convergence of the retrieved displacement towards the output of the SG filter, as predicted in Section 2.3. It also illustrates that PIB shows a large amplitude over areas where the gradient of the displacement has a large value (here, on the left of the u y displacement field). We can also see that large subsets give very smooth displacement fields, the corresponding lowpass SG filter having the same support as the subset (see [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF][START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]). Such large subsets are not used in practice. Small subsets are required to avoid a large Savitzky-Golay smoothing, but they also give a potentially large PIB.

Large speckle pattern

Figure 7 shows the results for a subset Ω of size 7 × 7 and zero-order shape functions (that is, a constant displacement is estimated over each subset). We can see that the seemingly random fluctuations in the retrieved displacement fields are actually caused by PIB (thus not by sensor noise) and are well predicted by our formulas. In particular, the cross-section plots of the displacement maps confirm that the seemingly random fluctuations of the estimated displacement along the SG filtering of the GT displacement are caused by PIB, the retrieved and predicted curves being superimposed. Interestingly, despite the null GT displacement along the x-direction, the x-component of the retrieved displacement is still affected by PIB. We can see that the amplitude of the PIB is quite large compared to the true displacement. It can be noticed that PIB is large where displacement gradient is large (that is, on the left-hand side of the displacement field). Since the retrieved displacement is a weighted mean of the true displacement in the considered subset, PIB is indeed likely to be larger if the displacement strongly varies within the subset. Figure 8 shows results for 13 × 13 subsets and first-order shape functions. It turns out that higher-order shape functions require larger subsets, since numerical issues affect smaller subsets (not shown here). With first-order shape functions, Λ has 6 components. The first two components correspond to the displacement φ(x) at the center of the subset Ω x and can thus be compared to the convolution of the true displacement with SG filter. The four other components correspond to partial derivatives of the displacement field. They also correspond to the output of (other) SG filters, as recalled in Section 2.3 3 . In all cases, we can see fluctuations of the estimated displacement around the output of the SG filter, these fluctuations being due to PIB. Because of the transfer function of the SG filter, the retrieved displacement field vanishes for certain values. It has even a wrong sign, as discussed in detail in [START_REF] Grédiac | A critical comparison of some metrological parameters characterizing local digital image correlation and grid method[END_REF]. It is quite surprising that the PIB give such large spurious measurements. Figures 9 and10 show results for 19 × 19 subsets and second-order shape functions. Even if second-order shape functions are rarely used in commercial DIC software programs for experimental mechanics, it can be seen that predictive formulas are still valid in this case. estimate the image gradients, we can see that PIB is still accurately predicted. With the same shape functions and subset sizes, we can see that the fine speckle pattern gives a PIB with a smaller amplitude, as discussed in Section 2.3.

Fine speckle pattern

A remark on checkerboard patterns

While random speckle patterns have a huge popularity in photomechanics, very recent papers show that checkerboard patterns give less random noise in the retrieved displacement, see [START_REF] Bomarito | Increasing accuracy and precision of digital image correlation through pattern optimization[END_REF][START_REF] Grédiac | On the optimal pattern for displacement field measurement: random speckle and DIC, or checkerboard and LSA?[END_REF]. The reason is that the average gradient norm within a subset has larger values than with any typical random speckle patterns. Figure 15 shows a checkerboard pattern of pitch equal to 6 (that is, it is made of juxtaposed black and white squares of width 3 pixels). Such a pattern has been used in the experimental assessment of [START_REF] Grédiac | Extracting displacement and strain fields from checkerboard images with the Localized Spectrum Analysis[END_REF][START_REF] Grédiac | On the optimal pattern for displacement field measurement: random speckle and DIC, or checkerboard and LSA?[END_REF]. In these papers, less spurious fluctuations have been observed in the displacement with checkerboard than with random speckle. This is confirmed by Figures 16 and 17 which show estimation of the Λ parameters with checkerboards patterns, with the same settings as in Figures 78with large speckle and 11-12 with fine speckle. While predictive formulas are not fully satisfied here, the curves being not perfectly superimposed, it can be noted that the amplitude of the spurious fluctuations is much smaller than with speckle patterns. The amplitude of the high-frequency spurious displacement on the x-component of the displacement is less than 5 10 -3 , an order of magnitude smaller than with speckles. This high-frequency phenomenon is caused by aliasing; it cannot be seen with large checkerboard pitches (not shown here). Moreover, with such a fine pattern, the numerical estimation of the gradient needed in the predictive formulas is certainly not consistent. While spurious fluctuations cannot be seen, contrary to the fluctuations caused by random speckle patterns, the PIB still plays a role, giving displacement fields which do not fit the output of the SG filter, as can be seen in the cross-section plots in Figures 16 and17 This experiment illustrates that a periodic and fine pattern such as a checkerboard gives a smaller difference with the output of the SG filter than the classic speckle patterns.

Towards a PIB-free SSD criterion?

This section discusses the extent to which it is possible to get rid of PIB, after Blanchet et al.'s approach [START_REF] Blanchet | Fattening free block matching[END_REF]. We can see from Section 2.3.3 that if φ minimizes the following GT-SSD criterion:

x i ∈Ωx 1 max(g 2 i cos 2 (θ i ), κ) I(x i ) -I (x i + φ(x)) 2 (42)
(with the notations of Section 2.3.3), then it is also the least-squares estimate of u over Ω. In this section, κ = 10 -2 for the large speckle pattern, and κ = 10 -4 for the fine one. Figure 18 is representative of results obtained with the large speckle pattern. We can see that the displacement retrieved with this criterion is much less impaired by the spurious fluctuations caused by PIB. Here, it is possible to estimate g 2 i cos 2 (θ i ) since the GT displacement is known. However, Figure 19 shows that this approach is much less efficient with the fine speckle pattern. The reason is that the GT-SSD criterion requires an estimation of the gradient g i and of the angle θ i , which is less accurate with the fine speckle than with the smooth softly-varying speckle pattern. For the same reason, similar results are obtained with the checkerboard pattern discussed in Section 3.2.4.

We now discuss two approaches to PIB-free estimation which do not require the knowledge of the true displacement u. In the present experiment, u is actually a 1-D displacement along the x-axis. It is thus possible to use Blanchet et al.'s approach to fattening-free block matching, and estimate φ over Ω x by minimizing (again with Matlab's fminunc function) the following 1D-SSD criterion:

x i ∈Ωx I(x i ) -I (x i + φ(x)) 2 max(| I x (x i + φ(x))| 2 , κ) (43)
where I x denotes the partial derivative along x of the interpolated image I , an initial guess of φ being given by the classic SSD criterion. The result is shown in Figures 20 and21 (to be compared to Figures 18 and19). PIB significantly decreases: 1D-SSD gives a displacement that roughly follows the output of the SG filter.

Nevertheless, realistic displacements in experimental mechanics are bidimensional. We also introduce and test the following 2D-SSD criterion:

x i ∈Ωx I(x i ) -I (x i + φ(x)) 2 max(|∇ I (x i + φ(x))| 2 , κ) (44)
Figure 22 is representative of results generally obtained: this approach only marginally allows us to decrease PIB, contrary to the case of 1D-SSD. A slight decreasing of the amplitude of PIB can be noticed. As noted earlier in the 1-D case, displacement information is lost only at points where the derivative of the image vanishes. On the contrary, in the 2-D case the aperture problems manifests itself at any pixel since the component of the displacement orthogonal to the gradient always vanishes. This probably leads here to an incomplete correction of the criterion.

Conclusion and open questions

This paper discusses several bias sources in image registration with local parametric estimation via a sum of squared differences criterion, with a focus on photomechanics. A predictive formula is proposed in Theorem 1 for biases caused by interpolation and undermatched shape functions. These sources of errors all depend on the gradient distribution in the underlying images. In particular, the pattern-induced bias (PIB), known as fattening effect in stereoscopic imaging, is caused by the gradient distribution and by the difference between the true displacement field and its local approximation by shape functions. In addition to these biases, the retrieved displacement is affected by spatially correlated random fluctuations caused by sensor noise propagation. This point is discussed in [START_REF] Sur | On biases in displacement estimation for image registration, with a focus on photomechanics -Extended version[END_REF].

Several results from the literature are extended or presented in a unifying way. In Section 2.3, we have also completed the contribution of [START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF] by establishing a rigorous link between the estimated displacement and the true displacement through the Savitzky-Golay (SG) filter, whose parameters depends on the order of the shape functions and on the size of the analysis subset. The link holds because of the very random nature of speckle images.

A numerical assessment of the predictive formulas is discussed as well. First, we have noticed that bicubic interpolation gives biases well below PIB in amplitude, in contrast to bilinear interpolation. Second, it is shown that PIB may have a large amplitude, either if zero-order (used in stereo-imaging) or first and second-order (used in DIC for experimental mechanics applications) shape functions are used. This bias term gives fluctuations around the output of the SG filter of the true displacement. It is striking to note that these fluctuations may have an amplitude of twice the true displacement. As mentioned in Section 2.3, the effect of PIB is equivalent to a spatially varying convolution. A fine understanding of this question would permit to go beyond the deconvolution procedure that is proposed in [START_REF] Grédiac | A robust-to-noise deconvolution algorithm to enhance displacement and strain maps obtained with local DIC and LSA[END_REF] in order to reduce the measurement bias.

In experimental mechanics applications, defining a marking pattern which is optimal with respect to relevant metrological criteria still remains an open question. Some guidelines are given in Section 3.2.5. In addition, first results discussed in [START_REF] Grédiac | On the optimal pattern for displacement field measurement: random speckle and DIC, or checkerboard and LSA?[END_REF] show that checkerboards give a smaller measurement bias than classic speckle patterns used in DIC, and suggest that checkerboards indeed give lower PIB than random speckle patterns. This was numerically verified and illustrated in Section 3.2.4.

Besides, the set of shape functions used to parameterize the local displacement also plays a role: we have shown that PIB involves the scalar product of the image gradient and the difference between the true and the retrieved displacements. It could be of interest to select shape functions minimizing PIB, by improving statistical criteria such as the ones discussed in [START_REF] Bouthemy | Robust model selection in 2D parametric motion estimation[END_REF].

Concerning modifications of the of squared differences criterion permitting to get rid of pattern-induced bias, we have verified that the approach proposed by Blanchet et al. [START_REF] Blanchet | Fattening free block matching[END_REF] is an effective method for 1-D displacements met in stereo-imaging. Our numerical experiments show that it is also a valid approach with higher-order shape functions than the

  ) , DI = I -I , D∇I = ∇I -∇ I , I denotes the interpolation of I at non-integer coordinates, and δ = uφ.

Figure 1 :

 1 Figure 1: Synthetic image in deformed state defined by Equation 39 (image in reference state defined by Equation 38 is not discernible to the naked eye).

Figure 2 :

 2 Figure2: For displacement a) (constant displacement of 0.2 pixel): 1) closed-form expression (first column), 2) bilinear interpolation (second column), and 3) bicubic interpolation (third column). In each of these three cases, the first row depicts, superimposed on the same graph: the ground truth u, the displacement φ(x) estimated from SSD, the displacement φ(x) estimated from SSD, the displacement predicted by Theorem 1, and Savitzky-Golay filtering of u. In each of the three cases, the second row depicts, superimposed on the same graph : the differences (bias estimations) between u on the one hand, and displacement retrieved with SSD, predicted displacement, displacement retrieved with SSD, output of the SG filter on the other hand.

Figure 3 :

 3 Figure3: For displacement b) (sine wave of amplitude 0.2 pixel and period 90 pixels): 1) closed-form expression (first column), 2) bilinear interpolation (second column), and 3) bicubic interpolation (third column). In each of these three cases, the first row depicts, superimposed on the same graph: the ground truth u, the displacement φ(x) estimated from SSD, the displacement φ(x) estimated from SSD, the displacement predicted by Theorem 1, and Savitzky-Golay filtering of u. In each of the three cases, the second row depicts, superimposed on the same graph: the differences (bias estimations) between u on the one hand, and displacement retrieved with SSD, predicted displacement, displacement retrieved with SSD, output of the SG filter on the other hand.

Figure 4 :

 4 Figure4: Effect of gray-level quantization. From left to right: 8-bit images, 10-bit images, 12-bit images, 14-bit images. First row: displacements. Second row: biases (difference between the predicted or retrieved displacements and the ground truth displacement).

Figure 5 :

 5 Figure 5: Two speckle images (size: 500×100 pixels) and y-component u y of the prescribed displacement field. There is no displacement along the x-direction (u x = 0). The first image shows the large speckle pattern discussed in Section 3.2.2, the second one the fine speckle pattern of Section 3.2.3.

Figure 6 :

 6 Figure 6: Fine speckle, zero-order shape functions. From left to right: 7 × 7, 19 × 19, 31 × 31, and 43 × 43 subset Ω. On the top, displacement following x-axis; on the bottom, displacement following y-axis (difference between retrieved displacement field and output of SG filter, and cross-section along the line y = 50).

Figure 7 :

 7 Figure 7: Large speckle, 7 × 7 subset Ω and 0-order shape functions. Displacement along x-(left) and y-axis (right), and corresponding cross-section plots.

Figure 8 :Figure 9 :

 89 Figure 8: Large speckle, 13×13 subset Ω and first-order shape functions. The cross-section plots show, from left to right and top to bottom, λ 1 , λ 3 , λ 5 , λ 2 , λ 4 , λ 6 , such that the xcomponent of the displacement over a subset (giving fields shown on the upper left) is λ 1 + λ 3 x + λ 5 y and the y-component (on the upper right) is λ 2 + λ 4 x + λ 6 y.

Figures 11 to 14 Figure 10 :

 1410 Figures 11 to 14 show the results of the same experiments as in the preceding section, but with the fine speckle pattern. Although such a pattern makes it difficult to reliably

Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 11121314 Figure 11: Fine speckle, 7 × 7 subset Ω and 0-order shape functions. Displacement maps along x-(left) and y-axis (right), and cross-section plots.

Figure 15 :

 15 Figure 15: Checkerboard pattern.

  . The difference caused by PIB reaches, however, values smaller than with random speckle patterns.

Figure 16 :

 16 Figure 16: Checkerboard, 7 × 7 subset Ω and 0-order shape functions. Displacement maps along x-(left) and y-axis (right), and cross-section plots.

Figure 17 :Figure 18 :

 1718 Figure 17: Checkerboard, 13 × 13 subset Ω and first-order shape functions. Displacement maps along x-(left) and y-axis (right). Cross-section plots, from left to right and top to bottom: λ 1 , λ 3 , λ 5 , λ 2 , λ 4 , λ 6 .

Figure 19 :

 19 Figure 19: Fine speckle, 19 × 19 subset Ω, zero-order shape functions, and displacement estimation (along x-on the left and y-axis on the right) with the GT-SSD criterion.

Figure 20 :

 20 Figure 20: Large speckle, 19 × 19 subset Ω, zero-order shape functions, and displacement estimation (along x-axis) with the 1D-SSD criterion. To be compared to Figure 18.

Figure 21 :

 21 Figure 21: Fine speckle, 19 × 19 subset Ω, first-order shape functions, and displacement estimation (along x-axis) with the 1D-SSD criterion. To be compared to Figure 19.

Figure 22 :

 22 Figure 22: Large speckle, 19 × 19 subset Ω, zero-order shape functions, and displacement estimation (along x-on the left and y-axis on the right) with the 2D-SSD criterion.

We use BSpeckleRender b which renders speckle images with patterns of intensity at pixel x varying as exp(-4|x-x0| 2 /R 2 ) with a center x0 given by a Poisson point process and a random radius R, instead of random black disks over a white background, so that the image gradient is a smooth function. Matlab software code is available at the following URL: https://members.loria.fr/FSur/software/BSpeckleRender/

It should be noted that these four components are the derivatives of the local displacement field φ x estimated over Ωx. In photomechanics, the derivatives of the displacement field, which are related to strain components, are rather computed from the derivatives of the global displacement φ.

This work was supported in part by the French National Research Agency (ANR) under grant ANR-18-CE08-0028-01 (ICAReS project).

constant local displacement considered in [START_REF] Blanchet | Fattening free block matching[END_REF]. We have shown that the resulting fatteningfree estimation approximates the output of a SG filter, which generalizes [START_REF] Blanchet | Fattening free block matching[END_REF]. However, the case of 2-D displacements is much more complicated because of the aperture problem, which discards displacement information at any pixel. A PIB-free estimation is still to be designed in this case.

Conflict of interest

The authors declare that they have no conflict of interest.