Shape determination of unidimensional objects : the virtual image correlation method

The proposed method, named Virtual Image Correlation, allows one to identify an analytical expression of the shape of a curvilinear object from its image. It uses a virtual beam, whose curvature field is expressed as a truncated mathematical series. The virtual beam width only needs to be close to the physical one; its gray level (in the transverse direction) is bell-shaped. The method consists in finding the coefficients of the series for which the correlation between physical and virtual beams is the best. The accuracy and the robustness of the method is shown by the mean of two examples. The first details a Young’s modulus identification from a cantilever beam image. The second is relative to a thermal plume image, that have a weak contrast and a lot of noise.


Introduction
"Unidimensional" objects, designated as beam by the mechanician, are of interest in various scientific fields: biology (human hair, blood vessels...), botany (plant stems...) and of course mechanics (beams in solid mechanics, flexible fibers in fluid mechanics...).These objects present a large aspect ratio (from 100 to 1000 in the examples).
The beam mechanics associates the flexural momentum to the curvature field.The method issued from the image processing field (such as skeletonizing, level set methods, ridge following methods [1] or Radon transform based methods [2]) generally provide a set of pixel as a result: due to its roughness, the second derivative can generally not be used as it.
Due to processors and camera improvements, the full-field measurement methods [3,4] recently took a major role in experimental mechanics.However, they cannot work on unidimensional objects whose width is close to the pixels size.Others aspects prevent to use them: the motion of flexible wires can be many orders larger than their width and the simple shape identification is not possible (because no reference image is available).
The present Virtual Image Correlation (VIC) method uses an image correlation algorithm close to the one used in full-field measurement methods.It consists in finding the best correlation between the physical beam image and a virtual beam.From an approximate initial shape, the virtual beam is gradually "deformed" until it perfectly matches (with respect to a quadratic distance) the physical beam image.

The virtual beam description
The virtual beam is described upon its mean line position x(s) (s is the curvilinear abscissa), its angles θ(s) (with respect to e 1 where (e 1 , e 2 ) is the reference frame) and curvatures γ(s).The curvature field is described by a truncated series: A n γn ( s). ( The tilde terms refer to dimensionless expressions: s ∈ [0, L] = s/L where L is the overall length of the beam.The functions γn are the basis functions of the series description (Lagrange or Fourier in the present case) and the A n are their respective weights.By successive integrations, this equation gives θ(s) and x(s).The virtual beam lightness l(r), where r ∈ [−R, R] and R is the radius of the virtual beam, varies like a shifted cosines function whose width R is chosen slightly larger than the width of the physical beam image.
The physical beam mean line is defined from the best correlation between the physical and virtual beam image.The method is detailed in the next section.The bell shape of l(r) gives, for the most of the case studies, a good determination of the mean line as soon as the physical beam lightness increases from the border to the center, even if its width is close to one pixel.Fig. 1 illustrates the method, in an uni-dimensional point of view (the sketch can be seen as a cross section of the beam).
The function f 0 models a physical lightness of the object: its symmetry axis is set at u 0 = 0.4 (pixels).
The function f 1 represents its discretization over the pixels of the camera.The virtual beam bellshaped function is g(x); it is discretized over a fine grid (here 0.1 pixel).Its shape and width are different from f 0 (x).The function f 2 is the cubic interpolation of f 1 over the fine grid of g.The best correlation ( f 2 (x) − g(x − u)) 2 dx is (numerically) obtained for the shift u = 0.378 pixels.Despite the raw discretization and the different functions f and g, the correlation method provides a good identification of the mean line, at the sub-pixel scale.

The optimization method
The researched parameters are the coefficients A k of the series and the two integration constants θ 0 and x 0 (the values of θ and x for the initial point).This last one is reduced to one of its component, here x 0,2 , in order to fix the curvilinear abscissa of the beam (else 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 as: 14th International Conference on Experimental Mechanics that is the correlation between the virtual (g) and physical ( f ) beam images.Writing the condition of minimum ∂Φ/∂V k = 0 (and considering g(X), where X(V k ) represents the current point of the virtual beam) gives: The first (boundary) term, in which 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 | ∂D g .Since g(±R) = 0 (see Eq. 2), using the divergence theorem leads to: As div(X) is constant (equals 2), the right hand term equals to zero and the optimization condition (Eq.4) reduces to: The analytical definition of g allow a simple calculus (and fast computation) of grad(g) and ∂X/∂V k that is detailed hereafter.Following [4], this equation is discretized in order to use a Newton iterative process: This, with Eq. ( 6), gives: wich is a simple 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.The computation of f −g requires the projection of f on the generally thinner and curvilinear frame of g; this is done by using a cubic interpolation.The term grad(g).∂X/∂Vk is involved in both the expressions of M kp and L k .From Eq. ( 2) and the beam geometry follows: grad(g) = l (r)ν, (9) and, considering that X = x + rν: The second term, collinear to τ, does not need to be computed as it is orthogonal to grad(g).Denoting as θn the primitive of γn the derivatives of ∂x/∂V k are: The displacement of the virtual beam between two steps is close to (∂X/∂V k )∆V k where the fields ∂X/∂V k represent unitary kinematic fields.Figure 2 shows an example of such fields (for clarity the fields are only represented at the mean line location, i.e. ∂x/∂V k ).They play the same role as the unitary displacement fields used in the DIC method [4].
The virtual image G is naturally discretized over the curvilinear frame (s, r) that does not correspond to the square grid of image F. To avoid any loss of information, the mesh size of the virtual image is much smaller.The computation of the second member of Eq. ( 8) requires one to project the luminance field f onto the mesh of G using, here, a cubic interpolation.
The overall length L of the beam is currently not straightforwardly obtained by the VIC method.The end of the fiber is detected from the sharp variation of ∂Φ/∂s which occurs at the point where the virtual beam exceeds the end of the fiber.The previous calculus (precisely the transformation from Eq. 4 to Eq. 6) supposes that the virtual beam already contains the physical beam image.A preliminary shape identification is done.It consists of a discrete version of the Virtual Image Correlation method, in which the virtual beam reduces to a straight segment (of typical length 4R) whose kinematic field is the sole rotation around its first extremity.The segment centers are retained as beam points.This method has proven to be robust and able to deal with loops: at the crossing points, the weight of the pixels belonging to the other branch branch tend to compensate each other.Fig 3 illustrates such preliminary identification.The image is extracted from a film of a thin fiber moved by a fluid flow in a transparent fracture [5].No preliminary image treatment was done; the image is used in inverse video.On this ill-defined image, the radius of the virtual beam was set at R = 1 pixel.The computing time is a few seconds.

Influence of the series order
From the rough identification are computed the corresponding parameters A k for an order N 0 that equals the number of segments (this calculus is not detailed here).The user can choose another one for the precise identification: initial values are obtained either by truncating A k if the chosen order N < N 0 or by filling it by zeros in the opposite case.Fig. 4 shows the result of the identification, for a Fourier series with N = 80: the virtual beam perfectly matches the physical one, even in the detail view.Due to the analytical definition of the virtual beam, this shape remains smooth (and infinitely derivable).The precision of the identification depends upon N: Fig. 5 shows the evolution of the correlation function Φ (Eq. 3) with respect to N (Eq.1).Orders too low (here N < 8) do not allow the virtual beam to

Exemple: identification of the Young modulus of a straight 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, slopes 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 reduces to the gravity, the beam theory writes as: in which M is the flexural moment, ρ = 2700 kg/m 3 the density, E the Young modulus and x 1 denotes the horizontal axis.In the large transformation framework, this problem accepts a numerical solution that is compared to the measurement obtained by our technique.The criterion retained for comparison is the least square distance between the 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 simple three point bending test that gave E = 72.6GPa.Fig. 7 show the discrepancy between analytical and experimental results.The ordinate discrepancy is in average of 2 pixels (1.2 mm).Another tests, in which the chuck was rotated of an half tour, shown that this was comparable to the imperfect straightness of the beam.The optical aberrations, other possible source of errors, were not compensated in that study.The identification was done with a Legendre series with N = 8.The choice of the order is helped by the pictures Φ(N) (such as Fig. 5).Another strong visual information is given by the image of the physical beam in the (unwrapped curvilinear) frame of the virtual one (Fig. 8): a correct identification theoretically leads to a symmetrical image, whatever the curviness of the beam.The figure (9a) represents a thermal plume picture (344x719 pixels): the fluid (sugar sirup) is heated by a circular device on the bottom of the reservoir.This experiment is similar to the thermal plume existing in the earth mantle.The markers are thermochromic liquid crystals which brighten at a given temperature (23.4 • C on this example).The analytical shape of the plume is required for further studies on this Rayleigh-Bénard convection mechanism.The external contour of the plume is represented by a discontinuous collection of pixels, whose luminance vary in a large amount; for example, the figure (9c) shows a domain with very weak luminosity and a discontinuity of the physical pixel line.Various image processing methods have been tried without success on this picture.The raw identification method detailed in Sec 4 proves here its robustness.
The identification has been done with a virtual radius R = 4 pixels and the order of the Fourier series is N = 168; the virtual beam image has (4248x25) pixels.The boundary of the virtual beam is depicted by the dashed line on detail figures (9b and 9b) (the dot spacing represents the virtual line 10004-p.7 discretization).Orders N such as N ≥ 40 give successful identifications; the corresponding correlation distances varies from Φ(N = 40) = 32.1% to Φ(N = 168) = 31.6%.

Conclusions
The proposed method is shown to be accurate and robust with respect to image noise and even loops.It does not require any preliminary image processing.The computation time remains of the order of the minute, within the Matlab environment, a current computer and high-resolution digital image.The method takes into account all the pixels of the physical beam image but does not use the rest of the pixels that constitutes, in general, a large amount of the image.Furthermore, this makes the method naturally insensible to large artifacts in the background (such as stones in Fig. 6).
A future development will consist in the determination of the local flexural stiffness of beam and wires by using two pictures that differ in the gravity direction.This should apply for beams whose free shape is unknown (on the contrary of the example in Sec 6).Other development may concern images sequences, in which the identified shape parameters of an image will be used as a first evaluation for the next image.Three dimensional development is also envisaged, with the use of two cameras.Finally, the method can be easily adapted for the contour identification problem by the use of a step-shaped virtual beam lightness.

Fig. 8 .
Fig. 8.The aluminium bar seen in the unwrapped frame of the virtual beam.