Identification of the shape of curvilinear beams and fibers

This work concerns the shape identification of curvilinear objects, for example bent beams or wires in mechanics. The beam’s digital picture is analyzed with the introduced Virtual Image Correlation method. This one consists in finding the optimal correlation between the beam’s image and a virtual beam, whose curvature field is described by a truncated series. The gray level and amplitude of the virtual beam does not need to reproduce exactly the ones of the physical beam image. The analytical form of the optimal shape allows one to derive mechanical properties: the identification of the Young’s modulus of a bar is given as an example. We will also show the robustness of the method with regards to the quality of the image.


Introduction
Full-field strain measurements recently took a major place in experimental mechanics [1,2]. However, these methods do not apply to uni-dimensional objects, such as wires or beams images, whose width can be of the order of one pixel. Furthermore such objects may show large deformations, even including the formation of loops.
Methods issued from the image processing field, such as skeletonizing, level set methods, ridge following methods [3] or Radon transform based methods provide numerical results consisting in a set of pixels. This kind of discrete result is unadapted for beam mechanics that requires the computation of the curvature field (the second derivative of the shape is very noisy in this case).
The proposed method consists in finding the best correlation between the physical beam image and a virtual beam. From an approximate initial shape consisting of a set of straight segments, the virtual beam is gradually "deformed" until it perfectly matches (with respect to a quadratic distance) the physical beam image.

Unidimensional sketch of the method
This section explains the method in a simplified 1D case in which x is the abscissa. The goal is to recover the symmetry axis of a physical image f 0 from the digitalized image f . In the present example, f 0 has a symmetry axis at x s = 1.4, ranges from 0 (black) at the border to 1 (white) at its center (Fig. 1a). The available information is the digital image f (Fig. 1b), where the gray level of each pixel is computed as the mean value of f 0 along the pixel (the pixel size corresponds here to an unit). Due to the choice of x s in between two pixels, the digital image, 4 pixel width, has no symmetry axis. The proposed method for the recovery of x s from f consists in finding the abscissa x 0 for which the quadratic distance between f and the virtual image g is minimum. This one, defined for any x ∈ [−R, R] is expressed as: It has a symmetrical bell shape (Fig. 1c, shifted by a trial value x 0 = 5) and its width and shape are a priori different from the one of f 0 . In order to maximize the precision of the method, g is discretized over a mesh (not represented in figure 1) thinner than the pixel size of f . The computation of the quadratic distance (Eq. 2) requires the computation of f * , the (cubic) interpolation of f over the fine mesh of g (Fig. 1d). Generally, f * is not similar to f 0 .
For the sake of simplicity, calculi over f * and g are expressed here as continuous integrals but are computed numerically (over the fine mesh of g). The figure 1e shows the variation of the quadratic distance φ as a function of the shift x 0 : it is minimum at x 0 = 1.396 pixels. Compared to the prescribed x s = 1.4 pixels, this result shows that the proposed method reaches a subpixel precision for this example. Extra tests with various values x s and square or triangle functions f 0 show quite similar accuracy as soon as f 0 is symmetric. Sub-pixel performance is usually observed with such algorithms used in the Digital Image Correlation field; mathematical explanations can be found in [2].
This simplified 1D analysis can be seen as a beam's mean line detection process along some axis x in a 2D image (Fig. 1f); function f corresponding to a cross section of the 2D image.

The virtual beam
We consider an image f of a physical beam (Fig. 1f, 3,5). Compared to the 1D case (previous section), the researched mean line of the beam is no more a position but a curve. This curve represents the mean line of the virtual beam g, defined from its curvatures γ(s), where s is the curvilinear abscissa. They are described by a truncated series: Advances in Experimental Mechanics VII defined, similarly to Eq. (1), from the distance r ∈ [−R, R] from the current point X to the virtual beam mean line then g is a sole function of r: The virtual beam g is not defined out of its definition domain D g (r ∈ [0, R], s ∈ [0, L]).

Optimization process
The searched parameters are the coefficients A k of the series and the two integration constants θ 0 and the ordinate x 0,2 of x 0 (abscissa x 0,1 is fixed because if not the problem is underdefined). These terms are joined together in a pseudo vector V k = {x 0,2 , θ 0 , A k }. The optimization consists in finding the minimum, over V k of the function Φ defined, similarly to Eq. 2, as: that is the quadratic distance between the virtual (g) and physical (f ) beam images, where D g is the definition domain of g and dS the surface element. For the sake of clarity this section uses continuous integrals and does not distinguish f from f * . Writing the condition of minimum ∂Φ/∂V k = 0 (and considering g(X), where X(V k ) is a current point of the virtual beam) gives: The first term, a line integral over the contour ∂D g of D g (n is the outer normal vector) is neglected on the two (small) ends of the beam at s = 0 and s = L. Its value over the lateral borders (r = ±R) can be also neglected if one supposes that f , at these points (in the background of the physical image), is a constant value f | ∂Dg : since g(±R) = 0 (see Eq. 4), using the divergence theorem leads to As div(X) is constant (equals to 2), the surface integral equals 2S where S is the surface of the virtual beam. As S does not change with respect to the shape parameters V k (a kinematic property of the beams of constant length), the right hand term equals to zero and the optimization condition (Eq. 6) reduces to the one proposed by Hild and Roux [2] for the DIC method: The analytical definition of g (Eq. 3, 4) allows a fast computation of grad(g) and ∂X/∂V k (not shown here). The problem is discretized in order to use a Newton iterative process:

Applied Mechanics and Materials Vols. 24-25 361
This, with Eq. (8), gives: which is a matrix equation, that can be rewritten as M kp ∆V p = L k . Its solution ∆V p defines the updated shape of the virtual beam. The iterative process stops with respect to a speed convergence criterion. From a numerical point of view, Eq. (10) is computed with discrete integrals. The mesh of g is obtained from a discretization of s and r; it is more refined than the pixel size (about 3 times in the examples). Computations require the value of f on the mesh of g: as presented in the first section, this is done with a cubic interpolation (the result was referred as f * in the first section). As the shape of the virtual beam changes at each step, this interpolation is reevaluated at each step.
Preliminary identification of the beam shape Fig. (1e) shows that, even in 1D, the quadratic distance is convex only close to the solution x 0 = x s . In this condition, the Newton algorithm requires the initial step to be close to the solution, in other words that D g contains the beam image in f since the first iteration.
For this reason, a preliminary step by step algorithm identifies the beam shape segment by segment. A segment correspond to a small straight virtual beam (of typical length 4R) whose gray level is defined by Eq. 4. Each new segment starts at the middle of the previous. Its angular orientation is defined from the smallest quadratic distance (Eq. 5) between it and the physical image f (the computation of Φ uses again the interpolation of f on the fine mesh of g).  This method has proven to be robust and able to deal with loops. Fig 3 illustrates such preliminary identification. The image is extracted from a film of a thin fiber transported by a fluid flow in a transparent fracture [4]. No preliminary image processing was done. Despite the low definition of the image, the identification is already smooth and accurate. The half-width of the virtual beam was set at R = 1 pixel and the computing time was a few seconds on a current computer.

Technical aspects
Once the preliminary identification is made, the user can choose the type (Fourier or Legendre for instance) and the order N of the series. The initial values of the parameters V k are obtained from this preliminary identification. The process starts and the shape of the virtual beam is gradually deformed until it fits the physical beam shape. definition of the virtual beam, this shape remains smooth (and infinitely differentiable). The precision of the identification depends upon the retained order of the series. Orders too low (here N < 8) do not allow the virtual beam to match the physical shape. The quadratic distance Φ decreases monotonically as N increases and asymptotically reaches a plateau whose value, generally nonzero, depends only upon the similarity between physical and virtual beams (but is not correlated to the precision of the identification).

Validation of the method: identification of the mechanical properties of a cantilever beam
A 2017-T4 aluminium straight bar (diameter 4.95 mm, length 2459 mm) has been fixed in the horizontal chuck of a milling machine in front of a black curtain. It bends under its own weight.
The proposed method provides the shape, angles and the curvatures of the bar. This last one has a strong mechanical meaning in the Timoshenko's beam theory as it is related to the bending moment. Considering that the external actions only reduce to the gravity, the beam theory writes as: in which M is the flexural moment, ρ the linear mass density, g the standard gravity, E the Young's modulus and x 1 denotes the horizontal axis. In the large transformation framework, this problem accepts a numerical solution. This solution is used to validate the precision of our technique. The criterion retained for the comparison is the least square distance between the ordinates x 2 . The best fit was obtained for E = 72 GPa, that corresponds to the value found in the literature. Furthermore, this value was confirmed by a three point bending test that gave E = 72.6 GPa. The mean discrepancy between the ordinates of the analytical shape and the one measured by the present method was only 2 pixels (1.2 mm). Compared to the beam size, this gives a relative error of 5.10 −4 . This identification was done with a Legendre series at an order N = 8. Another tests, in which the chuck was turned of an half turn, shown that this was comparable to the imperfect straightness of the beam. The optical aberrations were not compensated in the present study, taking them into account would increase the precision of the method. Another strong visual information is given by the image of the physical beam in the straightened frame of the virtual beam g (Fig. 6): a correct identification theoretically leads to a symmetrical image, whatever the beam is curved.