Software:YaDICs

YaDICs

Original author(s)	Coudert Sébastien, Seghir Rian, Witz Jean-françois
Initial release	January 2012 (2012-01)
Stable release	v04.14a / May 27, 2015 (2015-05-27)
Written in	C++
Operating system	Linux
Size	18.4MB
Type	Image processing
License	GPLv2 or later
Website	yadics.univ-lille1.fr/wordpress/

YaDICs is a program written to perform digital image correlation on 2D and 3D tomographic images. The program was designed to be both modular, by its plugin strategy and efficient, by it multithreading strategy. It incorporates different transformations (Global, Elastic, Local), optimizing strategy (Gauss-Newton, Steepest descent), Global and/or local shape functions (Rigid-body motions, homogeneous dilatations, flexural and Brazilian test models)...

Theoretical background

Context

In solid mechanics, digital image correlation is a tool that allows to identify the displacement field to register a reference image (called herein fixed image) to images during an experiment (mobile image). For example, it is possible to observe the face of a specimen with a painted speckle on it in order to determine its displacement fields during a tensile test. Before the appearance of such methods, researchers usually used strain gauges to measure the mechanical state of the material but strain gauges only measure the strain on a point and don't allow to understand material with an heterogeneous behavior. One can obtain a full in plane strain tensor by derivation of the displacement fields. Many methods are based upon the optical flow.

In fluid mechanics a similar method is used, called Particle Image Velocimetry (PIV); the algorithms are similar to those of DIC but it is impossible to ensure that the optical flow is conserved so a vast majority of the software used the normalized cross correlation metric.

In mechanics the displacement or velocity fields are the only concern, registering images is just a side effect. There is another process called image registration using the same algorithms (on monomodal images) but where the goal is to register images and thereby identifying the displacement field is just a side effect.

YaDICs uses the general principle of image registration with a particular attention to the displacement fields basis.

Image registration principle

YaDICs can be explained using the classical image registration framework:

Image registration general scheme

The common idea of image registration and digital image correlation is to find the transformation between a fixed image and a moving one for a given metric using an optimization scheme. While there are many methods to achieve such a goal, Yadics focuses on registering images with the same modality. The idea behind the creation of this software is to be able to process data that comes from a µ-tomograph; i.e.: data cube over 1000³ voxels. With such a size it is not possible to use naive approach usually used in a two-dimensional context. In order to get sufficient performances OpenMP parallelism is used and data are not globally stored in memory. As an extensive description of the different algorithms is given in.

Sampling

Contrary to image registration, Digital Image Correlation targets the transformation, one wants to extracted the most accurate transformation from the two images and not just match the images. Yadics uses the whole image as a sampling grid: it is thus a total sampling.

Interpolator

It is possible to choose between bilinear interpolation and bicubic interpolation for the grey level evaluation at non integer coordinates. The bi-cubic interpolation is the recommended one.

Metrics

Sum of squared differences (SSD)

The SSD is also known as mean squared error. The equation below defines the SSD metric:

[math]\displaystyle{ SSD(\mu,\mathcal{I_F},\mathcal{I_M})=\dfrac{1}{\left | \Omega_F \right |} \sum_{x_i \in \Omega_F} \left( \mathcal{I_F}(x_i) - \mathcal{I_M}({T}_{\mu}(x_i))\right)^{2}, }[/math]

where [math]\displaystyle{ \mathcal{I_F} }[/math] is the fixed image, [math]\displaystyle{ \mathcal{I_M} }[/math] the moving one, [math]\displaystyle{ \Omega_F }[/math] the integration area [math]\displaystyle{ \left | \Omega_F \right | }[/math] the number of pi(vo)xels (cardinal) and [math]\displaystyle{ {T}_{\mu} }[/math] the transformation parametrized by μ

The transformation can be written as:

[math]\displaystyle{ T_{\mu}(x)=x + \left \{\Phi(x)\right \}^{t} \left \{\mu \right \}. }[/math]

This metric is the main one used in the YaDICs as it works well with same modality images. One has to find the minimum of this metric

Normalized cross-correlation

The normalized cross-correlation (NCC) is used when one cannot assure the optical flow conservation; it happens in case of change of lighting or if particles disappear from the scene can occur in particle images velocimetry (PIV).

The NCC is defined by: [math]\displaystyle{ NCC(\mu,\mathcal{I_F},\mathcal{I_M})=\dfrac{\sum_{x_i \in \Omega_F} \left( \mathcal{I_F}(x_i) - \overline{\mathcal{I_F}}\right) \left( \mathcal{I_M}({T}_{\mu}(x_i)) - \overline{\mathcal{I_M}}\right)}{\sqrt{\sum_{x_i \in \Omega_F} \left( \mathcal{I_F}(x_i) - \overline{\mathcal{I_F}}\right)^{2}\sum_{x_i \in \Omega_F} \left(\mathcal{I_M}({T}_{\mu}(x_i))- \overline{\mathcal{I_M}}\right)^{2}}}, }[/math]

where [math]\displaystyle{ \overline{\mathcal{I_F}} }[/math] and [math]\displaystyle{ \overline{\mathcal{I_M}} }[/math] are the mean values of the fixed and mobile images.

This metric is only used to find local translation in Yadics. This metric with translation transform can be solved using cross-correlation methods, which are non iterative and can be accelerated using Fast Fourier Transform .

Classification of transformations

There are three categories of parametrization: elastic, global and local transformation. The elastic transformations respect the partition of unity, there are no holes created or surfaces counted several times. This is commonly used in Image Registration by the use of B-Spline functions and in solid mechanics with finite element basis. The global transformations are defined on the whole picture using rigid body or affine transformation (which is equivalent to homogeneous strain transformation). More complex transformations can be defined such as mechanically based one. These transformations have been used for stress intensity factor identification by and for rod strain by. The local transformation can be considered as the same global transformation defined on several Zone Of Interest (ZOI) of the fixed image.

Global

Several global transforms have been implemented:

• Rigid and homogeneous (Tx,Ty,Rz in 2D; Tx,Ty,Tz,Rx,Ry,Rz,Exx,Eyy,Ezz,Eyz,Exz,Exy in 3D)
• Brazilian (Only in 2D),
• Dynamic Flexion,

Elastic

First-order quadrangular finite elements Q4P1 are used in Yadics.

Local

Every global transform can be used on a local mesh.

Optimization

The YaDICs optimization process follows a gradient descent scheme.

The first step is to compute the gradient of the metric regarding the transform parameters [math]\displaystyle{ \begin{array}{lcl}\dfrac{\partial SSD(\mu,\mathcal{I_F},\mathcal{I_M})}{\partial \mu} &=&\dfrac{2}{\left | \Omega_F \right |} \sum_{x_i \in \Omega_F} \left( \mathcal{I_F}(x_i) - \mathcal{I_M}({T}_{\mu}(x_i))\right)\dfrac{\partial \mathcal{I_M}( {T}_{\mu}(x_i)}{\partial \mu} \\ &=& \dfrac{2}{\left | \Omega_F \right |} \sum_{x_i \in \Omega_F} \left( \mathcal{I_F}(x_i) - \mathcal{I_M}({T}_{\mu}(x_i))\right) \left( \dfrac{\partial {T}_{\mu}(x_i)}{\partial \mu} \right )^{t} \dfrac{\partial \mathcal{I_M}( {T}_{\mu}(x_i))}{\partial x} \\\end{array} }[/math]

Gradient method

Once the metric gradient has been computed, one has to find an optimization strategy

The gradient method principle is explained below:

[math]\displaystyle{ \mu_{k+1}=\mu_k+\alpha_k d_k }[/math]

The gradient step can be constant or updated at every iteration. [math]\displaystyle{ d_k=-\gamma_k \dfrac{\partial \mathcal{C}(\mu,\mathcal{I_F},\mathcal{I_M})}{\partial \mu} }[/math], [math]\displaystyle{ \gamma_k }[/math] allows one to choose between the following methods :

• [math]\displaystyle{ \gamma_k }[/math] [math]\displaystyle{ \Longrightarrow }[/math] steepest descent,
• [math]\displaystyle{ \gamma_k=\left [ \dfrac{\partial \mathcal{C}(\mu,\mathcal{I_F},\mathcal{I_M})}{\partial \mu}\dfrac{\partial \mathcal{C}(\mu,\mathcal{I_F},\mathcal{I_M})}{\partial \mu}^{t} \right ]^{-1} }[/math] [math]\displaystyle{ \Longrightarrow }[/math] Gauss-Newton.

Many different methods exist (e.g. BFGS, conjugate gradient, stochastic gradient) but as steepest gradient and Gauss-Newton are the only ones implemented in Yadics these methods are not discussed here.

The Gauss-Newton method is a very efficient method that needs to solve a [M]{U}={F}. On 1000³ voxels µ-tomographic image the number of degrees of freedom can reach 1e6 (i.e: on a 12×12×12 mesh), dealing with such a problem is more a matter of numerical scientists and required specific development (using libraries like Petsc or MUMPS) so we don't use Gauss-Newton methods to solve such problems. One has developed a specific steepest gradient algorithm with a specific tuning of the αk scalar parameter at each iteration. The Gauss-Newton method can be used in small problems in 2D.

Pyramidal filter

None of these optimization methods can succeed directly if applied at the last scale as the gradient methods are sensitive to the initial guests. In order to find a global optimum one has to evaluate the transformation on a filtered image. The figure below illustrates how to use the pyramidal filter to find the transformation.

Pyramidal process used in Yadics (and ITK).

Regularization

The metrics is often called image energy; people usually add energy that comes from mechanics assumptions as the Laplacian of displacement (a special case of Tikhonov regularization ) or even finite element problems. As one decided not to solve the Gauss-Newton problem for most of cases this solution is far from being CPU efficient. Cachier et al. demonstrated that the problem of minimizing image and mechanical energy can be reformulated in solving the energy image then applying a Gaussian filter at each iteration. We use this strategy in Yadics and we add the median filter as it is massively used in PIV. One notes that the median filter avoids local minima while preserving discontinuities. The filtering process is illustrated in the figure below :

See also

• Image registration
• Optical flow
• Displacement vector
• Particle Image Velocimetry