J.-E Pierré

J.-C Passieux

J.-N Périé

F Bugarin

L Robert

Unstructured Finite Element-based Digital Image Correlation with enhanced management of quadrature and lens distortions

Keywords: Digital Image Correlation, Finite Element, Distortion, Pin-hole

Like subset-based methods, the very first finite element versions of digital image correlation were closely related to the regular structure of images, as they were based on regular quadrilateral elements corresponding to an integer number of pixels. The use of unstructured meshes, to exploit the full potential of FE-DIC in structural mechanics, is now widespread. Most of the time, the formulation, the quadrature and the definition of the region of interest still rely on the pixels grid. In this paper, a formulation in the physical coordinate system and not in the image frame is proposed for 2D digital image correlation. In addition to a more precise definition of the region of interest, it allows the use of a more accurate quadrature rule. It is also shown that lens distortions can be successfully taken into account directly with such a formalism.

Introduction

Digital Image Correlation (DIC) is now extensively used in the field of experimental mechanics for identification and validation purposes. The basic principle of DIC methods is to look for a displacement field that ensures at best the grey level conservation between two images. Most industrial codes are based on subset based approaches [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF]. Since Digital Image Correlation was originally developed by the community of image analysis [START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF], the derived methods were (and are still) closely related to the regular structure of the image. Indeed, original methods rely on a piecewise description of the displacement in several square (or at least rectangular) subsets corresponding to an integer number of pixels, called Zones Of Interest (ZOI ) [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF].

Conversely, many recent developments rather rely on global descriptions of the unknown displacement. The latter is thus defined in the whole Region Of Interest (ROI ). A certain regularity of the displacement field can then be enforced a priori. Several types of decomposition of the displacement field have been proposed. For instance, the reader may refer to [START_REF] Sutton | Measurement of Crack Tip Opening Displacement and Full-field Deformations During Fracture of Aerospace Materials Using 2d and 3d Image Correlation Methods[END_REF][START_REF] Roux | Correlation image velocimetry: a spectral approach[END_REF][START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF][START_REF] Hild | Displacement measurement technique for beam kinematics[END_REF][START_REF] Réthoré | On the Use of NURBS Functions for Displacement Derivatives Measurement by Digital Image Correlation[END_REF][START_REF] Passieux | High resolution digital image correlation using proper generalized decomposition: PGD-DIC[END_REF][START_REF] Passieux | A Digital Image Correlation Method For Tracking Planar Motions Of Rigid Spheres: Application To Medium Velocity Impacts[END_REF], to name a few.

In this paper, we focus on a Finite Element approach to DIC (FE-DIC [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF]). It was shown that this kind of approach has many advantages in the context of continuum solid mechanics and especially for the experiment/simulation dialog [START_REF] Hild | Comparison of Local and Global Approaches to Digital Image Correlation[END_REF]. This kind of tool is particularly appealing in the context of a Finite Element Model Updating approach based on displacement fields (FEMU-U [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF][START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF][START_REF] Gras | Identification of the Outof-Plane Shear Modulus of a 3d Woven Composite[END_REF]). It is then tempting to use the (optimized) simulation mesh to perform the FE-DIC measurements [START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF]. However, the very first developments [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Hild | Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties -a Review[END_REF][START_REF] Fehrenbach | A fast algorithm for image registration[END_REF] were still based on regular quadrilateral elements (denoted Q4) corresponding to an integer number of pixels. The use of unstructured meshes is now widespread [START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF][START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3d kinematics[END_REF][START_REF] Baldi | Comparing Two Damage Models Under Shear Stress[END_REF][START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF][START_REF] Beeck | Quantification of Three-Dimensional Surface Deformation using Global Digital Image Correlation[END_REF][START_REF] Wittevrongel | A Self Adaptive Global Digital Image Correlation Algorithm[END_REF] to take advantage of the full potential of FE-DIC. Most of the previous references still rely on the pixels, since their integration scheme is based on pixels and since the definition of the ROI is also defined by a set of pixels.

Recently, a new formulation of the DIC problem in the physical space has been proposed in the context of stereo-DIC [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF][START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3d kinematics[END_REF]. It is based on a kinematic model (FE mesh, NURBS patch...) in the world coordinate system along with a projective model. This paper investigate the effect of such a formulation in the context of 2D FE-DIC. In practice, the 3D mesh of the specimen skin is used. Its position in space is defined so that its projection(s) on the reference image(s) matches at best different amers. The idea is then to search for the displacement, directly expressed in the physical space, minimising a continuous problem associated to the grey level conservation. Images are only considered as simple data. The algorithm compares, for any point in any element, the grey levels corresponding to its projection in different images. This paper shows that such a formulation is of interest for 2D FE-DIC for the two following reasons:

1. the assembly of the numerical operators requires an appropriate quadrature rule. In opposition with standard approaches, for which the reference image is used to define both the integration technique and the ROI, here, a dedicated surface quadrature technique in the world reference system is defined. The efficiency of those two integration methods will be compared; 2. as for stereo DIC algorithms, distortions can be taken into account during the measurement. Even if it seems less natural to resort to calibration in this context, [START_REF] Lava | Impact of lens distortions on strain measurements obtained with 2d digital image correlation[END_REF] recently shown that it is worth taking distortions into account even for 2D analysis.

The effect of lens distortions will be considered here but not the one due to out-of-plane motions [START_REF] Yoneyama | Lens distortion correction for digital image correlation by measuring rigid body displacement[END_REF][START_REF] Sutton | The effect of out-ofplane motion on 2d and 3d digital image correlation measurements[END_REF][START_REF] Pan | High-Accuracy 2d Digital Image Correlation Measurements with Bilateral Telecentric Lenses: Error Analysis and Experimental Verification[END_REF][START_REF] Pan | High-accuracy 2d digital image correlation measurements using lowcost imaging lenses: implementation of a generalized compensation method[END_REF]. In general, the latter cannot be taken into account in 2D-DIC without strong kinematic assumptions. For instance, it is possible when considering a rigid-body motion of the observed object [START_REF] Passieux | A Digital Image Correlation Method For Tracking Planar Motions Of Rigid Spheres: Application To Medium Velocity Impacts[END_REF][START_REF] Réthoré | 3d displacement measurements using a single camera[END_REF]. With a single camera, it is difficult to distinguish the effects of out-of-plane motions and strains as it is explained in [START_REF] Sutton | The effect of out-ofplane motion on 2d and 3d digital image correlation measurements[END_REF].

In this article, a synthetic image, subjected to an in plane (large) translation, is considered. The 2D displacement field is measured with a given mesh and a calibration model.

Section 2 presents the usual quadrature method. This method considers a Region of Interest (ROI ) which is simply a part of the reference image. Section 3 introduces a mesh based method quadrature, better suited for a Finite Element DIC Method. The two previous approaches are compared in Sect. 4 in terms of measurement uncertainties. At this stage, a simple projector is used. Finally, in Sect. 5, the way distortions can be taken into account into the mesh based integration is presented. A pin-hole model of camera is considered in order to evaluate the gain with respect to measurement uncertainties.

Quadrature at the pixels

Formulation

The optical flow problem between the reference state image f and the deformed state image g reads:

f (x) = g(y) (1)
where x and y are the pixels positions. The latter problem is usually written in a weak form over a region called Region of Interest (ROI ). The ROI is often defined as part of an image, and more precisely, a part of the reference state image f . In this equation, y = x + u(x) and the problem reads:

u = argmin u∈L 2 (ROI) ROI [f (x) -g(x + u(x))] 2 dx (2)
with g the deformed state image with corrected mean value and standard deviation, in order to take into account offset and scale changes in illumination [START_REF] Chambon | Dense matching using correlation: new measures that are robust near occlusions[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements -Basic Concepts,Theory and Applications[END_REF][START_REF] Pan | Equivalence of digital image correlation criteria for pattern matching[END_REF][START_REF] Passieux | A Digital Image Correlation Method For Tracking Planar Motions Of Rigid Spheres: Application To Medium Velocity Impacts[END_REF]. This problem being non-linear, a first order Taylor expansion of g is performed with respect to the displacement correction δu = u-u 0 given an initial estimate u 0 . Therefore, the linearised problem reads [START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF][START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF]

u = argmin u∈L 2 (ROI) ROI ∇f T δu(x) -r(x) 2 dx (3)
where r(x) = f (x) -g(x + u 0 (x)) is the gray level residual field and ∇f the gradient of the reference image f estimated by finite differences. This problem has got a variational form: Find δu ∈ L 2 (ROI), such that:

a(δu, v) = l(v) ∀v ∈ L 2 (ROI) (4)
with

       a(u, v) = ROI u (x)∇f (x)∇f (x)v(x) dx l(v) = ROI v (x)∇f (x) r(x) dx (5)
In this work, the interpolation of the unknown displacement δu is based on finite elements [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF] such that δu(x) = i N i (x) q i with q i the degrees of freedom and N i the shape functions. Then, the linearised DIC problem reads:

M DIC q = b DIC with: M DIC ij = a(N i , N j ) b DIC i = l(N i ) (6)
For a sufficiently regular image gradient, it can be shown that the correlation matrix M DIC is symmetric definite positive [START_REF] Fedele | Global 2d digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF]. Another interesting result is that M DIC only depends on the initial data, so its computation is done only once for each time step.

Quadrature

In this section, a classical quadrature rule based on the image is explained to compute M DIC and b DIC . Indeed, the integrals of equations ( 6) are approximated by:

ROI ψ(x)dx ≈ e∈Ne k∈Ie ω e k ψ(x e k ) (7)
where x e k and ω e k are the quadrature points and weight coefficients of element e respectively; N e being the number of elements. I e the set of quadrature points x e k is detailed here. As stated previously, the Region of Interest (ROI ) is given here by a set of pixels. The nodes of element e are first projected in the reference image and the element is reconstructed in the image plane. For instance, with linear elements, it means that the edges are assumed to remain straight in the image plane (this assumption is actually not true when considering lens distortions, which is the object of Sect. 5). In this case, I e is the set of pixels whose center belongs to element e and x e k is the center of the pixel. The associated weight coefficients are the area of one pixel, namely ω e k = 1. Since x e k is a pixel of the reference image f , this means that f is also the reference space for the integration. Finally, the computation of the left hand side of equation ( 6), for instance, reads:

a(N i , N j ) ≈ e∈Ne k∈Ie N i (x e k )∇f (x e k )∇f (x e k )N j (x e k ) (8)
The position of the quadrature points being known in the real coordinate system, another issue concerns the evaluation of the shape functions N(x e k ). Indeed, in expression (8), the shape functions N i has to be evaluated at pixel x e k , knowing its coordinates in the global system (x, y) (cf. Fig. 2). But with the Finite Element Method, an iso-parametric representation of the element is used, and the shape functions N are known in the reference coordinate system (ξ, η) [START_REF] Wittevrongel | A Self Adaptive Global Digital Image Correlation Algorithm[END_REF].

These basis functions are used to describe the unknown displacement field but also the geometry of the element. Indeed, knowing the real positions s i = (x i , y i ) of the vertices of element e, the mapping between the coordinates (ξ, η) of a point x 0 in the reference configuration and the coordinates (x, y) of the corresponding point x = φ(x 0 ) in the real configuration is given by [START_REF] Zienkiewicz | The finite element method[END_REF] (cf. Fig. 2).

x

= i N i (ξ, η)s i (9)
To compute the elementary DIC operators, the inverse of this mapping has to be determined. Indeed, the coordinates (x p , y p ) of a pixel x are known in the real configuration, and the coordinates (ξ p , η p ) of its corresponding point x 0 are searched for. This is performed, in practice, thanks to few iterations of a Newton algorithm on φ(ξ p , η p ) = x -N i (ξ p , η p )s i Remark. For the calculation of r(x k ), f (x k ) defines the gray level of pixel x k in the reference state image. In opposition, x k +u(x k ) may not be an integer number of pixels and thus requires to use a gray level interpolation (usually bilinear, bi-cubic or spline) for computing g(x k +u(x k )).

The limitations of such an approach are as follows: • the quadrature rule is not exact even for computing the quadrature of a constant function. This means that the area of an element is not exact (cf. Fig. 3). This yields measurement uncertainties;

• the geometric assumptions are required which may be relevant enough with telecentric lenses, but not in the general case. This is particularly important when quantitative comparison with numerical simulation are searched for.

In the following section, an alternative is presented which does not have these limitations.

Figure 3: With the classical method the workspace is part of an image and there are some approximations (the computed area, in red, is not exact).

3. Mesh Based Quadrature

Formulation

The idea is to compute the integration in the physical space defined by the mesh and not in the ROI positioned in the image. This method considers a projection model, denoted P , which maps the coordinates X (in meter) of a point in the physical space to the coordinates x (in pixel) of its projection onto the image plane. Different examples of projector will be detailed in Sections 4 and 5.

Assuming that P is known, it is possible to write the gray level conservation equation (1) with:

x = P (X)

y = P (X + U(X)) (10)
where the unknown is the displacement defined in the physical space U. Then, the weak form of the optical flow problem ( 1) is written as an integral over Ω which is the domain defined by the model in the physical space and not in the image:

U = argmin U∈L 2 (Ω) Ω [f (P (X)) -g(P (X + U(X)))] 2 dX (11)
Remark. A link between the ROI in the standard FE technique and the 3D domain Ω could be ROI = P (Ω). But, as it will be shown in Sections 4 and 5, in practice, with classical FE-DIC techniques, this equality may not be perfectly verified.

The problem is again linearised by considering a displacement increment δU = U -U 0 . This unknown is searched for in an approximation subspace defined in the physical space by the FE interpolation:

δU(X) = i N i (X)q i
The variational form thus becomes:

Find u ∈ L 2 (Ω) : a(u, v) = l(v) ∀v ∈ L 2 (Ω) with: (12)        a(u, v) = Ω u ∇P ∇f ∇f ∇P v dX l(v) = Ω v ∇P ∇f r(X) dX (13)
where the residual map r(X) and the gradients of the projector ∇P and the reference image ∇f read:

     r(X) = f (P (X)) -g(P (X + U 0 (X))) ∇P = ∇P (X + U 0 (X)) ∇f = ∇f (P (X)) (14)
And the correlation system to be solved iteratively is in the same form as before:

Mq = b with: M ij = a(N i , N j ) b i = l(N i ) (15)
It can be seen that the correlation matrix depends on the gradient of the projector (cf. eq. [START_REF] Gras | Identification of the Outof-Plane Shear Modulus of a 3d Woven Composite[END_REF]. A priori, this operator may depend on the approximation of the unknown U 0 (X). As a consequence, correlation matrix M may also depend on it, which implies that it may be rebuild at each iteration. A discussion on this point will be given in Sections 4 and 5 depending on the type of projector. Furthermore, the computation of each term of equation ( 14) requires a gray level interpolation of both reference and deformed images because of the projector: both P (X) and P (X + U 0 (X)) may not be integer numbers of pixels.

Quadrature

To compute the integrals of M and b in the physical space, a quadrature rule based on the mesh is required. As explained previously, the integral of, for instance, the left hand side of equation [START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF], is approximated by:

a(N i , N j ) ≈ e∈Ne k∈Ie N i (X e k )∇P (X e k )∇f (P (X e k )) ∇f (P (X e k ))∇P (X e k )N j (X e k ) (16)
where I e remains a set of quadrature points X k associated to weight coefficients ω e k . The difference with Sect. 2, is that the quadrature points are defined in the physical space. Following what is done in FE simulation, the integrals are estimated by:

Ω f (P (X))dX = φ -1 (Ω) f (P (φ(X 0 )))J φ dX 0 (17)
where X 0 is a point of the iso-parametric element and J φ the Jacobian determinant of the transformation between the iso-parametric and deformed elements in the physical space (cf. Fig. 2). In other words, it means that the mesh is no longer projected onto the reference image, as it was the case before. Only the quadrature points will be projected. The positions of those points are far more precise because the elements are not artificially rebuilt from the projection of the vertices of the elements. The interest quantities are searched for in the images but the computation stay in the physical space.

Up to now, the choice of the quadrature rule (X k and ω e k ) was not precised. In FE simulation, Gauss-Legendre quadrature is used since it is known to be optimal for the integration of polynomial functions. In opposition, in DIC, this choice is not relevant, because the gradient of the image ∇f is short varying. In the following, a more simple rectangle-like quadrature rule based on evenly placed points will be used. It has the advantage to, at least, be exact for the computation of the area of an element as it is illustrated in Fig. 4 which is not the case with classical pixel integration methods (cf. Fig. 3). In order to use exhaustively the image data, a sufficiently large number of physical integration points have to be considered in order to have, at least, (the projection of) one integration point per pixel in the image. For that, the following procedure was adopted: each node is projected in the reference image, then the maximal number of pixels between two nodes (the diameter of the element) is determined for each element. Then, the reference element is subdivided into nested homothetic sub-triangle or quadrangles. Thus, the quadrature points X are the center of gravity of these sub-triangles or quadrangles and the weight coefficients ω e k are the associated areas which are constant over a given element. (cf. Fig. 5). Remarks.

1. The positions and the number of integration points obviously depend on the mesh. If the mesh is refined, the number and the positions of the mesh could change. 2. The number of integration points (which is more or less the number of pixels for standard 2D-DIC) has a direct impact on the computational cost, in particular when using high definition sensors. In the present paper, since a front-parallel setup is considered, the number of integration points is almost the same than the number of pixels. This implies that the computational cost is merely unchanged. Nevertheless, if needed, some methods have been proposed to address the issue of numerical efficiency, such as [START_REF] Champagnat | Fast and accurate PIV computation using highly parallel iterative correlation maximization[END_REF][START_REF] Passieux | High resolution digital image correlation using proper generalized decomposition: PGD-DIC[END_REF][START_REF] Passieux | A Digital Image Correlation Method For Tracking Planar Motions Of Rigid Spheres: Application To Medium Velocity Impacts[END_REF][START_REF] Passieux | A dual domain decomposition method for finite element digital image correlation[END_REF].

In the next sections, it will be shown that this choice (using at least one integration point per pixel) is a good compromise between measurement uncertainties and computation cost. This formulation in the physical space and associated quadrature rule will be referred to as Mesh Based Formulation in the following.

A first perfect projector

Description of the projective model

In this section, a simple linear and constant projector is considered. It is based on a frontparallel camera model1 without distortions but with coincident spaces. Thus, the projector only depends on the scale factor k and is written:

x = P (X) = N f -kY kX ( 18
)
where N f denotes the number of pixel of the image in the vertical direction. The scale factor k corresponds to the number of pixels per meter: the scale transformation between the physical position and the coordinates in pixels. As explained above, the correlation operator M involves the spatial gradient of this projector. With such a simple model, the gradient is constant which means that the correlation matrix does not need to be rebuild at each iteration:

∇P = 0 -k k 0 (19)
This simple model of projection makes the Mesh based formulation very close to the classical pixel formulation. Only the quadrature rule differs in this case. In the next section, the effect of the quadrature on measurement uncertainties is analysed.

Effect on measurement uncertainties

In order to compare both methods, a real image (cf. Fig. 6) was taken as reference. A set of synthetic deformed state images were generated by prescribing a pure translation of magnitude 0.3, 0.4 and 0.5 pixel in the horizontal direction using image interpolations based on the shift / modulation property of the Fourier transforms. For this analysis, several unstructured meshes were used with different sizes of elements ranging from 4 to 30 pixels per element edges. Both linear quadrangles and triangles were used for comparison purposes. The same meshes are used to compare the quadratures rules. For each size and type of elements of the mesh, and for each deformed state images, a displacement measurement is performed. Then the standard deviation and the bias of the displacement error (defined by the difference between the prescribed and measured displacement fields) are computed and their values are averaged over the 3 shift magnitudes [START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF]. The results are plotted in Fig. 7. The displacement standard uncertainties are comparable even if it seems that the proposed quadrature rule gives slightly better results (cf. Fig. 7). In opposition, it is clear that the pixel integration yields a systematic error which is reduced by one order of magnitude when using the proposed more accurate quadrature rule.

Effect of the number of integration points

Previously, the number of integration points was computed in order to have at least one point per pixel. In this section, the choice of the number of quadrature points of the proposed mesh based formulation is discussed. The mesh is fixed with an mean element size of 20 pixels. In Fig. 8, the effect of the number of quadrature points on the measurement uncertainties is analysed as a function of the ratio of the average number of points to the average number of pixels by element in percent. It can be seen that the choice proposed previously to put approximately as many integration points as the number of pixels in each element is optimal. Indeed, it allows to consider the minimum set of integration points without affecting the measurement uncertainty. Taking more points would lead to useless extra computational costs, whereas using less points would lead to higher measurement uncertainties.

Application to an enhanced projector taking into account lens distortions

One of the main advantages of using such a formulation written in the physical space is that it is possible to use a more complex camera model, for instance including the actual distortions of the lens directly during measurement. Indeed, taking into account distortion is not new to DIC [START_REF] Zhang | Displacement/strain measurements using an optical microscope and digital image correlation[END_REF][START_REF] Sutton | Metrology in a scanning electron microscope: theoretical developments and experimental validation[END_REF][START_REF] Sun | AFM image reconstruction for deformation measurements by digital image correlation[END_REF][START_REF] Yoneyama | Lens distortion correction for digital image correlation by measuring rigid body displacement[END_REF][START_REF] Yoneyama | In-plane displacement measurement using digital image correlation with lens distortion correction[END_REF][START_REF] Sutton | Scanning Electron Microscopy for Quantitative Small and Large Deformation Measurements Part I: SEM Imaging at Magnifications from 200 to 10,000[END_REF][START_REF] Sutton | Scanning Electron Microscopy for Quantitative Small and Large Deformation Measurements Part II: Experimental Validation for Magnifications from 200 to 10,000[END_REF][START_REF] Lava | Impact of lens distortions on strain measurements obtained with 2d digital image correlation[END_REF].

In the next section, an example of a camera model is described [START_REF] Garcia | Mesure de formes et de champs de déplacements tridimensionnels par stéréocorrélation d'images[END_REF][START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements -Basic Concepts,Theory and Applications[END_REF][START_REF] Hartley | Multiple view geometry in computer vision[END_REF] and its use in the physical formulation of DIC is explained. Then, a virtual test based on synthetic images is analysed in order to estimate the benefits of including lens distortions.

An enhanced Model of a Camera Linear Model

A general projective transformation is a linear transformation mapping any arbitrary physical world point X w = (X, Y, Z) (denoted by X w with homogeneous coordinates in the projective space P 3 ) into a 2D point x = (u, v) (x in the projective space P 2 ). It consists of three consecutive transformations.

First, a rigid-body transformation (rotation + translation) T transforming a physical point X w in the world coordinate system R w into a physical point X c in the imaging sensor coordinate system R c :

X c ∼ = R t 0 1 T X w ( 20
)
with R the rotation matrix and t the translation vector. The parameters of this transformation T are called extrinsic parameters and represent the orientation and the position of the imaging system in the world. Second, a transformation P mapping a point in P 3 to a point in P 2 . This latter transforms the physical point X c into a 2D point in the retinal plane xr :

xr ∼ =   f 0 0 0 0 f 0 0 0 0 1 0   P X c ( 21
)
where f is the focal length.

Finally, an affine transformation A is applied to transform the 2D point xr into a 2D point x in the distortion-free image plane:

x ∼ =   α x α x cot θ c x + c y cot θ 0 αy sin θ cy sin θ 0 0 1   A xr ( 22
)
where α x and α y are the horizontal and vertical aspect ratio respectively. c x and c y are the horizontal and vertical coordinates of the intersection of the optical axis with the image plane respectively, also referred to as principal point. This latter point should match the center of the image. θ is the angle between the x-axis and the y-axis (cf. Fig. 10), also referred to as skew angle. Usually, the horizontal and vertical axis are almost perpendicular (i.e. θ ≈ π 2 ), then:

A =   α x 0 c x 0 α y c y 0 0 1   . ( 23
)
Note that the two last transformations ( 21) and ( 23) are almost always combined in only one transformation: AP. In order to introduce the parametric distortion model (as explain in the sequel), the latter transformation is usually replaced by a canonical projection Π and the transformation K both defined as following:

Π =   1 0 0 0 0 1 0 0 0 0 1 0   (24)
K =   f x 0 c x 0 f y c y 0 0 1   (25)
where f x = α x f and f y = α y f. K is called intrinsic matrix. Its parameters are called intrinsic parameters and represent the internal parameters of the imaging sensor. Note that the focal lens is encapsulated in K and not in the projection Π. Figure 9 summarizes the projection process from R w to the distortions-free image plane. Since, in this paper, only 2D-DIC is considered, all the images are taken from a front-parallel view. Consequently, -→ Z w (the z-axis of R w ) and -→ Z c (the z-axis of R c ) are parallel, so that Z coordinates can be considered as constant. Fig. 10 summarizes the projection process in the front-parallel case.

Note that the presented projection-perspective camera model, which does not include distortions, is linear in the projective space (i.e. from P 3 to P 2 ). However, due to the canonical projection Π, this model is not linear in the euclidian space (i.e. from R 3 to R 2 ).

Lastly, a tricky point consists in the images storage. Indeed the images f and g are stored as tables (or matrices) of pixels f i,j at row i and column j. Hence, the corresponding coordinate system (whose center is in the top left corner) requires an additional linear transformation composed of a 90 • -rotation and a translation. This latter transformation is given by: 90 •rotation and translation expressed by:

u v = 0 -1 1 0 u v + N f 0 . (26)
In the following, different approaches for taking lens distortions into account are presented.

Distortion Models

Distortion models can be classified in two categories: the parametric models and the nonparametric models. Despite their different approaches, the key point consists in their different locations in the global projection diagram of Fig. 9.

In the parametric model, denoted by D, only three types of distortions are considered: the radial distortion ∆ r (par-axial approximation), the decentralize distortion ∆ d (imperfect alignment of the optical centres of all lens) and the prismatic distortion ∆ p (bad parallelism of the lens between them and/or with the plane of the detector). For more advanced informations, the reader can refer to [START_REF] Garcia | Mesure de formes et de champs de déplacements tridimensionnels par stéréocorrélation d'images[END_REF][START_REF] Besnard | Caractérisation et quantification de surfaces par stéréocorrélation pour des essais mécaniques du quasi statique à la dynamique ultra-rapide[END_REF][START_REF] Bugarin | Vision 3d multi-images: contribution à l'obtention de solutions globales par optimisation polynomiale et théorie des moments[END_REF]. The expression of D, ∆ r , ∆ d and ∆ p are given by:

D : R 2 → R 2 xr → xd = xc + ∆ r (x c ) + ∆ d (x c ) + ∆ p (x c ) (27)
∆ r (x c ) =       x 3 i=1 κ i (x 2 c + y 2 c ) i y 3 i=1 κ i (x 2 c + y 2 c ) i       (28)
∆ d (x c ) = 2d 1 x c y c + d 2 (3x 2 c + y 2 c ) 2d 1 x c y c + d 2 (x 2 c + 3y 2 c ) (29) ∆ p (x c ) = p 1 (x 2 c + y 2 c ) p 2 (x 2 c + y 2 c ) (30)
In this approach [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF][START_REF] Brito | One-sided Radial Fundamental Matrix Estimation[END_REF][START_REF] Dufour | Integrated digital image correlation for the evaluation and correction of optical distortions[END_REF], lens distortions are considered between Π and K (cf. Fig. 11). More precisely, lens distortions are introduced after the projection on the retinal plane x r = (x r , y r ) and before applying the intrinsic parameters. Due to the good quality of the modern lenses, a development at the 3 rd order of the radial distortion ∆ d is enough for most of the camera lens with a focal length greater than 5 mm.

In the non-parametric model [START_REF] Cornille | Mesure précise de formes et de déformations 3d d'objets observés au Microscope Électronique à Balayage (MEB)[END_REF], the main idea is to recover the distortion shape from the remaining reprojection errors of the calibration of a pure projective model and fitting it thanks to an interpolation method (e.g. spline interpolation). This approach allows to represent any type of distortion field, even complex. However, this model, denoted by I, is included before applying the intrinsic matrix K (cf. Fig. 12).

Moreover, for the sakes of simplify and computational efficiency of the distortion removal process, the interpolation I can be replaced by its reciprocal I -1 , the distortion correction (cf. Fig. 13).

Choosing a distortion model is strongly related to the informations available. Indeed, if no calibration informations are available or in case of complex lens, a non-parametric model can be considered. However, if the calibration parameters are available, the parametric model will be more accurate.

Synthetic test case

In this section, the aforementioned non-linear projector is used in the mesh based DIC method described in Sect. 3 for the analysis of a synthetic test case. Like this, it is possible to filter the effect of lens distortions during the measurement.

In practice, a real physical experimentation is performed in order to feed the non-linear projector with realistic intrinsic parameter values. For that, the configuration is as followed: a CCD camera (Pike F-145b, 8 bit, 1388×1036 px) equipped with a Nikon Nikkor 50 mm f2 aperture 11 lens was set up, allowing a field-of-view of about 270×200 mm 2 and an optical resolution roughly equal to 5 px/mm. Then 110 images of a conventional target were taken. The calibration (via a projection-perspective first-order radial model sufficient for this kind of lens) was done using the commercial software Vic-3D TM from Correlated Solutions, Inc. [] for reliable results. Intrinsics were: u 0 = 706.053 px, v 0 = 503.725, k u = k v = 8381.24 px and r 1 = -0.0321253 with a reprojection error of 0.05 px. Note in addition that a second camera was also positioned (stereo pan angle of almost 15 • ) and the results from a stereo-rig calibration are close to those from a single camera calibration.

Then a set of synthetic images is created. It has to be noted that the synthetic images are created using the real distortion parameters that were previously characterised experimentally. The only way to quantify the efficiency of the proposed method in reducing the measurement uncertainty over a standard DIC, is to develop a synthetic test case as classically done in the literature ( [START_REF] Bornert | Assessment of Digital Image Correlation Measurement Errors: Methodology and Results[END_REF][START_REF] Schreier | Systematic errors in digital image correlation due to undermatched subset shape functions[END_REF]) in order to have a ground truth reference. Otherwise, the reference becomes unknown and it would be difficult to conclude. First a realistic synthetic speckle pattern is mapped onto a plane in the physical space, such that a gray level F p (X) is associated to each physical point X. Then, for any integer pixel position in the reference image plane x p , the coordinates X p of the corresponding point in the physical space are obtained by solving the non linear equation x p -P (X p ) = 0, thanks to a Newton algorithm. Then a reference image f is synthesised by estimating the gray level F p (X p ) associated to the physical point X p for each integer pixel position x p . Finally a deformed state image g is synthesised by estimating the gray level F p (X p + S) where S is a rigid body shift in the physical space.

In the following, a rather large shift S = (400, 200, 0) pixels is chosen in order to highlight the effect of distortions. Then a Pixel and a Mesh Based DIC between the distorted images f and g were performed with different element types and sizes. In addition, the effect of the number of quadrature points on the measurement uncertainties is again analysed as a function of the ratio of the average number of points to the average number of pixels by element in percent. Like previously, the measurement uncertainties are quantified as standard deviation and the bias (cf. Fig. 14,15). This was done in order to see that one integration point per pixel is the minimum, and that it is useless to have more information. As expected, all curves from Pixel Integration are widely above those from Mesh Based Integration because in the former case, lens distortions are not taking into account which is crippling for such large displacements. Consequently, effects of the elements size (cf. Fig. 14) or of the number of integration points (cf. Fig. 15) have become of second order in the errors.

It is interesting to see the displacement field measured on a square field of view of 300 × 300 px 2 for a given element size with both methods (cf. Fig. 16). With a classical Pixel Integration, the intrinsic parameters of the camera are not taken into account, which means that the measured displacement is not uniform and the iso-distortion circles can be seen on the displacement field. Thus, the Mesh Based Integration becomes clearly more precise especially if large displacements are considered.

Conclusion

In a classical 2D FE-DIC method, the integration of the gray-level difference is still closely related to the pixels of the reference image. Not only the ROI is an integer set of pixels (the element shapes are approximated), but the quadrature rule also uses pixels x as integration points and their area as weights.

In opposition, in this paper, it is chosen to minimise the gray-level difference in the physical space with a dedicated quadrature method. The region of interest is defined in the physical space by the mesh geometry. The quadrature rule is based on points X in the world coordinate system. A projective model is used to compute the corresponding gray-levels in both reference f (P (X)) and deformed g(P (X)) images. The number of quadrature points X is set such that in any pixel, one can find at least the projection P (X) of one point. It has to be noted that, with this method, the measured displacement is then directly expressed in the physical space.

At first, thanks to the proposed quadrature technique, it is shown that a systematic measurement error is avoided, even for a simple projection model (without lens distortions). Then, an enhanced camera model is used as a projector. In practice, a pin-hole model is implemented in order to consider lens distortions. Even if the latter have much less impact when dealing with small strains, it is shown that they have to be considered when measuring large displacements. Distortions should be taken into account for displacement-based FEMU approaches. The proposed method should allow performing 2D FE-DIC measurements in tricky optical setups, for example when using microscopes [START_REF] Efstathiou | Multiscale strain measurements of plastically deforming polycrystalline titanium: Role of deformation heterogeneities[END_REF], wide-angle lenses or portholes.

With the computed physical displacement, it becomes easier to use a displacement-based FEMU method or to perform a 2D identification [START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF]. Obviously, the proposed formulation will have even a greater interest in the context of Finite Element Stereo (or Multiview) Digital Image Correlation for the measurement of 3D surface displacement.

Figure 1 :

 1 Figure 1: Reference image with a point x (left) and deformed image with the corresponding point x + u(x) (right)

Figure 2 :

 2 Figure 2: Reference element (left) and real element on a mesh (right) for a quadrangle Q4

Figure 4 :

 4 Figure 4: With a mesh based quadrature, the workspace is the domain of definition of the finite elements (the computed area, in red, is exact)

Figure 5 :

 5 Figure 5: Subdivision of an element of an area A: (1) projection of the nodes of the element on the reference images, (2) computing of the diameter d of the element (here, d = 3 pixels), (3) subdivision of the iso-parametric element by subdividing each edges with the diameter, (4) computation of the gravity center of each sub-element. It creates d 2 sub-elements of the same area a = A/d 2 in the real element. This is exactly the same for quadrangles

Figure 6 :

 6 Figure 6: Reference image with an initial mesh with quadrangle (mean element size = 10 pixels)

Figure 7 :

 7 Figure 7: Evolution of the displacement standard uncertainty and the bias with respect to sub-pixel prescribed displacement for different sizes of elements. The results are computed for a Mesh Based Integration (mesh) and compared to the results of the Pixel Integration (pixel ) with quadrangles (qua) or triangles (tri )

Figure 8 :

 8 Figure 8: Evolution of the displacement standard uncertainty and bias with respect to sub-pixel prescribed displacement for different number of quadrature points (mean element size = 20 pixels). The results are computed for a Mesh Based Integration (mesh) and compared to the results of the Pixel Integration (pixel ) with quadrangles (qua) or triangles (tri )

Figure 9 :

 9 Figure 9: Diagram of the projection process from the world coordinate system to the distortions-free image plane.

Figure 10 :

 10 Figure 10: The three elementary transformations of the camera pin-hole model for a front-parallel view.

  Figure 11 summarizes the introduction of the parametric model into the projection diagram.

Figure 11 :

 11 Figure 11: Diagram of the projection process from the world coordinate system to the distortions image plane in the parametric case.

Figure 12 :

 12 Figure 12: Diagram of the projection process from the world coordinate system to the distortions image plane in the non-parametric case (direct interpolation).

Figure 13 :

 13 Figure 13: Diagram of the projection process from the world coordinate system to the distortions image plane in the non-parametric case (inverse interpolation).

Figure 14 :Figure 15 :

 1415 Figure 14: Evolution of the displacement standard uncertainty and the bias with respect to sub-pixel prescribed displacement for different elements sizes on a distorted image. The results are computed for a Mesh Based Integration (mesh) and compared to the results of the Pixel Integration (pixel ) with quadrangles (qua) or triangles (tri )

Figure 16 :

 16 Figure 16: Exact error of the vertical component of the displacement field in pixels measured on distorted images with the Mesh Based Integration (left) with a calibrated camera or with a classical Pixel Integration (right) without taking account the camera's intrinsic parameters.

this model is based on a pinhole model with d = 0 and T = I 4 cf. Sect. 5.1

Acknowledgments

This work was funded by the French "Agence Nationale de la Recherche" under the grant ANR-12-RMNP-0001 (VERTEX project).


Fabien Tarissan
email: tarissan@lix.polytechnique.fr

Camilo La Rota
email: camilo.larota@ens-lyon.fr

Inferring Update Sequences in Boolean Gene Regulatory Networks

Keywords: Mathematical programming, Inverse problems, Gene regulatory networks reconstruction

Introduction

This paper employs mathematical programming and mixed integer linear programming techniques for solving a problem arising in the study of genetic regulatory networks. More precisely, we solve the inverse problem consisting in the determination of the sequence of updates in the digraph representing the gene regulatory network (GRN) of Arabidopsis thaliana in such a way that the generated gene activity is as close as possible to the observed data. Differences among cells of different tissues depend on the specific set of genes that are active in each tissue. Therefore, one usually assumes that the different steady states of a GRN dynamics correspond to the different possible cell fates ( [START_REF] Thomas | Multistationarity, the basis of cell differentiation and memory. i. structural conditions of multistationarity and other nontrivial behavior[END_REF]). This leads to explain the changes observed during the development of the organisms by the fact that perturbations on specific elements of the network make the system switch from one steady state to another. Some hypothesis can be made about these perturbations, which are then treated as initial conditions for the new tissue being formed. However, an important unknown is (are) the update sequence(s) of the gene activity that let the system evolve from a given set of initial conditions to the set of steady states. Indeed, different update sequences determine different sets of basins of attraction of the GRNs. However, the steady states remain the same under any sequence. Usually, a specific update sequence is assumed to rule the dynamics of the GRNs [START_REF] Demongeot | Robustness in Regulatory Networks: a Multi-Disciplinary Approach[END_REF][START_REF] Gershenson | Classification of random boolean networks[END_REF]. The present study differs from this approach in that we sought to infer the update sequence from the biological observations. It also differs from our previous paper as we focus here on asynchronous sequences whereas in [START_REF] Tarissan | Network reconstruction: a mathematical programming approach[END_REF] the updates were synchronous.

The problem

Given a directed graph G = (V, E), a discrete set T of time instants (which we suppose to be an initial contiguous proper subset of N) and the following functions:

• a function α : E → {+1, -1} called the arc sign function; A gene regulatory network (GRN) is a 8-tuple (G, T, α, ω, θ, χ, ι, γ) such that:

∀v ∈ V χ(v, 0) = ι(v) (1) ∀v ∈ V, t ∈ T {0} χ(v, t) =      H(v, t -1) if γ(v, t) = 1 χ(v, t -1) otherwise ( 2
)
where H is the Heaviside function defined for v ∈ V and t ∈ T by 1)-( 2)-(3) together are called the evolution rules of the GRN. For any particular t ∈ T , χ(•, t) : V → {0, 1} is called a configuration. Since the evolution rules relate a configuration at time t with a configuration at time t-1, χ(•, t) is called a fixed configuration (or fixed point) if it remains invariant under the application of one complete cycle of updates encoded by γ. Furthermore, as long as the evolution rules are purely deterministic (as is modelled above), a fixed point of a GRN is determined by its initial configuration.

H(v, t) =        1 if u∈δ -(v) α(u, v)ω(u, v)χ(u, t) ≥ θ(v) 0 otherwise, (3)
with δ -(v) = {u ∈ V | (u, v) ∈ E} for all v ∈ V . Eqns. (
In this paper we deal with an inverse problem related to the estimation of update sequence in GRNs. More precisely, we address the following.

Update Sequence Estimation in GRNs (USEGRN). Given a digraph G, a time instant set T , an arc sign function α, an arc weight function ω, a threshold function θ and a set I of initial configurations, find an update function γ with the property that for all ι ∈ I there exists a gene activation function χ such that (G, T, α, ω, θ, χ, ι, γ) are GRNs whose fixed points are at a minimum distance to observed data.

In other words, we attempt to estimate the sequence of updates in a GRN from the knowledge of the digraph topology in such a way that (a) the GRN evolution rules are consistent with respect to a certain set of initial configurations and (b) the fixed points induced by the estimated values are as close as possible to the observed ones.

As the reader might notice, the problem strongly depends on the modelling of the update sequence encoded by γ. In [START_REF] Gershenson | Classification of random boolean networks[END_REF], the authors proposed to describe such a sequence by means of periods and delays parameters for each gene. Assuming p v and d v to be such values for gene v, we can reformulate Equation 2in the previous modellisation according to the following relation:

∀t ∈ T γ(v, t) = 1 ⇐⇒ ∃n ∈ N s.t. t = np v + d v 3
The mathematical programming formulation

The methodology we shall follow is that of modelling the USEGRN by means of a mathematical programming formulation:

min x f (x) subject to g(x) ≤ 0,     
where x ∈ R n are the decision variables and f : R n → R is the objective function to be minimized subject to a set of constraints g : R n → R m which may also include variable ranges or integrality constraints on the variables. The primary concern in solving the USEGRN is thus modellistic rather than algorithmic. One of the foremost difficulties is that of employing a static modelling paradigm -such as mathematical programming -in order to describe a problem whose very definition depends on time. Another important difficulty resides in describing the necessary and sufficient conditions for a configuration to be a fixed point in a mathematical form. We solve this difficulty by introducing two decision variables: a binary variable s stating that the network has been stable for at least two successive time steps; a binary variable y that will indicate the first time the network is stable. The last difficulty concerns the proper modelling of the update sequence as proposed in [START_REF] Gershenson | Classification of random boolean networks[END_REF]. The solution relies on the use of two binary variables π and δ for each gene and indexed over the possible values for the periods and the delays. Then, π v (p) (resp. δ v (d)) is set to 1 if the period (resp. delay) of v is p (resp. d). We provide above such a formulation:

• Sets: V of genes in the network, E of edges in the network, T of time instants, P of periods values, D of delay values and R of regions. • Parameters:

• ι : R × V → {0, 1} is the initial configuration of the network (vector of boolean values affected to the genes) for each region. • α : A → {+1, -1} is the sign of the arc weights; • w : V → R + is the arc weight function;

• θ : V → R is the threshold function;

• φ : V × R → {0, 1} is the targeted fixed configuration for region r.

• Variables:

• for all r ∈ R, v ∈ V , t ∈ T , x t r,v ∈ {0, 1} is the activation state of gene v at time t in region r;

• for all r ∈ R, v ∈ V , t ∈ T , h t r,v ∈ {0, 1} is the projection of state of gene v at time t in region r according to Heaviside function; • s : R × T → {0, 1} is a decision variable indicating that the network is stable during at least two successive time steps in region r. • y : R × T → {0, 1} is a decision variable that indicates the first time the network reaches a stable state in region r. • for all v ∈ V , p ∈ P , π v,p ∈ {0, 1} is a decision variable that indicates that the periodicity of gene v is p.

• for all v ∈ V , d ∈ D δ v,d ∈ {0,
min r∈R t∈T \{1} (y t-1 r -y t r ) v∈V |x t r,v -φ r,v | .
• Constraints:

• Heaviside function computation rule (for all t ∈ T \ {1}, v ∈ V, r ∈ R) :

θ v h t r,v -|V |(1 -h t r,v ) ≤ u∈δ-(v) α uv w uv x t-1 r,u ≤ (θ v -1)(1 -h t r,v ) + |V |h t r,v • state transition rules (for all r ∈ R, v ∈ V , p ∈ P , d ∈ D): x 0 r,v = ι r,v ∀t ∈ T \ {1} s.t. t = np + d π v,p δ v,d x t r,v = π v,p δ v,d x t-1 r,v ∀t ∈ T \ {1} s.t. t = np + d π v,p δ v,d x t r,v = π v,p δ v,d h t-1 r,v π v,p δ v,d d ≤ p • fixed point conditions (for all r ∈ R, t ∈ T \ {1}): v∈V |x t r,v -x t-1 r,v | ≤ V s t r v∈V |x t r,v -x t-1 r,v | ≥ s t r y t r f t r = 0 u>t s u r = f t r (1 -y t r ) ≤ f t r (|P | + |D|) 2 ≤ τ ∈T y τ r 4

Reformulations and solutions

The above problem is a nonconvex Mixed-Integer Non-Linear Problem that can be reformulated exactly to a Mixed-Integer Linear Problem using the techniques proposed in [START_REF] Liberti | Reformulations in mathematical programming: A computational approach[END_REF]. After standard mathematical manipulations, all the nonlinearities reduce to product terms of binary and/or integer variables, which can be reformulated by adding new auxiliary variables and constraints as follows:

xy terms (x, y : binary) xz terms (x : binary, z : integer)

η ≥ 0 ζ ≥ z L x η ≤ y ζ ≤ z + (|z L | + |z U |)(1 -x) η ≤ x ζ ≤ z U x η ≥ x + y -1 ζ ≥ z -(|z L | + |z U |)(1 -x)
where z L and z U stand for the boundaries of z and η and ζ are the new variables that replace the products in the equations.

We solved to optimality a few real-life instances from the GRN of Arabidopsis thaliana using AMPL [START_REF] Fourer | The AMPL Book[END_REF] to model the problem and CPLEX [START_REF] Ilog | ILOG CPLEX 10.1 User's Manual[END_REF] to solve it. The size of the GRNs involved were such that CPLEX obtained the optimal solution in a matter of minutes.

  • a function ω : E → R + called the arc weight function; • a function χ : V × T → {0, 1} called the gene state function; • a function ι : V → {0, 1} called the initial configuration; • a function θ : V → R called the threshold function; • a function γ : V × T → {0, 1} called the updating function.

  1} is a decision variable that indicates that the delay of gene v is d.• Objective function:

Acknowledgements

This work was supported by EU FP6 FET Project Morphex. We are deeply grateful to L. Liberti (LIX) for his feedback on the mathematical formulation.


Fabien Tarissan
email: tarissan@lix.polytechnique.fr

Matthieu Latapy
email: matthieu.latapy@lip6.fr

Christophe Prieur
email: prieur@liafa.jussieu.fr

Efficient Measurement of Complex Networks Using Link Queries

Complex networks are at the core of an intense research activity. However, in most cases, intricate and costly measurement procedures are needed to explore their structure. In some cases, these measurements rely on link queries: given two nodes, it is possible to test the existence of a link between them. These tests may be costly, and thus minimizing their number while maximizing the number of discovered links is a key issue. This is a challenging task, though, as initially no information is known on the network. This paper studies this problem: we observe that properties classically observed on real-world complex networks give hints for their efficient measurement; we derive simple principles and several measurement strategies based on this, and experimentally evaluate their efficiency on real-world cases. In order to do so, we introduce methods to evaluate the efficiency of strategies. We also explore the bias that different measurement strategies may induce.

I. PRELIMINARIES

Complex networks, modeled as large graphs, are everywhere in science, society, and everyday life: in transportation (airlines, roads); communication (internet, web, file or email exchanges); social life (collaborations, friendship, economical exchanges); life sciences (interactions between proteins or genes, dependencies between species); language analysis (synonymy, co-occurrences of words); etc. As a consequence, much effort is devoted to the analysis and modeling of such networks, leading to key insight on various key topics like epidemy or information spreading, algorithm and protocol design, resilience to failures and attacks, etc.

However, it must be clear that most real-world complex networks are not directly available: collecting information on their structure generally relies on intricate and expensive measurement procedures. Conducting such a measurement often is a challenge in itself, and is an important part of the work needed to study a complex network.

In general, complex network measurements consist in a combination of a few simple measurement primitives. In several cases, this primitive consists in testing the existence of a link, which we call a link query: given two nodes u and v, a measurement operation makes it possible to decide whether there is a link between them or not. This simple test may be expensive (regarding the needed resources or time, or the load it induces on the network, for instance) and so conducting measurements with as few calls to the measurement primitive as possible is a key issue.

For instance, in online social networks like Facebook or Flickr 1 , privacy concerns and reduction of server load often lead to limitations in the queries that one is allowed to perform to explore networks between users. Link queries are however allowed in most cases. Likewise, measurements of real-world social networks often rely on interviews, in which link queries play a central role [START_REF] Wellman | The network is personal: Introduction to a special issue of social networks[END_REF]. In biological networks like protein interactions or gene regulatory networks, link queries also play a key role [START_REF] Bouvel | Combinatorial search on graphs motivated by bioinformatics applicatitons: a brief survey[END_REF], [START_REF] Guelzim | Topological and causal structure of the yeast transcriptional regulatory network[END_REF].

In all these contexts, and others, link queries are very expensive: they have a significant load on server running online social network software and their number is generally bounded; they have a significant cost for interviewers and participants in sociological studies ; or they require costly biological experiments, depending on the case.

In this paper, we formalise this problem as follows: given a graph G = (V, E), we want to define strategies (ordered lists of link queries) which lead to the discovery of as many links of the network as possible. In other words, we want to minimize the number of link queries while maximizing the number of observed links, i.e. the number of positive answers to these tests 2 .

In order to do so, we will rely on simple intuitions derived from statistical properties observed on most real-world complex networks, which we discuss in Section II. We then propose several measurement strategies in Section III based on these principles. We also need a way to compare and evaluate measurement strategies, see Section IV. We finally use this to experimentally evaluate proposed strategies in Section V.

Before entering in the core of this paper, we give the needed formalism and notations, and discuss related work.

A. Formalism and notations

In all the paper, we will consider an undirected 3 graph G = (V, E), with n = |V | nodes and m = |E| links. We suppose that all the nodes are known, and focus on link discovery only. In other words, we know V but know nothing about E (although we will make some statistical assumptions in accordance with classical empirical observations in the field, see Section II).

We will denote by N (v) the set of neighbors of v ∈ V :

N (v) = {u ∈ V, (u, v) ∈ E} and by d(v) its degree: d(v) = |N (v)|.
A measurement consists in a series of link queries, i.e. tests of the existence of link (u, v) for two nodes u and v in V . At a given stage in such a measurement, one has already discovered a set of links, which we will denote by E ′ ⊆ E. The set of extremities of links in E ′ will be denoted by V ′ ⊂ V . Notice that, although we know V , in general V ′ = V . We will also denote by n ′ the number of nodes in V ′ and m ′ the number of discovered links so far:

n ′ = |V ′ | and m ′ = |E ′ |. We also define N ′ (v) = N (v) ∩ V ′ and d ′ (v) = |N ′ (v)| for all v ∈ V ′ . Notice that both V ′ , E ′ , n ′ , m ′ , N ′
and d ′ vary during a measurement; however, the context will make it clear which value we consider.

B. Related work

This work belongs to the fields of complex network metrology, which mostly focused on the specific case of the internet topology until now, see for instance [START_REF] Lakhina | Sampling biases in ip topology measurements[END_REF]- [START_REF] Dall'asta | A statistical approach to the traceroute-like exploration of networks: theory and simulations[END_REF]. This area of research aims mainly at evaluating the relevance of collected complex network samples and properties observed on them, and correcting these observations. Viewing the measurement as the combination of many instance of a simple primitive (link queries, here) which we want to optimize is new, and is an important contribution of this paper.

Another related problem is the one of link prediction: given a network in which new links may appear, one wants to predict which new links will appear in the future based on currently existing ones [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF], [START_REF] Liben-Nowell | The link prediction problem for social networks[END_REF]. In this context, authors use properties of the known network to infer probable future link, which is similar to what we do below in the measurement context.

II. UNDERLYING PRINCIPLES

Our goal is to design measurement strategies based on link queries (test of the existence of a link between two given nodes) which will minimize the number of such queries and maximize the number of discovered links (i.e. the number of positive answers to these tests). In order to do so, we will rely on some simple statistical properties which are observed on most real-world complex networks [START_REF] Watts | Collective dynamics of small-world networks[END_REF].

A. Properties of complex networks

First, we will suppose that G is sparse: its density δ = 2.m n.(n-1) is very small. In other words, the probability that a link exists between two randomly chosen nodes is very small, i.e. a random link query will fail with high probability.

The second key property is the fact that most complex networks have a very heterogeneous degree distribution (often close to a power law). Since the degree of a node is the number of links attached to it, this means that there is a high variability between the number of links of each node (many nodes have very few links, but some have more, and even many more).

Finally, another key property is the local density: although randomly chosen nodes have a very low probability to be linked, two nodes which have a neighbor in common are linked with a much higher probability. This is generally captured by the clustering coefficient or the transitivity ratio [START_REF] Watts | Collective dynamics of small-world networks[END_REF]- [START_REF]Approximating clustering coefficient and transitivity[END_REF], defined by:

cc(G) = v ∆(v) ∨(v) n tr(G) = 3.∆(G) ∨(G)
where, for each v ∈ V , ∆(v) denotes the number of triangles (sets of three nodes with three links) to which v belongs;

∨(v) = d(v).(d(v)-1)
2 denotes the number of pairs of neighbors of v;

∆(G) = v ∆(v); and ∨(G) = v ∨(v).
A classical observation in complex network studies is that both these quantities are high, at least compared to the density. In other words, if one chooses a random pair of links with an extremity in common (transitivity ratio) or a random node and two of its neighbors (clustering coefficient) then the probability that the third possible link exists is high.

B. Consequences on measurements

The properties above, observed on most real-world complex networks, have a strong impact on measurements and will play a key role here.

First, the low density of complex network implies that randomly choosing two nodes and testing the presence of a link between them is very inefficient. Notice however that, when only link queries are possible, one has no choice but to begin with a series of such random measurements. However, it must be clear that exploring a large complex network with such a strategy only is not reasonable.

Instead, the existence of nodes with degree much larger than the average may be useful for efficient measurement. Suppose that we test a random pair (u, v). The probability that it is positive (i.e. the link (u, v) exists) is proportional to the degree of u (resp. v). Therefore, if it exists then one may guess that u (resp. v) has a high degree, and so testing all pairs (u, w) (resp. (v, w)) for any w will probably lead to the discovery of many links. Notice that u and v play a symmetric role in this reasoning. We will call this observation the degree principle.

Likewise, the high local density may be used for efficient measurement: when we know that two nodes u and v have a neighbor w in common then testing pair (u, v) certainly makes sense as this link exists with high probability. We call this the triangle principle.

We may now turn to the definition of measurement strategies based on these principles.

III. MEASUREMENT STRATEGIES

First notice that when one starts a measurement in our framework, no link is known and we have no way to distinguish between vertices. Therefore, there is no choice but to test random pairs of nodes. We call this null strategy random k . As soon as some links are discovered, though, one may try to design more efficient strategies. The triangle principle indicates that, when a ∨ pattern is discovered one may test the missing link in the triangle. This leads to the following strategy.

Strategy 2: ∨-random k with k an integer.

while m ′ < k do Test a random untested pair (u, v) if (u, v) exists then Test all untested pairs (v, w), for any w in N ′ (u)
Test all untested pairs (u, w), for any w in N ′ (v)

Applying directly the degree principle would lead to a strategy in which we test the pairs (u, v) for all v as soon as a random test led to the discovery of a link of u. However, the degree principle becomes stronger if one waits until several links of a node are found. We therefore propose a strategy in which a series of tests (performed according to another strategy) is followed by a use of the degree principle on nodes for which we discovered many links.

Strategy 3: (∨-)Complete Simple -cs k (resp. ∨-cs k ) with k an integer. Apply random k (resp. ∨-random k ) foreach u ∈ V ′ in decreasing order of d ′ (u) do
Test all untested pairs (u, v), for any v ∈ V This strategy may be improved by using the links it discovers for choosing the next link queries to perform. This leads to the following strategy.

Strategy 4: (∨-)Complete -c k (resp. ∨-c k ) with k an integer. Apply random k (resp. ∨-random k ) Let X = V ′ while X is nonempty do Let u in X with d ′ (u) maximal Remove u from X Test all untested pairs (u, v), for any v ∈ V if (u, v) exists and is the first link of v discovered then Add v to X
One may try to use an even stronger version of the degree principle by noticing that the probability of a link between two nodes is even larger if both have a high degree. Therefore, link queries between nodes for which we already discovered many links have an even higher probability of positive outcome. This leads to the following strategy.

Strategy 5: (∨-)Test-Between-Found -tbf k (resp. ∨-tbf k ) with k an integer. Apply random k (resp. ∨-random k ) foreach (u, v) ∈ V ′ × V ′ in decreasing order of d ′ (u) + d ′ (v) do Test (u, v) if it was untested
Finally, one may try to combine the strategies above in order to improve their efficiency. Indeed, some of them use complementary principles which both help in discovering more links with less link queries. One may therefore expect even better results with combinations of them. We will therefore consider the following strategy.

Strategy 6: (∨-)TBF-Complete -tbfc k (resp. ∨-tbfc k ) with k an integer. Apply tbf k (resp. ∨-tbf k ) Apply c 0
It must be clear that many variants and improvements of the strategies above are possible. Probably, completely different strategies may also be defined. Our goal here however is to evaluate the relevance of the degree principle and triangle principle in the design of measurement strategies. We therefore focus on these relatively simple strategies, which we consider as a natural first set of strategies derived from these basic principles.

IV. EVALUATION METHODOLOGY

For any measurement strategy S, let us define m ′ S (q) as the number of links discovered with q link queries with strategy S 4 . It must be clear that our goal, for a given q, is to design a strategy S that maximises m ′ S (q). Conversely, one may want to discover a given number x of links and ask for the strategy S that will minimize the q such that m ′ S (q) = x. However, given two numbers of queries q and r it is possible that a given strategy S discovers more links with q tests than another strategy T , while T discovers more with r tests (we will observe such a situation in Section V-B). As a consequence, it makes no sense to say that S is better than T , nor the converse; this depends on the allowed number of link queries.

Going further, one may notice that if S and T discover the same number of links after a given number q of tests, but if S discovers more links than T for any number r < q of test, then it seems natural to consider that S surpasses T (it discovers the same number of links, but faster). 4 Notice that, in practice, it is in general impossible to reach a situation where we test all pairs of nodes: q = n.(n-1) 2

, or conversely where we discovered all existing links: m ′ S (q) = m.

A simple way to formalise these intuitions is to define the efficiency of a strategy S for a given number of queries q as the (discrete) integral of the function m ′ S from 0 to q: E q (S) = q i=1 m ′ S (i). Notice that the obtained value will depend on the considered graph, and on q. It seems difficult to avoid this, as the efficiency of strategies do indeed depend on the graph under concern, and on the number of allowed link queries. We will therefore always compare strategies ran on the same graph and with the same number of link queries here.

Another weakness of this definition is that it may give any positive value for the efficiency of a strategy, making it hard to evaluate how far from the worst or best solution we are. In order to avoid this we introduce the normalised efficiency: E q (S) = Eq(S)-Eq(min)

Eq(max)-Eq(min) where min and max stand for the worst and best strategies, i.e. the ones with minimal and maximal efficiencies.

Notice that strategies min and max are easy to determine: min consists in testing pairs of nodes with no links between them as long as possible, thus n.(n-1)

2

m times, and then performing the positive tests; conversely max consists in performing first the m positive tests. As a consequence, we can compute easily E q (min) and E q (max) for any q, and thus obtain the normalized efficiency of any strategy.

The notion of normalized efficiency however remains insufficient. Indeed, as we consider sparse graphs, there are only very few positive link queries, and thus one may expect to be much closer to the min strategy than to the max. As a consequence, the efficiency of any strategy will be very low.

A solution to this problem consists in comparing strategies to the random one, denoted by ran, which consists in performing link queries on random untested pairs of nodes. The expected efficiency of this strategy is easy to compute, as the probability of success of a link query is exactly the density δ; we obtain: E q (ran) = q i=1 i • δ = q.(q+1)

2

• δ. Finally, we introduce the relative efficiency, which indicates how a given strategy S performs compared to the random one (and the minimal and maximal ones) after q link queries: R q (S) = Eq(S) Eq(ran) . Notice that the relative efficiency does not give a value between 0 and 1 and therefore does not have the advantage of being relatively independent from the context. However, we cannot normalize it as we would loos the benefit of the comparison to the random strategy. We will therefore use both the normalized efficiency and the relative efficiency to discuss efficiency of strategies below, and keep in mind that in any case the efficiency of a strategy depends on the graph under concern and on the number of link queries allowed. Only the full m ′ S () function can describe the efficiency of strategy S entirely, on a given graph.

V. EXPERIMENTAL EVALUATION

In this section, we present experiments aimed at illustrating the differences between the proposed measurement strategies, and how they may be evaluated. We first present the dataset we used, which is a typical real-world case. We then examine a typical situation and discuss the observations. We deepen this by observing the impact of the initial random period of measurement; and finally we discuss the bias that measurement strategies may induce on observed properties.

A. Dataset

We use here data on an online social network which we consider as a typical example of complex networks studied in the literature. This social network comes from the Flickr site, which provides facilities for publishing online photos, sharing them with others, discuss them, etc. Users may also subscribe to various interest groups and have lists of other users known as their contacts.

Here we used a complete measurement of Flickr conducted in August 2006 [START_REF] Prieur | The stength of weak cooperation: A case study on flickr[END_REF]. We considered the largest of the 72 875 groups observed then 5 , which contained 31 523 members.

We then defined three different networks among these 31 523 users: a posted a comment on a photo from b and b posted a comment on a photo from a. One may also define a symmetric-contact graph in which two users a and b are linked if both a is a contact of b and b is a contact of a. In order to save space, we will not consider it here. Likewise, we do not detail the features of these networks; the key point here is that they are sparse, have heterogeneous degree distributions and high clustering coefficient and transitivity ratio. To this regard, they are similar to most real-world complex networks, and so the principles discussed in II apply.

B. A typical example

Let us first try all our strategies with the same parameter k = 1 000 and on the contact graph. We represent in Figure 1 the number m ′ S (q) of links discovered by each strategy S as a function of the number q of link queries performed, for q between 0 and Q = 4.10 6 . The obtained plot is representative of what is obtained on other graphs.

This plot shows clearly that measurement strategies perform very differently, and that trying to optimize them is relevant. This is a first important result in itself. Moreover, both the degree principle and the triangle principle are useful in doing so: strategies based on each of them perform significantly better than the random strategy. However, the degree principle seems to be much stronger: while the improvement of ∨-random remains quite low, the improvement obtained by the complete strategy is huge. This is probably due to the fact that, although the clustering coefficient and transitivity ratio are much larger than the density, they remain quite small; instead, the largest degrees in the graphs are relatively close to its number of nodes.

The best final results (the largest number of links discovered at the end of the measurement) are obtained with mixed strategies, namely tbf-complete and ∨-tbf-complete, which succeed in discovering between 22 and 23% of existing links by performing only 1% of possible link queries. They slightly outperform complete, which was expected as they are more subtle (though a stronger improvement may have been expected).

Notice that, although these strategies finally outperform complete, they discover links later than this strategy. In this sense, they may therefore be considered as less efficient, which is captured by our notion of efficiency, see Table I

C. Impact of the initial phase

In order to test the impact of the initial phase on the efficiency of the strategies, we conducted a similar experiment in which we increased the parameter k to 1 500. We present the results in Figure 2. The change of the k value implies in particular that the random phase will last longer than in the previous runs as it looks for a larger set of discovered links. This induces a delay before the beginning of the second phase of the strategies which should in turn decrease their efficiency as they have less queries to test the existence of the links. Surprisingly though, the efficiency of the strategies does not seem to be affected. The amount of discovered links after the same number of queries is for instance comparable in the contact network case (around 21% of the existing links). This can be explained by the fact that while searching for the 1 500 links, the random phase has improved the partial knowledge of the network topology. It is very likely then that the highly connected nodes have emerged more significantly during this phase. Thus the ordering used by the elaborated strategies, based on the degree principle, is in turn more pertinent.

The plots based on the comment and symmetric-comment networks also show that the behaviour of the strategies can be very similar in some cases. This suggests to investigate other criteria to sort out their efficiency.

D. Measurement bias

Until now, we focused on our ability to discover many links with as few link queries as possible. However, different strategies discover different links, which may have consequences on the properties of the obtained samples: they may be biased by the measurement strategy, and biased differently depending on the strategy we use. This can be observed visually in Figure 3 for instance, and confirmed by the statistics given in Table II. of the samples obtained by each measurement strategy with k = 1500 applied on the symmetric-comment network. We also display the properties of the actual network (first row), for comparison.

These experiments clearly show that the observed properties are biased by the measurement (they are not the same as the ones of the actual network), and moreover that different strategies lead to different bias.

One can notice for instance that the complete and the test-between-found strategies induce a very different bias on the properties. It is likely to be due to the fact that the number of involved nodes (m ′ ) is very low for the second strategy but all the possible links between them have been tested. This means in particular that all the possible triangles have been discovered which leads naturally to an over-evaluation of the clustering coefficient. The strict opposite happens in the complete case since many nodes are discovered but the links between them are not directly tested.

In some specific cases though, the values are correctly evaluated by the strategies. Strategy ∨-random, for instance, gives a correct value of the transitivity ratio. This is well explained by the strategy itself that tests the existence of the third link of a triangle as soon as two nodes appear to have a common neighbor.

It is also worth noticing that the mixed strategies have a better evaluation of the clustering coefficient than other strategies. This can be explained by the fact that, as the name suggests, they mix the effects of the different strategies. In particular, the over-evaluation of this property given by the test-between-found phase seems to be compensated by the under-evaluation of the complete phase.

These observations suggest to put in perspective the quantitative assessments of the runs and to try integrating the qualitative point of view in the evaluation of the efficiency of the strategies.

VI. CONCLUSION AND PERSPECTIVES

In this paper, we studied the problem of measuring large complex networks when the measurement operation consists in testing the existence of a link between two nodes. We proposed different strategies for ordering the link queries in order to minimize their number while maximizing the number of discovered links. Those strategies rely on the expected statistical properties of the network in order to predict the existence of the links and we tested this approach on several real-world networks based on the Flickr database.

The empirical results confirmed that the principles underlying the development of the strategies are relevant in this measurement context. The experiments showed that the elaborated strategies made a huge improvement compared to the random approach. But they also raised the question of accounting for the bias they induce on the extracted samples. It turned out that the different strategies gave different evaluations of the statistical properties of the original networks. This result suggests to try combining them in order to compensate those unwanted effects, which is what we plan to investigate more specifically in the future.

Going further, we also intend to extend the kind of realworld networks on which test the strategies, add measurement properties to the list of statistical properties considered (such as the assortativity and the degree-degree correlation) and try to adapt the strategies to the directed graphs.

Strategy 1 :

 1 random k with k an integer. while m ′ < k do test a random untested pair

•

  contact: two users a and b are linked if a is a contact of b or b is a contact of a; • comment: two users a and b are linked if a posted a comment on a photo from b or b posted a comment on a photo from a; • symmetric-comment: two users a and b are linked if both

Fig. 1 .

 1 Fig. 1. Number of links (vertical axis) discovered by each strategy as a function of the number of link queries performed (horizontal axis) in a typical case (contact network, k = 1 000). The tbf, tbf-complete and ∨-tbf-complete strategies are indicated as a unique curve (named TBFs) in the plot as the three curves overlap eachother.

Fig. 3 .

 3 Fig. 3.Drawings of samples obtained with the complete (left) and tbf-complete (right) strategies after the 20 000 link queries. The position of the nodes is the same in the two drawings (it is obtained by a classical graph drawing algorithm ran on the actual network), which makes it possible to observe visually that the links discovered by each strategy are not the same.

Fig. 2 .

 2 Fig.2. Number of links (vertical axis) discovered by each strategy as a function of the number of link queries performed (horizontal axis) for each of our three graphs (from left to right: contact, comment and symmetric-comment), with k = 1 500. The tbf, tbf-complete and ∨-tbf-complete strategies are indicated as a unique curve (named TBFs) in the two last plots as the three curves overlap eachother.

  .

		m ′	% tested	% found	E	R
	random	9 609	1.04	1.03	0.006	0.99
	∨-random	21 030	1.04	2.25	0.010	1.64
	c 1000	209 485	1.04	22.4	0.142	24.2
	tbf 1000	68 874	0.46	7.36	0.048	15.6
	tbfc 1000 ∨-tbfc 1000	218 448 214 175	1.04 1.04	23.4 22.9	0.131 0.134	22.3 22.7

TABLE I

 I Efficiency of each strategy after 4.10 6 links queries on the contact network: the number m ′ of discovered links; the percentage of tested pairs of nodes; the percentage of existing links found; and the efficiency coefficients E and R.

http://www.facebook.com/ and http://www.flickr.com/

Notice that, whereas we suppose that link queries are very expensive, the computational cost of each strategy is not our concern here; we consider it as negligible compared to measurement costs, which fits most real-world cases.

This means that we make no difference between (u, v) and (v, u), for any u and v.

FlickrCentral, http://flickr.com/groups/central/