<figure><img src="image_1.png" /></figure>

HAL
open science

# Curve and boundaries measurement using B-splines and virtual images

Julien Réthoré, M.L.M. François

► To cite this version:

Julien Réthoré, M.L.M. François. Curve and boundaries measurement using B-splines and virtual images. Optics and Lasers in Engineering, 2013, 52, pp.145-155. 10.1016/j.optlaseng.2013.06.018. hal-01006711

HAL Id: hal-01006711

https://hal.science/hal-01006711

Submitted on 7 Jul 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Curve and boundaries measurement using B-splines and virtual images

Julien Réthoré<sup>a,*</sup>, Marc François<sup>b</sup>

Laboratoire de Mécanique des Contacts et des Structures (LaMCoS), Universités de Lyon/INSA Lyon/CNRS
UMR 5259, 20 Avenue des Sciences, F-69621 Villeurbanne Cedex, France

<sup>b</sup>LUNAM Université, Université de Nantes-Ecole Centrale Nantes, GeM, Institute for Research in Civil and Mechanical Engineering, CNRS UMR 6183, 2, rue de la Houssinière BP 92208, 44322 Nantes Cedex 3, France

This paper presents a digital image correlation method able to detect and measure a contour or silhouette in a digital image, with a sub-pixel precision. A virtual image that describes a simple contour or silhouette roughly similar to the physical object to detect is created. The matching with the physical image of the object is then performed based on digital image correlation. We look for the displacement field which transforms the physical image which is the closest to the virtual one. In order to keep a problem of finite dimension, this field is decomposed on a B-spline basis. The algorithm is close to the digital image correlation method; an efficient implementation is proposed that makes use of levelsets. The full method is tested over both synthesized and real images, one of a star shaped toy and one on an ill-defined image of a high temperature metallic specimen on which no other boundary measurement was possible. The influence of the arbitrary choice of the virtual image contour or silhouette arbitrary parameters is discussed with respect to the precision of the method. This one is shown to be sub-pixel in any case.

## 1. Introduction

Detection of contour or silhouette boundary is of interest in various scientific domains, for example, map making, physics, chemistry (detection of fronts), automotive engineering (computer aided driving), production engineering (control), biology (measurement of cell boundaries, filament behavior), reverse mechanical engineering, etc.

Most of the existing methods focus on the capability to detect ill-defined boundaries. Ridge enhancing filters [10,22,23], as well as efficient phase congruency methods [12], emphasize points belonging to a ridge. Numerous approaches of ridge detection are based on the computation of the Hessian of the Gaussian filtered image [7,1,21]. The minimal path methods apply for the detection of ridges in complex images [4,20,5,14]. Finally, the beamlet theory [6] is particularly efficient to identify ridges in noisy images. The precision of such methods varies typically from 1 to 1/10 of the pixel size.

The present method applies for smooth, continuous and non-branched curve detection. This strong hypothesis allows us to obtain a better precision which has already been proven to be close to 1/1000 pixel in [18] and is still relevant for numerous applications (measurement of manufactured objects boundaries,

regular fronts in physics or chemistry, etc.). The Virtual Image Correlation (VIC) method consists in the creation of a virtual image, roughly similar to the contour or silhouette present in the physical image, which is deformed until it matches as best as possible with the contour (or silhouette) present in the "physical" image. The virtual image is based on a curve in the vicinity of which the gray level evolves smoothly with the distance to the curve, according to a bell-shaped function for the contour detection (or to a smooth step for the silhouette detection). In the former version of the VIC [9,18] the curve is described by a function series of the curvature (whose optimal weights are researched). However, due to a computation at any point of the (sub-pixel) grid associated to the curvilinear frame of the curve, this method requires a large amount of CPU time.

In order to improve numerical efficiency, many improvements have been made as follows:

* Levelsets are used to create a curvilinear co-ordinate system within this narrowband of pixels. The first one corresponds to the signed distance of the contour or silhouette and the second one is proportional to the curvilinear abscissa along the curve. The computation of the virtual image gray level, virtual image gradient, and B-splines basis functions is straightforwardly derived for these levelsets.

* The curve (of the virtual image) is obtained from a displacement field applied (in the normal direction) to a simple initial curve (a circle or segment).

<footnote>* Corresponding author. Tel.: +33 472 438 787.
E-mail addresses: julien.rethore@lamcos.insa-lyon.fr, Julien.Rethore@insa-lyon.fr (J. Réthoré), marc.francois@univ-nantes.fr (M. François).</footnote>
* This displacement field is defined by a set of B-spline functions [11] whose weights are the researched parameters. The local effect of the variation of one of these parameters helps to the convergence speed of the method.

This strategy leads to some similarity to the Digital Image Correlation (DIC) technique as developed by [2,17,16] in which an initial image is deformed towards a final one with respect to some deformation field basis. However, major differences between DIC and VIC exist: at first, the virtual image cannot exactly match the physical one thus the optical flow conservation hypothesis cannot be fulfilled; at second, the gradient of the virtual (theoretical) image is perfectly smooth and analytically defined; and, at third, the virtual image is not defined over the full size of the physical one but only exists in the vicinity of the contour or silhouette of interest.

The paper starts by detailing the generation of virtual images. Then the radial displacement field based on B-spline functions is defined. The correlation method that is elaborated is then introduced as well as its implementation based on levelsets. A performance analysis is carried out to estimate the error and uncertainty of the method on synthesised images. Last two examples based on real images are proposed. The first example is aimed at demonstrating the ability of the method to capture tortuous curves with a high geometrical accuracy. It consists in detecting the contour of a star shaped toy made of modeling clay. The second example shows how the proposed technique can be used as an alternative to the full field displacement measurements when this one cannot be used. Indeed, in this high temperature testing, it was almost impossible to apply a speckle texture onto the specimen surface whereas it is usually easy to get images of the silhouette of the specimen. Herein, the contour of the steel specimen under tension loading at 900 °C is detected and its variations due to strain localization and necking are estimated using the proposed methodology.

## 2. The virtual image

Fig. 1 shows an example of a problem of contour detection. The virtual image (Fig. 2) consists in a simple sketch of the contour or silhouette to be identified. It is based on a parameterized curve $x(s)$, where $s$ is the curvilinear abscissa. This curve can be a circle in case of closed contour or silhouette, or a straight segment in case of opened contour or silhouette. With respect to the physical image, this curve must be sized and placed in order to intersect the shape of interest (on the physical image) at least one time. Furthermore, other curves, closer to the shape of interest, can also be used. The virtual image gray level $g(r)$ depends upon the distance $r \in (-R, R)$ from the current point $X$ to the curve

$$X(s, r) = x(s) + r\nu(s) \qquad (1)$$

where $(\tau, \nu)$ denote the local unit tangent and normal vectors of the curve. The gray level function $0 \le g(r) \le 1$ (in this study the gray levels of the physical image are also ranging between 0 and 1) is defined by a sinusoidal bell-shaped function for the contour detection and by a sinusoidal smooth step for a silhouette detection:

$$g^c(r) = \frac{1}{2} \left( 1 + \cos \frac{\pi r}{R} \right) \qquad (2)$$

$$g^s(r) = \frac{1}{2} \left( 1 - \sin \frac{\pi r}{2R} \right) \qquad (3)$$

where suffixes $g^c$ and $g^s$ denote respectively contour and silhouette cases (see Fig. 2). The choice of these functions revealed to be a good compromise between simplicity and efficiency, however,

others have been tested and can be used for example linear variation or more complex ones such as Gaussian functions. It will be shown in Section 7.2 that the best strategy (in terms of precision) consists in beginning with a wide band (a large value of $R$) and reducing it towards, at an ideal, the length of the dark to light gray transition of the edge of the physical image (which can be one pixel in case of perfectly sharp images).

## 3. The displacement field

The virtual image is deformed until it best matches the curve (contour or silhouette) in the physical image. This is done by finding the displacement field $u$ which transforms the virtual image as close as possible to the physical one. It will be seen in Section 5 that the displacement field evolves along the gradient of the virtual image. This imposes the direction of the displacement field to be normal to the curve since the virtual image's gradient is null along its tangent. Let $u(x)$ be the displacement field of the virtual curve, defined at any point $x$ of the curve at the initial state:

$$u(x) = u(s)\nu(s) \qquad (4)$$

where $u(s)$ is the intensity of the displacement, $s$ the curvilinear abscissa and $\nu$ the normal vector to the curve. Fig. 3 shows an example of such a field and corresponding deformed virtual image.

For the correlation process, the displacement field needs to be defined at every point of the virtual image. Since Eq. (4) only defines it at the points which belong to the curve, this field is extended to the neighborhood by considering that a current point $X$ has the same displacement as its projection $x$ (Eq. (1)) onto the curve, i.e. $u(X) = u(x)$.

One can observe in Fig. 3b that such a radial transformation does not keep either the width or the length constant. However, it will be verified in Section 7 that this has a weak influence onto the precision of the measurement.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Example of “physical” image: a picture of a star shaped modeling clay molding. The (closed) silhouette is researched.</figcaption></figure>
### 4. The B-spline basis functions

The function $u(s)$ represents the unknown of the problem whose dimension is thus infinite. For this reason, the field is decomposed over a field basis $P_{i,p}$ which reduces the problem to the determination of the $n$ intensities $\Lambda_i$ thus a finite dimensional problem is as follows:

$$u(s) = \sum_{i=1}^{n} \Lambda_i P_{i,p}(s) \qquad (5)$$

Among the large choice of possible basis (for example Fourier series, polynomials, etc.), the B-spline functions have been retained for their high level of continuity and their local support, which will lead to band matrixes in the correlation calculus which is shown in Section 5. These are classically obtained from the Cox-De Boor recursive formula

$$P_{i,q}(s) = \frac{s-s_i}{s_i-q-s_i} P_{i,q-1}(s) + \frac{s_{i+q+1}-s}{s_{i+q+1}-s_{i+1}} P_{i+1,q-1}(s), \qquad (6)$$

$$P_{i,0}(s) = \begin{cases} 1 & \text{if } s_i \le s < s_{i+1}, \\ 0 & \text{otherwise.} \end{cases} \qquad (7)$$

in which $s_i$ is the knot vector with $i \in \{1, 2, ..., n + p + 1\}$ the knot index and $n$ is the number of basis functions $P_{i,p}$ of order $p$. Furthermore, for the sake of simplicity, the knot vector is chosen uniformly: $s_{i+1} - s_i = s_i - s_{i-1}$. In this context, the maximum continuity $C^{p-1}$ is obtained and an element size which is the step between two consecutive knots can be defined. The elements are defined only over the “active” subset of the segment between the knots $i=p+1$ and $i=n+1$ (those over which $p+1$ B-spline functions can be computed). The number of elements is thus defined as $n_{el} = n - p$. Conversely, the number of B-spline functions of degree $p$ defined over $n_{el}$ elements is $n = n_{el} + p$. Note that the maximum continuity is $C^{p-1}$ when the first and last knots have unit multiplicity. The opposite configuration is when all the internal knot have $p$ multiplicity. In this case $p$-finite element functions ($C^0$) are constructed. The main advantages of using $C^{p-1}$ B-spline functions are that they have high continuity and a low “refinement” cost. Indeed, for $C^{p-1}$ functions, the number of degrees of freedom (the number of B-spline functions) is $n_{el} + p$, whereas it is $n_{el} \times p + 1$ for $C^0$ functions, $n_{el}$ being the number of elements. As a consequence increasing the degree from $p$ to $p+1$ costs one degree of freedom for $C^{p-1}$, whereas it costs $n_{el}$ degrees of freedom for $C^0$. As the optical flow equation used in the following for image correlation leads to an ill-posed problem, the lower the number of degrees of freedom the better the problem is constrained and the less oscillations appear in the solution. Further, the support of B-spline functions enlarges when the level of continuity is increased. Thus, a larger amount of information (gray level gradient) is used to solve the problem for each function.

Finally, joining Eqs. (4) and (5) gives a decomposition of the displacement field at any point $X$ of the domain of definition $\Omega$ of

the virtual image over a basis $u_i$:

$$\mathbf{u}(\mathbf{X}; \Lambda_i) = \sum_{i=1}^{n} \Lambda_i \mathbf{u}_i(s) \qquad (8)$$

$$\mathbf{u}_i(s) = P_{i,p}(s)\nu(s) \qquad (9)$$

### 5. The correlation method

The gray levels are denoted as $f(\mathbf{X})$ for the physical image and $g(\mathbf{X})$ for the virtual one. The best match between the two images is searched by finding the displacement field $\mathbf{u}$ for which the sum of the square difference,

$$\Phi = \iint_{\Omega} [f(\mathbf{X} + \mathbf{u}) - g(\mathbf{X})]^2 d\mathbf{X}, \qquad (10)$$

between them is minimum. In this expression, $\Omega$ denotes the domain of definition of the virtual image, supposed to be fully embedded in the physical one. From Eq. (5), this field and thus the function $\Phi$ depends upon the researched parameters $\Lambda_i$. At this step one can remark that, contrary to the DIC theory, this equation does not accept, in general, zero as a solution since the virtual and physical images are different (in contradiction to the optical flow conservation considered in DIC), however, the problem still consists in minimizing it. Hence, the computation is very similar to the one presented by [2]. From Eq. (8), the equation can be transformed into

$$\Phi(V_k) = \iint_{\Omega} [f(\mathbf{X} + \Lambda_k \mathbf{u}_k) - g(\mathbf{X})]^2 d\mathbf{X}, \qquad (11)$$

in which, and since, the Einstein convention is used. We look for a small variation $\Delta\Lambda_k$ of the field parameters for which this equation becomes stationary. Thus it is linearized around the current parameter value $\Lambda_k$:

$$\Phi(\Lambda_k + \Delta\Lambda_k) = \iint_{\Omega} [f(\mathbf{X} + \Lambda_k \mathbf{u}_k) + \Delta\Lambda_k \nabla f(\mathbf{X} + \Lambda_k \mathbf{u}_k) \cdot \mathbf{u}_j - g(\mathbf{X})]^2 d\mathbf{X} \qquad (12)$$

<figure><img src="image_4.png" /><figcaption>Fig. 3. A displacement field between initial (circular) and final (bean-shaped) shapes. (a) The field, (b) the circular virtual image of Fig. 2a deformed by this field.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Typical virtual images in case of (a) closed contour (b) closed silhouette (c) open contour (d) open silhouette. The virtual image is undefined in hatched region.</figcaption></figure>
where $\nabla$ denotes the spacial derivation. This equation is stationary as soon as the $n$ derivatives $\partial\Phi/\partial\Delta\Lambda_i$ are null. This leads to a linear matrix-vector system of dimension $n \times n$:

$$M_{ij}\Delta\Lambda_j = b_i, \qquad (13)$$

$$M_{ij} = \iint_{\Omega}(\nabla f(\mathbf{X} + \Lambda_k \mathbf{u}_k) \cdot \mathbf{u}_i)(\nabla f(\mathbf{X} + \Lambda_k \mathbf{u}_k) \cdot \mathbf{u}_j) d\mathbf{X}, \qquad (14)$$

$$b_i = \iint_{\Omega}[g(\mathbf{X}) - f(\mathbf{X} + \Lambda_k \mathbf{u}_k)]\nabla f(\mathbf{X} + \Lambda_k \mathbf{u}_k) \cdot \mathbf{u}_i d\mathbf{X}. \qquad (15)$$

To preserve computational cost, $\nabla f(\mathbf{X} + \Lambda_k \mathbf{u}_k)$ is substituted by $\nabla g(\mathbf{X})$ which is calculated once for all since $g$ is the virtual image. Furthermore, being analytically defined, $\nabla g$ is perfectly smooth. The only term to be updated between two iterations is the field $f(\mathbf{X} + \Lambda_k \mathbf{u}_k)$. This computation is done only in the domain $\Omega$, the narrowband in the immediate neighborhood of the contour or silhouette where the gradient of the virtual image is non-zero. In practical computation, this band is only few pixel width, thus the computation time remains very brief and the method is not sensible to the presence of artifacts of any brightness which does not intersect the curve of interest. Furthermore, the B-spline functions $P_{i,p}$ thus by Eq. (9) the fields $\mathbf{u}_i$ are very localized: they overlap (this means they are simultaneously non-zero) only the $p$ previous and following functions. As a consequence, the band width of the matrix $M_{ij}$ is only $2p+1$, where $p$ is the order of the B-spline, leading to fast resolution of the system in Eq. (13). Finally, one can remark that the noise inherent to any imaging process

(thermal noise of the camera sensor for example) leads to rough gradients $\nabla f$. This difficulty is often overcome by a filtering operation on $f$ which induces a loss of information. On the contrary, the present method only deals with the “raw” information contained in $f$.

<figure><img src="image_6.png" /><figcaption>Fig. 6. Synthesized 8 bits, 512 × 512 pixel disc image.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Finite element description of levelsets for a line. (a) Finite element mesh of the narrowband. (b) Levelset $r$ and virtual contour in black. (c) Levelset $s$ and virtual contour.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Finite element description of levelsets for a closed contour. (a) Finite element mesh of the narrowband. (b) Levelset $r$ and virtual contour in black. (c) Levelset $s$ and virtual contour.</figcaption></figure>
## 6. Levelset description

In the previous section, the displacement field has been shown to depend upon the curvilinear abscissa $s$ (Eq. (5)) and that the gray level depends upon the distance $r$. One thus needs a curvilinear coordinate system ($s,r$) within $\Omega$ (the narrowband around the virtual contour). Levelsets [19] provide an efficient way to

compute them. One is associated to $r$, its iso-0 contour being the curve itself (in the virtual image), and the other is associated to $s$. To preserve computational costs, a regular mesh with quadrangle elements is built in the narrowband. The typical size of these elements is set to 10 pixels and the levelset field is then interpolated at the pixel location. However, this allows for a reduction by a factor of about 100 of the numerical cost to

<figure><img src="image_9.png" /><figcaption>Fig. 7. Mean (left) and RMS (right) errors for the detection of the disc image for different element size and interpolation degree $p$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Mean (left) and RMS (right) errors for the detection of the disc with respect to the different curve width $R=R_f=R_g$ and element size.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Mean (left) and RMS (right) errors for the detection of the disc with respect to the different curve width ratios $R_f/R_g$ and element size.</figcaption></figure>
compute the values of the levelsets from the initial virtual shape geometry. The two levelsets are written as

$$s(\mathbf{X}) = \sum_{k \in \mathcal{N}_{\text{band}}} s_k N_k(\mathbf{X}) \qquad (16)$$

$$r(\mathbf{X}) = \sum_{k \in \mathcal{N}_{\text{band}}} r_k N_k(\mathbf{X}), \qquad (17)$$

where $\mathcal{N}_{\text{band}}$ is the set of nodes in region within which the levelsets are computed. This region has to include the narrowband $\Omega$. It is, however, extended slightly from the close vicinity of the curve. $N_k$ are finite element shape functions and $r_k$ and $s_k$ are the corresponding degrees of freedom for $r$ and $s$, respectively.

The nodal absolute values $r_k$ of $r$ are given by the computation of the minimum distance between the node of interest and any point of the curve (which can be discretized at the pixel size or finer, since it consists in an analytical curve). The sign of $r_k$ is the one of the scalar product $\nu(\mathbf{X}-\mathbf{x})$, where in case of silhouette detection, $\nu$ has to be chosen on the right or left hand (depending whether the silhouette is black on white or the contrary).

<figure><img src="image_12.png" /><figcaption>Fig. 10. Average over the element size of the mean error for the detection of a disc as a function of the curve width ratio $R_f/R_g$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. The virtual image corresponding to the analysis of the star image (Fig. 1) at (a) initial and (b) final stages of optimization computed over the support of the levelsets finite element mesh. The smooth transition between black and white is not visible because $R=5$ pixels.</figcaption></figure>

From this levelset $r$, $s$ is computed so that its gradient (normal to the iso-contour) is orthogonal to the gradient of $r$. A minimization of

$$\Psi = \int_{\Omega} \|\nabla s - \nabla r \times \mathbf{e}_z\|^2 \, d\Omega, \qquad (18)$$

is performed, where $\mathbf{e}_z$ is the unit vector normal to the image plane and $\times$ represents the cross product. Using Eqs. (16) and (17), this minimum is attained when the partial derivatives of $\Psi$ with respect to the researched nodal values $s_i$ are null and one obtains

$$(\int_{\Omega} \nabla N_i \cdot \nabla N_j \, d\Omega) s_j = \sum_{k \in \mathcal{N}} r_k \int_{\Omega} (\nabla N_i \times \nabla N_k) \cdot e_z \, d\Omega, \qquad (19)$$

which is a linear matrix-vector system whose straightforward resolution gives the values of $s_j$. In addition to previous orthogonality conditions, a boundary condition must supplement the system. First, a node is chosen arbitrarily to have a 0 value. Then the minimum value of $s$ is subtracted to all the nodal values so that $s$ starts at 0. Fig. 4 shows the result obtained for a sinusoidal line. Fig. 4(a) depicts the mesh that is constructed around the virtual contour. Levelsets $r$ and $s$ are plotted on Figs. 4(b) and (c). For closed contours, a discontinuity in $s=0$ is expected. The finite element decomposition of $s$ is thus enriched using the eXtended Finite Element Method [13]. The discretization of $s$ enriched with discontinuous functions reads as follows:

$$s(\mathbf{X}) = \sum_{k \in \mathcal{N}} s_k N_k(\mathbf{X}) + \sum_{k \in \mathcal{N}_{\text{enr}}} \bar{s}_k N_k(\mathbf{X})(\mathcal{H}(\mathbf{X}) - \mathcal{H}(\mathbf{X}_k)), \qquad (20)$$

where $\mathcal{N}_{\text{enr}}$ is the set of enriched nodes (those that have their support cut by the discontinuity line), $\mathcal{H}$ the Heaviside function (it values 1 on one side of the discontinuity line, 0 on the other side), $\mathbf{X}_k$ the position of node $k$ and $\bar{s}_k$ its corresponding enriched degree of freedom. The same methodology is applied as for the open contour. The only difference is the numerical integration of the enriched basis functions within the elements cut by the discontinuity. One may refer e.g. to [3] or [13]. Fig. 5 shows the results for the levelset description of the contour of a star. One can easily see the discontinuity line in the bottom branch. After the two levelsets have been initialized, the virtual image is computed from the values $r$ interpolated at the pixels from its nodal values $r_k$ using Eq. (2) or (3). The displacement field basis is computed using $s$ evaluated at the pixels from the nodal values $s_k$ and Eq. (8).

<figure><img src="image_14.png" /></figure>
# 7. Performance analysis

In this section, in order to evaluate the metrological performance of the method, a physical image is generated by the same algorithm as used for the virtual image. Thus *f* is noiseless and its dynamic ranges from 0 to 1. However, in order to take into account the digitalization, the gray levels are converted to 8-bit integers. The size of the image is 512 × 512 pixel. The test consists in detecting the circular disk boundary of known theoretical diameter from the image in Fig. 6.

The width of the narrowbands are denoted as $R_f$ for *f* and $R_g$ for *g*. At each curvilinear abscissa *s*, the normal positioning discrepancy $d(s)$ corresponds to the distance between the theoretical circle used for the "physical" image *f* and the deformed curve ($x = u$) of the virtual (recovered) image *g*, measured along the normal $v$. Retained indicators are its mean *D*, for the systematic error, and Root Mean Square (RMS) $\sigma(d)$, for the uncertainty (i.e. the fluctuations of the positioning).

## 7.1. Compared influence of the element size and the degree of interpolation

Using the Finite Element Method or B-splines functions, the question of the element size arises. This one has been varied over one decade. Among other parameters which may influence the quality of the solution, the degree of interpolation *p* is chosen. The latter has been varied from *p*=1 to 5. Other parameters are fixed in this study: narrowband width is set to $R_g = R_f = 10$ pixels. Fig. 7 shows the results of this study.

It is observed that the mean error does not depend either on the element size or on the B-spline degree. The average mean error is about $5.5 \times 10^{-3}$ pixel. Concerning the RMS error, increasing the element size decreases the uncertainty: the lower the number of unknown the better the conditioning of the problem and thus the lower the fluctuations of the solution. The interpolation degree *p* has again no significant influence of the results. Indeed, in $C^{p-1}$ B-spline functions, increasing the interpolation degree from *p* to *p*+1 only costs 1 degree of freedom thus it deteriorates the conditioning of a really small amount. RMS error between $1.0 \times 10^{-3}$ to $4.0 \times 10^{-4}$ pixel is obtained.

## 7.2. Compared influence of the narrowband width and finite element size

In a real picture the narrowband corresponds to the width of the object image boundary. Due to the finite size of the pixels, even if the image is perfectly neat, a pixel in which the boundary passes through has an intermediate gray level. Furthermore, any defocusing or other optical effect (diffraction for example) will generate some shading and thus a wider narrowband. In this section the ratio are varied between "physical" and virtual images narrowband width (respectively $R_g$ and $R_f$) in order to check the influence of the arbitrary choice of $R_g$.

At first, the widths are varied simultaneously, i.e. $R_g=R_f$, from 6 to 30 pixels. The interpolation degree is set to $p=3$. Fig. 8 shows the results. It is observed that the curve width does not modify the behavior of the error when the element size is varying. The mean error is nearly constant and the RMS error decreases. However, it is difficult to conclude on the positive or negative influence of the curve width as the evolution of the mean error and the RMS error seems to have an erratic evolution when $R_g=R_f$ goes from 6 to 30 pixels.

At second, the ratio $R_f/R_g$ is varied. The interpolation degree is set to $p=3$ and the curve width in the physical image is set to $R_f=10$ pixels. The results are shown in Fig. 9. The mean error and the RMS error are plotted as functions of the element size for different curve width ratios $R_f/R_g$. Again the tendency of the error estimates is not changed by the studied parameter. However, mean errors from $7.0 \times 10^{-4}$ to $4.0 \times 10^{-2}$ pixel are obtained.

If one plots the averaged (over the element size) mean error as a function of $R_f/R_g$ relevant conclusions arise. Fig. 10 shows the results. It is clearly concluded that there exists an optimal value for the curve width ratio and obviously this optimal value is 1. The RMS error is lower influenced by the curve width ratio, only a factor of about 10 between the lower and the higher error level for the same element size. Further, it appears that the lower the curve width ratio the lower the RMS error. It is thus concluded that large $R_f/R_g$ ratios lead to higher systematic error and fluctuations whereas lower $R_f/R_g$ ratios give lower uncertainty and systematic error of the same order of magnitude except for the optimal ratio

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Detected contour after 0, 10, 20, 30, 40 and 50 iterations for scale 3 when the initial virtual shape is a circle.</figcaption></figure>
of 1. Indeed, when the curve width ratio is large, i.e. the curve width of the actual image $R_f$ is large compared to the curve width in the virtual image $R_g$, the position of the contour of the virtual shape is not well constrained within the thickness of the physical curve width. Conversely, low curve width ratio, i.e. when the curve width in the virtual image is larger than the curve width in the actual image, the boundary of the shape is constrained within the actual curve width $R_f$ thus leading to lower uncertainty levels. Concerning the systematic error (the mean error) whether the curve width ratio is lower or greater than 1 it always leads to increasing errors (see Fig. 9). The averaged mean error evolution plotted in Fig. 10 also appears to be of the same nature when $R_f/R_g$ is decreasing or increasing from 1. This result is of practical importance as it shows that $R_g$ can first be chosen larger than $R_f$ to obtain an estimation of the shape with low uncertainty levels. Once this first estimation is obtained, the curve width $R_g$ can be decreased until it reaches approximately the value of $R_f$ which is not known in practice. This will allow to decrease the systematic error whereas a low fluctuation level is preserved.

From this a priori performance analysis, the following conclusions arise:

* As for usual Digital Image Correlation analysis, the larger the number of degrees of freedom (or the smaller the window size) the larger is the error measurement both in terms of systematic error and in RMS error (see e.g. [8]).

* However, as mentioned in Section 5 and contrary to the usual Digital Image Correlation method which are searching for a vector displacement field, in the present formulation, only the component of the displacement orthogonal to the shape is tracked. This leads to a reduction of the number of unknown by a factor of 2 and thus an increase of the performance in terms of uncertainty.

* Using $C^{P-1}$ B-splines functions as a shape descriptor increases the performance of the detection algorithm compared to other polynomial functions [15,8]. The local influence of each degree of freedom, but with large support than using e.g. higher order $C^0$ finite element functions, provides for higher robustness with

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Comparison at different places of the contour detected using the segment by segment initialization procedure (blue +) or a circle (red ×). (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)</figcaption></figure>
respect to the local discrepancies. The low “cost” degree increasing allows for improving the detection resolution (finer details can be described but keeping the measurement uncertainty almost constant). Indeed the performance analysis carried out in this section shows that high order $C^{p-1}$ B-spline can be used without increasing the measurement uncertainty.

* A methodology can be elaborated for practical use: the idea is to use several “scales” for which the curve width increases and the descriptor is coarsened. Then the coarse details of the curve, (i.e. large wave length) are first captured and the ability of the technique to describe finer details is progressively restored from coarse “scales” to the finest “scale”.

The advantages of the proposed method for detecting a curve using the coarse to fine strategy above will be used in the next section that is dedicated to actual examples.

## 8. Examples

In this section, the capabilities of the proposed technique are illustrated through two examples with real digital images. The first one is dedicated to the detection of a star shaped toy and the second one to the analysis of a mechanical test at high temperature.

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Detected contours (blue lines) for the six images of the high temperature experiment. On the left side of the specimen, one can see the wires of the thermocouples (seven pairs of thin light gray lines) that were used to measure the temperature of the specimen. Note that these data are not used herein, only the evolution of the specimen contour is tracked out. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)</figcaption></figure>
## 8.1. Contour of a star shaped toy

The modeling clay star image (Fig. 1) in JPEG format (1704 × 2272 pixels) was initially coded with 16-bit RGB encoding which then converted into 8-bit grayscale format. Using the strategy elaborated above, the analysis is performed at three “scales” which correspond to the use of a decreasing width of the virtual image R=20, 10, 5 pixels and, respectively, to the increase of the number of B-spline elements with 25, 50, 100 elements (all of them of degree 3). Finer details of the actual shape are thus progressively resolved. The analysis was performed twice depending on the initial virtual shape: first using a circle of radius 500 pixels and second using the “segment by segment” initialization strategy proposed in [18]. Fig. 11 (a) shows the initial virtual image with this 500 pixels circle. Because R is set to 5 pixels only, the gray level variations within the narrow band are not clearly visible on the picture. Based on the finite element description of the levelsets r and s, the white part of the virtual image appears with a step shaped contour. When the initial virtual shape is a circle the virtual image deforms intensively to reach the actual shape. Fig. 12 shows how the latter analysis converges at scale 3 (R=20 pixels and 25 elements). The computed solution after 0, 10, 20, 30, 40 and 50 iterations (when convergence is obtained) is depicted. The final virtual image is shown in Fig. 11(b). Using the “segment by segment” initiation of [18] gave very similar identification. This can be seen in Fig. 13 on which the contour detected when using the two different initial conditions is plotted. However, it is difficult to quantify the gap between the two solutions as the curvilinear abscissa of the two analysis do not coincide.

## 8.2. Analysis of a high temperature experiment

Performing Digital Image Correlation displacement measurement is not easy for high temperature experiments. First the surface of the sample emits infrared light that must be filtered. If the temperature varies strongly during the experiment different wave lengths may be filtered depending on the actual temperature of the specimen surface. Consequently, it is impossible to set an infrared filter adapted to the entire experiment. Further when very high temperatures are reached (about 1000 °C), it is also difficult to find the appropriate technique to apply a speckle onto the sample surface. However, it is almost always possible to have images of the sample in white on a black background. In this context, the Virtual Image Correlation technique proposed in the paper allows to measure the shape of the sample edge thus giving information on the specimen deformation during mechanical and/or temperature loading.

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Comparison of the detected contours for the high temperature experiment.</figcaption></figure>

Here, a tension test on a cylindrical specimen of 16MND5 steel at 900 °C is analyzed. Fig. 14 shows six images of test over which we have superimposed the contour detected for the specimen. A total of 50 elements of degree 5 are used as a shape descriptor and the analysis is run using three scales. The initial virtual shape is a straight line located on the right of the specimen. A very good agreement is obtained and quantitative estimation of the sample section variation is plotted in Fig. 15. It is observed that for steps 1 (which is the initial geometry of the sample) to 4, the section variation is almost homogeneous whereas for step 5 strain localization starts. For step 6, necking is well captured as a strong reduction of the specimen section is measured. However, one may have to reduce the continuity level at the necking point so that a discontinuity in the profile tangent could be captured. This could have been obtained by inserting multiple knots at the corresponding position in the knot vector describing the curve. Improving the geometrical description of the curve in such cases is one of the perspectives of this work.

Another perspective of this work consists in generating NURBS meshes for iso-geometric analysis [11] directly from the image of an object. Here, the two edges of the unloaded specimen are detected using only one second order B-spline element. Based on the B-splines description of the edges of the solid, a NURBS mesh is then generated over the bulk. The mesh presented in Fig. 16 has been obtained for NURBS function of degree 3 in the vertical direction degree 1 in the horizontal direction. Further, it has been h-refined four times and two times in the vertical and horizontal directions, respectively, compared to the mesh used for the edges' detection.

## 9. Conclusions

In this paper, a curve detection/measurement technique is proposed based on the matching between a real image and a

<figure><img src="image_19.png" /><figcaption>**Fig. 16.** NURBS mesh built from the contour of the specimen.</figcaption></figure>
virtual one. Contrary to the existing methods usually focussing on their ability to detect ill-defined boundaries, the proposed approach can detect curve with sub-pixel accuracy and the proposed algorithm is shown to have uncertainty measurement as low as $10^{-3}$ pixel.

Finite element based levelsets are used for an efficient implementation of the method. The local influence of the shape parameters due to the B-spline functions used as a basis for displacements also acts for fast and robust computation, on the contrary of former version of the VIC method [9,18] in which the series used made each parameter have a global influence. Another advantage of the B-spline basis is their high continuity order and their support in terms of element compared to those of usual finite element shape function which is larger. This leads to better performances in the context of image matching as demonstrated in [15,8].

Using the proposed method, it has been shown in the last example that a NURBS mesh is generated directly from the image of a solid. This opens a wide way for building NURBS meshes for solids that cannot be modeled using CAD software. For example, using IRM images of organs NURBS meshes could be designed directly. Among the possible applications of the method, this is all the more so interesting that its extension to 3D is straightforward.

Another possible development of the method itself consists in the use of the B-spline curves (and not only the use of their basis functions) as a shape descriptor, the control point of them being the researched data. It will apply for the conversion of any existent drawing in B-spline curves, again with sub-pixel precision.

## References

[1] Aylward SR, Bulitt E. Initialisation, noise, singularities, and scale in height ridge traversal for tubular object. IEEE Trans Med Imaging 2002;21(2):61-75.

[2] Besnard G, Hild F, Roux S. 'Finite-element' displacement fields analysis from digital images: application to Portevin-Le Châtelier bands. Exp Mech 2006;46(6):789-803.

[3] Black T, Belytschko T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45:601-20.

[4] Caselles V, Kimmel R, Sapiro G. Geodesic active contours. Int J Comput Vision 1997;22(1):61-79.

[5] Deschamps T, Cohen LD. Fast extraction of minimal paths in 3D images and application to virtual endoscopy. Med Image Anal 2001;5(4):281-99.

[6] Donoho D, Huo X, Beamlets and multiscale image analysis. In Multiscale and multiresolution methods, Lecture notes in computational science and engineering, vol. 20. Berlin: Springer-Verlag; 2002.

[7] Eberly D, Gardner R, Morse B, Pizer S, Scharlach C. Ridges for image analysis. J Math Imaging Vision 1994;4(4):353-73.

[8] Elguedj T, Réthoré J, Buteri A. Isogeometric analysis for strain field measurements. Comput Meth Appl Mech Eng 2011;200(1-4):40-56.

[9] François MLM, Semin B, Auradou H. Identification of the shape of curvilinear beams and fibers. Appl Mech Mater 2010;24-25:359-64.

[10] Haralick R. Ridges and valleys on digital images. Comput Vision Graphics Image Process 1983;22(10):28-38.

[11] Hughes T, Cottrell J, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Meth Appl Mech Eng 2005;194:4135-95.

[12] Kovesi P. Image features from phase congruency. Videre: J Comput Vision Res 1999;1(3):2-27.

[13] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46(1):133-50.

[14] Mueller D, Maeder A. Robust semi-automated path extraction for visualising stenosis of the coronary arteries. Comput Med Imaging Graphics 2008;32(6):463-75.

[15] Réthoré J, Elguedj T, Simon P, Coret M. On the use of NURBS functions for displacement derivatives measurement by digital image correlation. Exp Mech 2009;50(7):1099-116.

[16] Réthoré J, Hild F, Roux S. Extended digital image correlation with crack shape optimization. Int J Numer Meth Eng 2008;73(2):248-72.

[17] Roux S, Hild F. Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches. Int J Fract 2006;140(1-4):141-57.

[18] Semin B, Auradou H, François MLM. Accurate measurement of curvilinear shapes by virtual image correlation. Eur Phys J Appl Phys 2011;56:1-10.

[19] Sethian J. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and material science. Cambridge: Cambridge University Press; 1996.

[20] Sethian JA. Level set methods: evolving interfaces in geometry, fluid mechanics, computer vision and material sciences. Cambridge University Press; 1999.

[21] Steger C. An unbiased detector of curvilinear structures. IEEE Trans Pattern Analysis Mach Intell 1998;20(2):113-25.

[22] Toft PA. Using the generalized Radon transform for detection of curves in noisy images. Proc IEEE ICASSP 1996;4:2219-22.

[23] Zhang Q, Couloigner I. Accurate centerline detection and line width estimation of thick lines using the Radon transform. IEEE Trans Image Process 2007;16(2):310-6.