Strain analysis by a total generalized variation regularized optical flow model

ABSTRACT

In this paper, we deal with the problem of estimating the local strain tensor from a sequence of micro-structural images realized during deformation tests of engineering materials. Since the strain tensor is defined via the Jacobian of the displacement field, we propose to compute the displacement field by a variational model which takes care of properties of the Jacobian of the displacement. In particular, we are interested in areas of high strain. The data term of our variational model relies on the brightness invariance property of the image sequence. As prior we choose the second order total generalized variation of the displacement field. This splits the Jacobian into a smooth and a non-smooth part. The latter reflects the material cracks. An additional constraint is incorporated to handle physical properties of the non-smooth part for tensile tests. We prove that the resulting convex model has a minimizer and show how a primal-dual method can be applied to find a minimizer. The corresponding algorithm has the advantage that the strain tensor is directly computed within the iteration process. It is further equipped with a coarse-to-fine strategy to cope with larger displacements. Numerical examples with simulated and experimental data demonstrate the very good performance of our algorithm. In comparison to state-of-the-art engineering software, our method can resolve local phenomena much better.

KEYWORDS:

• Strain analysis
• optical flow
• total generalized variation
• image processing
• variational method
• materials science
• tensile tests
• convex optimization

AMS SUBJECT CLASSIFICATIONS:

• 68U10
• 94A08
• 74R99

1. Introduction

The (Cauchy) strain tensor plays a fundamental role in mechanical engineering for deriving local mechanical properties of materials. In this paper, we are interested in estimating the strain tensor from a sequence of microstructural images of a certain material acquired during in situ deformation tests. An experimental setup of a tensile test is shown in Figure 1.

Figure 1. Illustration of an in situ tensile test. (a) Experimental setup for the tensile test inside a scanning electron microscope. Left: Test material within the microscope. Right: Part of the material with the regions of interest, where a sequence of images (micrographs) is taken for increasing load and (b) Stress–strain curve with three selected micrographs taken under increasing load. The stress $\sigma =F/A_0$ is defined as the applied force F divided by the initial cross section of the specimen $A_0$. The (global) strain is defined as the relative elongation of the specimen measured by the extensometer.

Let $u\text{≔}(u_1,u_2{)}^{\mathrm{T}}$ be the displacement field that describes how each point moves from one image to another. For deformations of a continuum body, the strain tensor ϵ is defined via the Jacobian of the displacement field u by (1) $\begin{array}{rl}\epsilon & =\left(\begin{array}{cc}{\epsilon }_{11}& {\epsilon }_{12}\\ {\epsilon }_{12}& {\epsilon }_{22}\end{array}\right)\text{≔}\frac{1}{2}(\mathrm{\nabla }u+\mathrm{\nabla }{u}^{\mathrm{T}})\\ & =\left(\begin{array}{cc}{\mathrm{\partial }}_x{u}_1& \frac{1}{2}({\mathrm{\partial }}_y{u}_1+{\mathrm{\partial }}_x{u}_2)\\ \frac{1}{2}({\mathrm{\partial }}_y{u}_1+{\mathrm{\partial }}_x{u}_2)& {\mathrm{\partial }}_y{u}_2\end{array}\right).\end{array}$(1)

State-of-the-art software packages for strain analysis like Veddac [1], VIC [2], NCorr [3], and [4–6] apply correlation-based methods to estimate the displacement field and use the result to compute the strain tensor. Roughly speaking, correlation- based methods compare certain windows around each pixel in a predefined search window. Due to the extent of the window around the pixels and since the displacement is often only computed on a coarser grid to reduce the computational effort, the local resolution of these methods is limited. This is especially a drawback when we are interested in the local strain behaviour and the appearance of microstructural damage.

In this paper, we propose to use a variational model to compute the displacement field u from a sequence of images f. Variational models are composed of a data term $E_{\mathrm{Data}}(u;f)$, which incorporates the information given by the image sequence, and a prior $E_{\mathrm{Reg}}\left(u\right)$, which has to be chosen in such a way that its minimizer reflects known or desirable properties of the displacement field: $\underset{u}{\arg \min }E\left(u\right),E\left(u\right)\text{≔}{E}_{\mathrm{Data}}(u;f)+\lambda E_{\mathrm{Reg}}\left(u\right).$

In image processing, the apparent displacement u between image frames is known as optical flow. Variational methods for optical flow estimation go back to Horn and Schunck [7], and there is a vast number of refinements and extensions of their approach. We refer to [8] for a comprehensive overview. In particular, optical flow models with priors containing higher order derivatives of the flow were successfully used, e.g. in [9–17].

Although it seems natural to apply ideas from variational optical flow models also for strain analysis, such methods have rarely been addressed in the literature. The papers [9,11] aim at computing derivatives simultaneously to the optical flow field but are not related to engineering applications. The computation of the (Lagrangian) strain tensor by a variational method was addressed in [18]. There, the authors proposed a smooth fourth order optical flow model which directly computes the strain tensor from an image sequence obtained in a biaxial tensile test with an elastomer. In contrast to our work, they were interested in the macro-scale behaviour and compute the minimizer of their smooth energy function by solving the corresponding Euler-Lagrange equations.

We propose a data term $E_{\mathrm{Data}}(u;f)$ based on the brightness invariance assumption of the image sequence, which is usual in optical flow estimation. The more interesting part is the choice of the prior, where we have in mind that we are interested in local variations of the Jacobian of the displacement field, i.e. of the strain. In the conference paper [19], we built up on the assumption that the displacement field u is additively composed of a non-smooth part v and a smooth part w. Then, we penalized the Jacobian of v and the second order derivative of w by an infimal convolution prior. The infimal convolution (IC) of the first and the second order derivatives was introduced in imaging by Chambolle and Lions [20]. For a recent generalization called ICTV, we refer to [10]. In this paper, we will see that it is better for our task to split the Jacobian of u instead of u itself. To this end, we apply the second order total generalized variation (TGV) of the displacement field as prior. For tensile tests, a physical prior on the non-smooth part of the strain is added. The TGV of functions was introduced by Bredies et al. [21] for image restoration tasks and has meanwhile found many applications. For a discrete variant, see also [22,23]. In particular, TGV regularization was successfully used in connection with optical flow in [15–17]. We mention that there is a software package corresponding to [16] available since 2015. Our first code in [19] was written before any public software package with TGV or infimal convolution prior was available. Our applications show the potential of the method in an impressive way.

The resulting variational model is convex and a minimizer can be computed with primal–dual methods, which are meanwhile standard in image processing. The corresponding algorithm computes the strain tensor directly within the iteration process. This is an advantage compared to standard engineering software, which estimates the displacement field by correlation-based methods first and then computes its Jacobian. Our algorithm is further equipped with a coarse-to-fine strategy to cope with larger displacements as proposed, e.g. in [24–27].

We demonstrate by artificial and real-world examples that the proposed model leads to very good results, which cannot be obtain by state-of-the-art software packages. With our extensive numerical results, we aim at convincing engineers of the advantages of the proposed method.

2. Model

We start with the continuous model and turn to the discrete model for digital images afterwards. The reason is that the discrete notation appears more clumsy and it is easier to get the clue from the continuous notation.

2.1. Continuous model

Suppose that we are given grey-valued images $f_1,f_2:{\mathbb{R}}^2\supset \mathrm{\Omega }\to \mathbb{R}.$ We are interested in the flow field $u=(u_1,u_2):\mathrm{\Omega }\to {\mathbb{R}}^2$, which describes how $f_1$ transforms into $f_2$. Here, we focus on the brightness invariance assumption, which reads with $\mathit{x}\text{≔}(x,y)$ as (2) $f_1\left(\mathit{x}\right)-f_2(\mathit{x}+u(\mathit{x}\left)\right)=f_1\left(\mathit{x}\right)-f_2(x+u_1(\mathit{x}),y+u_2(\mathit{x}\left)\right)\approx 0.$(2) A first order Taylor expansion around an initial optical flow field $\overline{u}=({\overline{u}}_1,{\overline{u}}_2)$ gives (3) $f_2(\mathit{x}+u(\mathit{x}\left)\right)\approx f_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)+⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2\\ {\mathrm{\partial }}_y{f}_2\end{array}\right)(\mathit{x}+\overline{u}(\mathit{x}\left)\right),u\left(\mathit{x}\right)-\overline{u}\left(\mathit{x}\right)⟩.$(3) Plugging (3(3) $f_2(\mathit{x}+u(\mathit{x}\left)\right)\approx f_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)+⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2\\ {\mathrm{\partial }}_y{f}_2\end{array}\right)(\mathit{x}+\overline{u}(\mathit{x}\left)\right),u\left(\mathit{x}\right)-\overline{u}\left(\mathit{x}\right)⟩.$(3) ) into (2(2) $f_1\left(\mathit{x}\right)-f_2(\mathit{x}+u(\mathit{x}\left)\right)=f_1\left(\mathit{x}\right)-f_2(x+u_1(\mathit{x}),y+u_2(\mathit{x}\left)\right)\approx 0.$(2) ), we obtain $0\approx f_1\left(\mathit{x}\right)-f_2(\mathit{x}+\overline{u}\left(\mathit{x}\right))-⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2\\ {\mathrm{\partial }}_y{f}_2\end{array}\right)(\mathit{x}+\overline{u}(\mathit{x}\left)\right),u\left(\mathit{x}\right)-\overline{u}\left(\mathit{x}\right)⟩.$ Later, we apply a coarse-to-fine scheme [24,25] and use the result from one scale as an initialization for the next one. Applying a non-negative increasing function $\varphi :\mathbb{R}\to {\mathbb{R}}_{\ge 0}$ results in the data term (4) ${\int }_{\mathrm{\Omega }}\varphi \left(-⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)\\ {\mathrm{\partial }}_y{f}_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)\end{array}\right),u\left(\mathit{x}\right)-\overline{u}\left(\mathit{x}\right)⟩-f_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)+f_1\left(\mathit{x}\right)\right)\mathrm{d}\mathit{x}.$(4) For the task at hand, we need a prior which recognizes local changes in the displacement field. Therefore, it seems to be useful to consider the smooth global and non-smooth local behaviour of the displacement field. For $z\in {\mathbb{R}}^{d_1,d_2}$ we denote by $|z|$ the square root of the sum of the squared components of z, i.e. its Frobenius norm.

Let us first consider the following Banach spaces of scalar functions $u:\mathrm{\Omega }\to \mathbb{R}$:

• functions of bounded variation $BV\left(\mathrm{\Omega }\right)$ with norm $\parallel u{\parallel }_{\mathrm{BV}}\text{≔}\parallel u{\parallel }_{L_1}+\mathrm{TV}\left(u\right)$, where $\mathrm{TV}\left(u\right)\text{≔}\underset{\varphi \in C_c^{\mathrm{\infty }}(\mathrm{\Omega },{\mathbb{R}}^2),|\varphi |\le 1}{sup}{\int }_{\mathrm{\Omega }}u\mathrm{div}\varphi \mathrm{d}x,$ where $C_c^{\mathrm{\infty }}(\mathrm{\Omega },{\mathbb{R}}^2)$ denotes the space of compactly supported, smooth functions on Ω mapping to ${\mathbb{R}}^2$. The distributional first order derivative Du is a vector-valued Radon measure with total variation $|Du|\left(\mathrm{\Omega }\right)=\mathrm{TV}\left(u\right)$. In particular for $u\in W^{1,1}\left(\mathrm{\Omega }\right)$ the Sobolev space of functions with absolutely integrable weak first order partial derivatives, the regularizer becomes $\mathrm{TV}\left(u\right)={\int }_{\mathrm{\Omega }}\mid \mathrm{\nabla }u\mid \mathrm{d}x,$

• functions of bounded Hessian $BH\left(\mathrm{\Omega }\right)$ with norm $\parallel u{\parallel }_{\mathrm{BH}}\text{≔}\parallel u{\parallel }_{\mathrm{BV}}+{\mathrm{TV}}_2\left(u\right)$, where ${\mathrm{TV}}_2\left(u\right)\text{≔}\underset{\varphi \in C_c^2(\mathrm{\Omega },{\mathbb{R}}^{2,2}),|\varphi |\le 1}{sup}{\int }_{\mathrm{\Omega }}u{\mathrm{div}}^2\varphi \mathrm{d}x$ and ${\mathrm{div}}^2\varphi \text{≔}{\mathrm{\partial }}_{xx}{\varphi }_{11}+{\mathrm{\partial }}_{xy}{\varphi }_{12}+{\mathrm{\partial }}_{yx}{\varphi }_{21}+{\mathrm{\partial }}_{yy}{\varphi }_{22}$. If $u\in W^{2,1}\left(\mathrm{\Omega }\right)$ the term becomes ${\mathrm{TV}}_2\left(u\right)={\int }_{\mathrm{\Omega }}|{\mathrm{\nabla }}^{2}u|\mathrm{d}x,{\mathrm{\nabla }}^{2}u=\left(\begin{array}{cc}{u}_{xx}& u_{xy}\\ u_{yx}& u_{yy}\end{array}\right),$

• of functions of total generalized variation of order two ${\mathrm{BGV}}_{\lambda }^2\left(\mathrm{\Omega }\right)$, $\lambda =({\lambda }_1,{\lambda }_1)\in {\mathbb{R}}_{>0}^2$ [21] with norm $\parallel u{\parallel }_{{\mathrm{BGV}}_{\lambda }^2}\text{≔}\parallel u{\parallel }_{L_1}+{\mathrm{TGV}}_{\lambda }^2\left(u\right),$ where ${\mathrm{TGV}}_{\lambda }^2\left(u\right)\text{≔}\underset{\stackrel{\varphi \in C_c^{\mathrm{\infty }}(\mathrm{\Omega },\mathrm{Sym}({\mathbb{R}}^2\left)\right),|\varphi |\le {\lambda }_2,}{}{\left((\mathrm{div}{\varphi }_1{)}^2+(\mathrm{div}{\varphi }_2{)}^2\right)}^{1/2}\le {\lambda }_1}{sup}{\int }_{\mathrm{\Omega }}u{\mathrm{div}}^2\varphi \mathrm{d}x.$ Here $\mathrm{Sym}\left({\mathbb{R}}^2\right)$ denotes the space of symmetric $2\times 2$ matrices, i.e. we have ${\varphi }_{12}={\varphi }_{21}$. In [28], it was shown that ${\mathrm{TGV}}_{\lambda }$ can alternatively be characterized as $\mathrm{TGV}\left(u\right)=\underset{a\in \mathrm{BD}(\mathrm{\Omega },{\mathbb{R}}^2)}{inf}\{{\lambda }_1|Du-a|(\mathrm{\Omega })+{\lambda }_2|\mathcal{E}a|(\mathrm{\Omega }\left)\right\},$ where $\mathrm{BD}(\mathrm{\Omega },{\mathbb{R}}^2)$ is the space of functions with bounded deformation, and $\mathcal{E}a$ is the weak symmetrized derivative of a. Under sufficient smoothness assumptions, the term reduces to $\mathrm{TGV}\left(u\right)=\underset{a}{inf}\left\{{\lambda }_1{\int }_{\mathrm{\Omega }}|\mathrm{\nabla }u-a|\mathrm{d}x+{\lambda }_2{\int }_{\mathrm{\Omega }}|\tilde{\mathrm{\nabla }}a|\mathrm{d}x\right\},$ where $|\tilde{\mathrm{\nabla }}a|\text{≔}(a_{1,x}^2+\frac{1}{2}(a_{1,y}+a_{2,x}{)}^2+a_{2,y}^2{)}^{1/2}$.

In our task we deal with vector fields $u:\mathrm{\Omega }\to {\mathbb{R}}^2$ and define $\mathrm{TV}\left(u\right)$, ${\mathrm{TV}}_2\left(u\right)$ and $\mathrm{TGV}\left(u\right)$ componentwise.

In the previous paper [19], we assumed that u=v+w can be decomposed additively into a (componentwise) non-smooth part $v\in \mathrm{BV}(\mathrm{\Omega },{\mathbb{R}}^2)$ and a smooth part $w\in \mathrm{BH}(\mathrm{\Omega },{\mathbb{R}}^2)$ and suggested under the above smoothness assumptions the prior (5) $\begin{array}{rl}\mathrm{IC}\left(u\right)& \text{≔}\underset{v+w=u}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }v|+{\lambda }_2|{\mathrm{\nabla }}^{2}w|\mathrm{d}\mathit{x}\right\}\\ & =\underset{w}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }u-\mathrm{\nabla }w|+{\lambda }_2|{\mathrm{\nabla }}^{2}w|\mathrm{d}\mathit{x}\right\},\end{array}$(5) which is the infimal convolution of $|\mathrm{\nabla }\cdot |$ and $|{\mathrm{\nabla }}^2\cdot |$.

In this paper, we decompose $\mathrm{\nabla }u=a+\tilde{a}$ instead of u and use the TGV of the displacement field u as a prior, which reads under the corresponding smoothness assumptions as (6) $\begin{array}{rl}\mathrm{TGV}\left(u\right)& \text{≔}\underset{a+\tilde{a}=\mathrm{\nabla }u}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\tilde{a}|+{\lambda }_2|\tilde{\mathrm{\nabla }}a|\mathrm{d}\mathit{x}\right\}\\ & =\underset{a}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }u-a|+{\lambda }_2|\tilde{\mathrm{\nabla }}a|\mathrm{d}\mathit{x}\right\}.\end{array}$(6) Of course, if a is a conservative vector field, i.e. the Jacobian of some function w, then $\mathrm{IC}$ and $\mathrm{TGV}$ coincide. However, the latter is in general not the case. In our applications, cracks are reflected by the non-smooth part $\tilde{a}$ of $\mathrm{\nabla }u$. Here, it is important to notice the following difference: the IC prior (5(5) $\begin{array}{rl}\mathrm{IC}\left(u\right)& \text{≔}\underset{v+w=u}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }v|+{\lambda }_2|{\mathrm{\nabla }}^{2}w|\mathrm{d}\mathit{x}\right\}\\ & =\underset{w}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }u-\mathrm{\nabla }w|+{\lambda }_2|{\mathrm{\nabla }}^{2}w|\mathrm{d}\mathit{x}\right\},\end{array}$(5) ) enforces sparsity of the Jacobian $\mathrm{\nabla }v$ of the non-smooth part v of u; while the TGV prior enforces sparsity of the non-smooth part $\tilde{a}$ of the Jacobian $\mathrm{\nabla }u$ of u.

2.2. Discrete model

In practice, we are concerned with images $f_1,f_2:\mathcal{G}\to \mathbb{R}$ defined on a two-dimensional rectangular grid $\mathcal{G}\text{≔}\{1,\dots ,N_1\}\times \{1,\dots ,N_2\}$. Besides this notion, we use the alternative representation of an image f as a vector of length $N=|\mathcal{G}|=N_1{N}_2$. In the task at hand, we are given micrographs $f_1$ and $f_2$, e.g. from a tensile test, corresponding to different loads as depicted in Figure 1. Considering only grid points $\mathit{x}=\mathbf{j}\in \mathcal{G}$ and assuming for a moment that also the partial derivatives of $f_2$ are given, the data term in (4(4) ${\int }_{\mathrm{\Omega }}\varphi \left(-⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)\\ {\mathrm{\partial }}_y{f}_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)\end{array}\right),u\left(\mathit{x}\right)-\overline{u}\left(\mathit{x}\right)⟩-f_2(\mathit{x}+\overline{u}(\mathit{x}\left)\right)+f_1\left(\mathit{x}\right)\right)\mathrm{d}\mathit{x}.$(4) ) becomes (7) $\sum _{\mathbf{j}\in \mathcal{G}}\varphi \left(-⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)\\ {\mathrm{\partial }}_y{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)\end{array}\right),u\left(\mathbf{j}\right)-\overline{u}\left(\mathbf{j}\right)⟩-f_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)+f_1\left(\mathbf{j}\right)\right).$(7) This term is only well defined if $\mathbf{j}+\overline{u}\left(\mathbf{j}\right)\in \mathcal{G}$. Here, we use bilinear interpolation to compute the required values of $f_2$ between grid points. If points lie outside the grid, we assume mirror extension of the image at the boundary. By ${\mathrm{\nabla }}_{x}f$ we denote the forward differences in x-direction with Neumann boundary conditions defined as $\left({\mathrm{\nabla }}_{x}f\right)\left(\mathbf{j}\right)\text{≔}\left\{\begin{array}{ll}f(\mathbf{j}+(1,0\left)\right)-f\left(\mathbf{j}\right)& \mathrm{if}\mathbf{j}+(1,0)\in \mathcal{G},\\ 0& \mathrm{otherwise},\end{array}\right.\mathbf{j}\in \mathcal{G}.$ Again, we use the notation ${\mathrm{\nabla }}_{x}f$ both for the function and the vectorized version. Furthermore, we denote also the corresponding difference matrix by ${\mathrm{\nabla }}_x$. The differences in y-direction can be defined analogously. We replace the partial derivatives of $f_2$ in (7(7) $\sum _{\mathbf{j}\in \mathcal{G}}\varphi \left(-⟨\left(\begin{array}{c}{\mathrm{\partial }}_x{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)\\ {\mathrm{\partial }}_y{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)\end{array}\right),u\left(\mathbf{j}\right)-\overline{u}\left(\mathbf{j}\right)⟩-f_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right)+f_1\left(\mathbf{j}\right)\right).$(7) ) by the forward differences ${\mathrm{\nabla }}_x{f}_2,{\mathrm{\nabla }}_y{f}_2$ and apply once more bilinear interpolation to compute values between grid points which results in the data term $\sum _{\mathbf{j}\in \mathcal{G}}\varphi \left(-⟨\left(\begin{array}{c}{\mathrm{\nabla }}_x{f}_2(\mathbf{j}+\overline{u}\left(\mathbf{j}\right))\\ {\mathrm{\nabla }}_y{f}_2(\mathbf{j}+\overline{u}\left(\mathbf{j}\right))\end{array}\right),u\left(\mathbf{j}\right)-\overline{u}\left(\mathbf{j}\right)⟩-f_2(\mathbf{j}+\overline{u}\left(\mathbf{j}\right))+f_1\left(\mathbf{j}\right)\right).$ For our task, we use the function $\varphi \left(t\right)\text{≔}\mid t\mid$ because it is well-known that the ${\ell }_1$-norm reduces the influence of outliers. Then the data term can be rewritten in the convenient matrix–vector notation (8) $E_{\mathrm{Data}}(u;f_1,f_2)\text{≔}\parallel \left(AB\right)u+c{\parallel }_1=\parallel A{u}_1+B{u}_2+c{\parallel }_1,$(8) where $\begin{array}{rl}A& \text{≔}\mathrm{diag}\left(\right({\mathrm{\nabla }}_x{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right){)}_{\mathbf{j}\in \mathcal{G}}),B\text{≔}\mathrm{diag}(\left({\mathrm{\nabla }}_y{f}_2\right(\mathbf{j}+\overline{u}\left(\mathbf{j}\right)\left){)}_{\mathbf{j}\in \mathcal{G}}\right),\\ c& \text{≔}-\mathrm{diag}\left(\right({\mathrm{\nabla }}_x{f}_2(\mathbf{j}+\overline{u}(\mathbf{j}\left)\right){)}_{\mathbf{j}\in \mathcal{G}}){\overline{u}}_1-\mathrm{diag}(\left({\mathrm{\nabla }}_y{f}_2\right(\mathbf{j}+\overline{u}\left(\mathbf{j}\right)\left){)}_{\mathbf{j}\in \mathcal{G}}\right){\overline{u}}_2\\ & +\left(f_2\right(\mathbf{j}+\overline{u}\left(\mathbf{j}\right)){)}_{\mathbf{j}\in \mathcal{G}}-f_1.\end{array}$

For the discrete versions of the TGV prior, we also need the backward differences ${\tilde{\mathrm{\nabla }}}_{x}u$ in x-direction, i.e. $\left({\tilde{\mathrm{\nabla }}}_{x}u\right)\left(\mathbf{j}\right)=\left\{\begin{array}{ll}u\left(\mathbf{j}\right)-u(\mathbf{j}-(1,0\left)\right)& \mathrm{if}\mathbf{j}\pm (1,0)\in \mathcal{G},\\ 0& \mathrm{otherwise},\end{array}\right.\mathbf{j}\in \mathcal{G}$ and analogously for the y-direction. We set $\mathrm{\nabla }\text{≔}\left(\begin{array}{c}{\mathrm{\nabla }}_x\\ {\mathrm{\nabla }}_y\end{array}\right),\tilde{\mathrm{\nabla }}\text{≔}\left(\begin{array}{cc}{\tilde{\mathrm{\nabla }}}_x& \\ \frac{1}{2}{\tilde{\mathrm{\nabla }}}_y& \frac{1}{2}{\tilde{\mathrm{\nabla }}}_x\\ & {\tilde{\mathrm{\nabla }}}_y\end{array}\right),$ and use the tensor products $\mathbf{\nabla }u\text{≔}(I_2\otimes \mathrm{\nabla })u,\stackrel{\mathbf{~}}{\mathbf{\nabla }}u\text{≔}(I_2\otimes \tilde{\mathrm{\nabla }})u.$ Denoting by $\parallel \cdot {\parallel }_{2,1}$ the mixed norm defined for $x\in {\mathbb{R}}^{dN}$ as $\parallel x{\parallel }_{2,1}\text{≔}\sum _{i=1}^N\parallel (x_{i+jN}{)}_{j=0}^{d-1}{\parallel }_2,$ the discrete counterpart of the TGV prior in (6(6) $\begin{array}{rl}\mathrm{TGV}\left(u\right)& \text{≔}\underset{a+\tilde{a}=\mathrm{\nabla }u}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\tilde{a}|+{\lambda }_2|\tilde{\mathrm{\nabla }}a|\mathrm{d}\mathit{x}\right\}\\ & =\underset{a}{inf}\left\{{\int }_{\mathrm{\Omega }}{\lambda }_1|\mathrm{\nabla }u-a|+{\lambda }_2|\tilde{\mathrm{\nabla }}a|\mathrm{d}\mathit{x}\right\}.\end{array}$(6) ) reads $\mathrm{TGV}\left(u\right)\text{≔}\underset{a}{min}\left\{{\lambda }_1\right(\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\left)t\right\}.$

In summary, our model is given by (9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) From a physical point of view, it makes sense to incorporate additional constraints on the displacement field or the strain. In the TGV model, we split the strain into parts representing the global and local features. With this interpretation, it is a reasonable assumption that for tensile tests the strain corresponding to local phenomena is positive in the direction of the applied force. This is motivated by the fact that cracks can only open or widen during a tensile test. Mathematically, if a force is applied in x-direction, we make the restriction ${\mathbf{\nabla }}_x{u}_1-a_1\ge 0,$ so that the model becomes (10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10)

Remark 1

Let us emphasize that strain in materials science is only defined for continuous deformation. In this sense, it is not defined for cracks. In accordance with the definition (1(1) $\begin{array}{rl}\epsilon & =\left(\begin{array}{cc}{\epsilon }_{11}& {\epsilon }_{12}\\ {\epsilon }_{12}& {\epsilon }_{22}\end{array}\right)\text{≔}\frac{1}{2}(\mathrm{\nabla }u+\mathrm{\nabla }{u}^{\mathrm{T}})\\ & =\left(\begin{array}{cc}{\mathrm{\partial }}_x{u}_1& \frac{1}{2}({\mathrm{\partial }}_y{u}_1+{\mathrm{\partial }}_x{u}_2)\\ \frac{1}{2}({\mathrm{\partial }}_y{u}_1+{\mathrm{\partial }}_x{u}_2)& {\mathrm{\partial }}_y{u}_2\end{array}\right).\end{array}$(1) ) in the continuous setting, we refer to its discrete counterpart $\mathrm{\nabla }u$ as the strain. Nevertheless, we will see later that strain in the classical sense is sometimes better resembled by the smooth part a.

By the next proposition, the functions (9(9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) ), resp. (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ) possess minimizers under a mild assumption which was in particular fulfilled for all our practical examples.

Proposition 2

Let $\ker \left(\right(AB\left)\right)\cap \ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }\right)=\left\{\mathbf{0}\right\}$. Then there exist minimizers of (9(9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) ) and (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ).

Proof.

It is sufficient to show that the continuous, convex function E in (9(9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) ) is coercive, which implies the assertion for (9(9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) ). Regarding (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ), the assertion follows immediately since the function remains coercive after adding the constraint.

Assume in contrast that E is not coercive. Then, for some $C\in \mathbb{R}$, there exists a sequence $(u^{\left(r\right)}{)}_{r\in \mathbb{N}}$ with $\parallel u^{\left(r\right)}{\parallel }_2\to \mathrm{\infty }$ and $E\left(u^{\left(r\right)}\right)\le C$. We split $u^{\left(r\right)}=u_a^{\left(r\right)}+u_b^{\left(r\right)}$ with $u_a^{\left(r\right)}\in \ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }\right)$ and $u_b^{\left(r\right)}\in \ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{)}^{\perp }$, where $\ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{)}^{\perp }$ denotes the orthogonal complement of $\ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }\right)$.

Suppose that $\parallel u_b^{\left(r\right)}{\parallel }_2\to \mathrm{\infty }$ as $r\to \mathrm{\infty }$. Then we have $\begin{array}{rl}\mathrm{TGV}\left(u^{\left(r\right)}\right)& =\underset{a}{min}\left\{{\lambda }_1\right(\parallel \mathbf{\nabla }{u}^{\left(r\right)}-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\left)\right\}\\ & =\underset{z}{min}\left\{{\lambda }_1\right(\parallel \mathbf{\nabla }{u}_b^{\left(r\right)}-z{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}z+\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{u}_a^{\left(r\right)}{\parallel }_{2,1}\left)\right\}\\ & =\underset{z}{min}\left\{{\lambda }_1\right(\parallel \mathbf{\nabla }{u}_b^{\left(r\right)}-z{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}z{\parallel }_{2,1}\left)\right\},\end{array}$ where $z\text{≔}a-\mathbf{\nabla }{u}_a^{\left(r\right)}$. Let $\mathbf{\nabla }{u}_b^{\left(r\right)}=w^{\left(r\right)}+{\tilde{w}}^{\left(r\right)}$, where $w^{\left(r\right)}\in \ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}{)}^{\perp }$ and ${\tilde{w}}^{\left(r\right)}\in \ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\right)$. Then $\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{u}_b^{\left(r\right)}=\stackrel{\mathbf{~}}{\mathbf{\nabla }}{w}^{\left(r\right)}$ and since $\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }$ is injective on $\ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{)}^{\perp }$ we have that $\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }{u}_b^{\left(r\right)}{\parallel }_2=\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}{w}^{\left(r\right)}{\parallel }_2\to \mathrm{\infty }$ as $r\to \mathrm{\infty }$. Now $\stackrel{\mathbf{~}}{\mathbf{\nabla }}$ is injective on $\ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}{)}^{\perp }$ such that $\parallel w^{\left(r\right)}{\parallel }_2\to \mathrm{\infty }$ as $r\to \mathrm{\infty }$. Let $z^{\left(r\right)}=z_1^{\left(r\right)}+z_2^{\left(r\right)}$ be a minimizer of $\mathrm{TGV}\left(u^{\left(r\right)}\right)$, where $z_1^{\left(r\right)}\in \ker (\stackrel{\mathbf{~}}{\mathbf{\nabla }}{)}^{\perp }$ and $z_2^{\left(r\right)}\in \ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\right)$. Then, by the equivalence of norms, $\parallel w^{\left(r\right)}-z_1+{\tilde{w}}^{\left(r\right)}-z_2{\parallel }_{2,1}$ can only be bounded as $r\to \mathrm{\infty }$ if $\parallel w^{\left(r\right)}-z_1^{\left(r\right)}+{\tilde{w}}^{\left(r\right)}-z_2^{\left(r\right)}{\parallel }_2^2=\parallel w^{\left(r\right)}-z_1^{\left(r\right)}{\parallel }_2^2+\parallel {\tilde{w}}^{\left(r\right)}-z_2^{\left(r\right)}{\parallel }_2^2$ is bounded, which is only possible if $\parallel z_1^{\left(r\right)}{\parallel }_2\to \mathrm{\infty }$. But then we get $\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}{z}_1^{\left(r\right)}{\parallel }_{2,1}\to \mathrm{\infty }$ as $r\to \mathrm{\infty }$. Thus, we obtain the contradiction $\mathrm{TGV}\left(u^{\left(r\right)}\right)\to \mathrm{\infty }\mathrm{as}r\to \mathrm{\infty }.$

Assume that $\parallel u_b^{\left(r\right)}{\parallel }_2\nrightarrow \mathrm{\infty }$ as $r\to \mathrm{\infty }$. Then there exists a subsequence $(u_b^{\left(r_j\right)}{)}_{j\in \mathbb{N}}$ with $\parallel u_b^{\left(r_j\right)}{\parallel }_2\le \tilde{C}$ for some $\tilde{C}\in \mathbb{R}$ and it follows $\begin{array}{rl}\parallel \left(AB\right)u^{\left(r_j\right)}+c{\parallel }_1& \ge \parallel \left(AB\right)u_a^{\left(r_j\right)}+c{\parallel }_1-\parallel \left(AB\right)u_b^{\left(r_j\right)}{\parallel }_1\\ & \ge \parallel \left(AB\right)u_a^{\left(r_j\right)}+c{\parallel }_1-\tilde{C}.\end{array}$ Since by assumption $\left(AB\right)$ is injective on $\ker \left(\stackrel{\mathbf{~}}{\mathbf{\nabla }}\mathbf{\nabla }\right)$, the right-hand side goes to $+\mathrm{\infty }$ as $j\to \mathrm{\infty }$, which is a contradiction.

3. Algorithm

In this section, we use the primal–dual hybrid gradient method with modified dual variable (PDHGMp) [29,30] to compute a minimizer of (9(9) $E_{\mathrm{TGV}}\left(u\right)=\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}.$(9) ). The minimization of (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ) follows similarly.

We rewrite the minimization problem as (11) $\begin{array}{rl}& \underset{u,a,s,t}{min}\{\parallel (AB)u+c{\parallel }_1+{\lambda }_1\parallel s{\parallel }_{2,1}+{\lambda }_2\parallel t{\parallel }_{2,1}\}\\ & \mathrm{such}\mathrm{that}\left(\begin{array}{cc}\mathbf{\nabla }& -I_{4N}\\ 0& \stackrel{\mathbf{~}}{\mathbf{\nabla }}\end{array}\right)\left(\begin{array}{c}u\\ a\end{array}\right)=\left(\begin{array}{c}s\\ t\end{array}\right).\end{array}$(11) The basic PDHGMp algorithm for this problem is given in Algorithm 1. For the special form, we refer to [31, Alg. 8].

We have to comment on the proximal steps within the algorithm. The update for a is straightforward. Regarding s,t, observe that the problems can be separated into N subproblems, e.g. for s they have the form $\hat{s}=\underset{s\in {\mathbb{R}}^4}{\arg \min }\left\{\parallel s{\parallel }_2+\frac{{\tau }_2}{2{\lambda }_1}\parallel s-x{\parallel }_2^2\right\},$ where $x=(b_1^{\left(r\right)}+\mathbf{\nabla }{u}^{(r+1)}-a^{(r+1)}{)}_{i+jN,j=0,\dots ,3}\in {\mathbb{R}}^4$, $i\in 1,\dots ,N$. The solution of these problems is given by a grouped or coupled shrinkage, see the appendix. If we extend the model by the positivity constraint, i.e. (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ), the update step for s becomes a bit more involved since we need to compute N problems of the form $\hat{s}=\underset{s\in {\mathbb{R}}^4}{\arg \min }\left\{\parallel s{\parallel }_2+{\iota }_{\ge 0}\left(s_1\right)+\frac{{\tau }_2}{2{\lambda }_1}\parallel s-x{\parallel }_2^2\right\}$ with x as above. For $x_1\ge 0$, the term ${\iota }_{\ge 0}\left(s_1\right)$ can be neglected and we end up with the usual coupled shrinkage. For $x_1<0$, we have ${\hat{s}}_1=0$ and for the remaining three components the usual coupled shrinkage can be applied.

Due to the diagonal structure of A and B, the proximal step to get u can be separated into N subproblems of the form (12) $\hat{y}=\underset{y\in {\mathbb{R}}^2}{\arg \min }\{|\alpha y_1+\beta y_2+\gamma |+\frac{1}{2}\parallel y-x{\parallel }_2^2\},$(12) where $\alpha ,\beta ,\gamma \in \mathbb{R}$. The solution is a generalized soft shrinkage of x explained in the appendix.

Remark 3

Computation of the strain ϵ within the primal–dual algorithm

The primal–dual method uses the Lagrangian of (11(11) $\begin{array}{rl}& \underset{u,a,s,t}{min}\{\parallel (AB)u+c{\parallel }_1+{\lambda }_1\parallel s{\parallel }_{2,1}+{\lambda }_2\parallel t{\parallel }_{2,1}\}\\ & \mathrm{such}\mathrm{that}\left(\begin{array}{cc}\mathbf{\nabla }& -I_{4N}\\ 0& \stackrel{\mathbf{~}}{\mathbf{\nabla }}\end{array}\right)\left(\begin{array}{c}u\\ a\end{array}\right)=\left(\begin{array}{c}s\\ t\end{array}\right).\end{array}$(11) ), which contains the summand $⟨(\mathbf{\nabla }u-a)-s,b⟩$, with the dual variable b. Hence, the algorithm computes the desired strain tensor directly within the iteration process and no subsequent computation of the derivative of the optical flow is required.

We updated the basic PDHGMp algorithm by a coarse-to-fine scheme as proposed, e.g. in [24–27] to cope with large displacements. Moreover, we applied the common trick median filtering on every scale to improve the results, see also [26,27]. For more information and the pseudo code for the coarse-to-fine scheme we refer to [32]. To ensure convergence of the algorithm, we use 3000 iterations of the PDHGMp algorithm in each level of the coarse-to-fine scheme. The high number of iterations leads to a computation time of about 20 min which is not competitive to state-of-the-art optical flow methods. We want to emphasize that there is still some potential to optimize the computation time but for our applications this was not necessary since the computation time is negligible in comparison to the time demand of the experiments.

4. Numerical examples

In this section, we present and analyse numerical results of our approaches.

First, we consider artificial data to demonstrate the behaviour of different optical flow models and to underline why the TGV regularized model appears to be appropriate for the described practical tasks. Then, we deal with displacements in aluminium matrix composites (AMC) during tensile tests, where we focus on the detection of local damage and crack propagation. Finally, we demonstrate the flexibility of our method by showing results for different experimental settings and materials.

The algorithm is implemented in C++. Unless stated otherwise, we set (13) ${\lambda }_1\text{≔}0.2\mathrm{and}{\lambda }_2\text{≔}10$(13) for all real-world examples. The positivity constraint (10(10) $\begin{array}{rl}{\tilde{E}}_{\mathrm{TGV}}\left(u\right)& =\parallel \left(AB\right)u+c{\parallel }_1+\underset{a}{min}\{{\lambda }_1\parallel \mathbf{\nabla }u-a{\parallel }_{2,1}+{\lambda }_2\parallel \stackrel{\mathbf{~}}{\mathbf{\nabla }}a{\parallel }_{2,1}\}\\ & \text{subject to}{\mathrm{\nabla }}_x{u}_1-a_1\ge 0.\end{array}$(10) ) on the strain is used for the real-world examples of tensile tests with AMCs. For the artificial examples and the compression and fatigue tests, we do not involve the constraints.

The data in Figure 9 was provided by K. Lichtenberg from the ‘Institute for Applied Materials (IAM)’ at the Karlsruhe Institute of Technology (KIT).

4.1. Artificial examples

To explain the differences between various regularization terms, we use the following simple example: a segment of $100\times 100$ pixels of one exemplary micrograph of an aluminium composite is taken and the simulated displacement field in Figure 2 consisting of the sum of a piecewise constant and a linear part is applied to warp this image. Then, the warped and the initial image are used as the input to reconstruct the simulated displacement field and its strain tensor via variational models with the same data term (8(8) $E_{\mathrm{Data}}(u;f_1,f_2)\text{≔}\parallel \left(AB\right)u+c{\parallel }_1=\parallel A{u}_1+B{u}_2+c{\parallel }_1,$(8) ) and the following priors: besides the TGV prior, we used the

1. ${\mathrm{H}}^1$ regularizer, e.g. used in the Horn–Schunck model [7], (H1) $\lambda {\mathrm{H}}^1\left(u\right)\text{≔}\lambda \parallel \mathbf{\nabla }u{\parallel }_2^2,$(H1)

2. $\mathrm{TV}$ regularizer proposed in [33] and used for optical flow, e.g. in [25], (TV) $\lambda \mathrm{TV}\left(u\right)\text{≔}\lambda \parallel \mathbf{\nabla }u{\parallel }_{2,1},$(TV)

3. ${\mathrm{TV}}_2$ regularizer with second order differences (TV2) $\lambda {\mathrm{TV}}_2\left(u\right)\text{≔}\lambda \parallel {\mathbf{\nabla }}^{2}u{\parallel }_{2,1},$(TV2)

4. $\mathrm{TV}\text{-}{\mathrm{TV}}_2$ regularizer, e.g. in [34], (TV-TV2) ${\mathrm{TV}}_{1\text{-}2}\left(u\right)\text{≔}{\lambda }_1\parallel \mathbf{\nabla }u{\parallel }_{2,1}+{\lambda }_2\parallel {\mathbf{\nabla }}^{2}u{\parallel }_{2,1},\text{-}{\mathrm{TV}}_2$(TV-TV2)

5. $\mathrm{IC}$ regularizer (IC) $\mathrm{IC}\left(u\right)\text{≔}\underset{u=v+w}{min}{\lambda }_1\parallel \mathbf{\nabla }v{\parallel }_{2,1}+{\lambda }_2\parallel {\mathbf{\nabla }}^{2}w{\parallel }_{2,1}+{\lambda }_3\parallel w{\parallel }_2^2,$(IC) where we added a small quadratic term to avoid the ambiguity introduced by the difference matrices. We set ${\lambda }_3\text{≔}5\times {10}^{-5}$ . Note that such a term in relation with the difference operators was not necessary for the $\mathrm{TGV}$ regularizer since we do not decompose u, but the strain.

Figure 2. Original displacement $u_1$ in pixels, additively composed of a piecewise constant and a linear part and ${\mathrm{\nabla }}_x{u}_1$.

The minimizers $u_1$ of the various models and the strain component ${\epsilon }_{11}={\mathrm{\nabla }}_x{u}_1$ are shown in Figure 3(a,b), respectively. The parameters stated in the caption are optimized with respect to the best visual impression. The models with ${\mathrm{H}}^1$ and ${\mathrm{TV}}_2$ priors cannot find a sharp displacement field boundary, while the $\mathrm{TV}$ model introduces additional boundaries due to the so-called staircasing effect. The $\mathrm{TV}\text{-}{\mathrm{TV}}_2$ regularized model is better, but it is clearly worse than the IC and TGV methods. Clearly, for this example both methods IC and TGV give good results since $u_1$ is indeed additively composed of a non-smooth and a smooth component. The next example does not have this property.

Figure 3. Results for the simulated example from Figure 2 for various regularization terms. Parameters: $\lambda =50$ for (${\mathrm{H}}^1$), $\lambda =0.1$ for (TV), $\lambda =0.1$ for (${\mathrm{TV}}_2$), ${\lambda }_1=0.1$, ${\lambda }_2=0.02$ for ($\mathrm{TV}\text{-}{\mathrm{TV}}_2$), ${\lambda }_1=0.1$, ${\lambda }_2=1$, ${\lambda }_3=0.5\times {10}^{-5}$ for ($\mathrm{IC}$) and ${\lambda }_1=0.1$, ${\lambda }_2=2$ for (TGV). (a) Displacement $u_1$ in pixels using various regularization terms and (b) Derivative ${\mathrm{\nabla }}_x{u}_1$ using various regularization terms.

Our next artificial example in Figure 4 shows a displacement field consisting of two parts. The lower part contains a jump whereas we have a purely linear transition in the upper part. In terms of a tensile test, this might be seen as a crack in the lower part and a purely elastic deformation in the upper part of the image. As the displacement field is a transition between the non-smooth upper part and the smooth lower one, it can not be split additively in an appropriate way, whereas it is actually possible to split the strain. Therefore, the TGV model clearly outperforms the IC model. The parameters for the TGV model are those in (13(13) ${\lambda }_1\text{≔}0.2\mathrm{and}{\lambda }_2\text{≔}10$(13) ). Compared to the TGV model, we chose smaller parameters for the IC model, namely ${\lambda }_1=0.1$, ${\lambda }_2=1$, ${\lambda }_3=0.5\times {10}^{-5}$ , since otherwise the resulting displacement was too smooth. We also added the result of the large displacement optical flow (LDOF) model by Brox and Malik [35] using their implementation with standard parameters. However, it is clearly visible that the result contains too much structures in the smooth upper part and at the same time is not able to reconstruct the jump in the lower part.

Figure 4. Results for a simulated example with a linear transition of the displacement field in the upper half and a jump in the lower half of the image. (a) Displacement $u_1$ in pixels using IC, TGV regularization, and the method in [35] and (b) Derivative ${\mathrm{\nabla }}_x{u}_1$ using IC, TGV regularization, and the method in [35].

4.2. Tensile tests with AMC

Next, we deal with tensile tests of aluminium matrix composites (AMC). An AMC material is highly suitable for lightweight applications in aerospace, defense and automotive industry. Compared to monolithic materials, composite materials have advantageous properties such as higher ultimate tensile strength or a higher stiffness to density ratio. Here, we mainly focus on calculating local strains of silicon carbide particle reinforced AMCs from a sequence of scanning electron microscope (SEM) images acquired during tensile tests, where the specimen is pulled in x-direction and elongates with increasing force. Figure 1 illustrates the experimental setup and the resulting image sequence schematically. We are interested in the local deformation behaviour of the composite material on a microscale. Therefore, it is necessary to perform SEM monitored tensile tests to study deformation and crack initiation due to the inhomogeneous microstructure. We will especially focus on ${\epsilon }_{11}$, which describes the change in displacement $u_1$ for the x-direction. A positive value indicates tension and a negative one compression.

Unless stated otherwise, we use the image under no or low load as the reference image $f_1$ and compute the displacement between the reference image and images under higher load. Hence, the displacement and strain are given in the coordinates of the reference image and we overlay the results with the reference image.

The real-world example in Figure 5 illustrates again the difference between the TGV and IC regularized models. In this example, the material has a few cracks in its initial state and they widen up, but the surrounding aluminium matrix deforms smoothly. Whereas the computed displacement $u_1$ looks roughly the same for both methods, the differences are visible in its derivative ${\mathrm{\nabla }}_x{u}_1$. The TGV model only shows some structure where the cracks actually are, but the IC model leads to structures that pass through the whole image since it is, in analogy to the previous artificial example, impossible to split the displacement into two parts for this example. The TGV model outperforms the IC model in this setting. Nevertheless, for the cracks themselves, both methods lead to equally good results as the lower half of the enlarged region in Figure 5 shows. In summary, both methods show the main structures, but the IC model introduces wrong structures. Therefore, we focus on the TGV model in the rest of our examples.

Figure 5. Comparison of IC and TGV models for a real-world example where the cracks widen up but the surrounding aluminium matrix deforms smoothly. (a) Microstructure images of the specimen under low and high load and (b) Results of the TGV and IC models showing the displacement $u_1$ in $\text{μ}\mathrm{m}$ ($1\text{μ}\mathrm{m}=4\mathrm{px}$) and the strain ${\mathrm{\nabla }}_x{u}_1$.

In our next example, we are interested in the detection of different kind of cracks. Figure 6 shows other real-world data arising from a tensile test. Cracks correspond to peaks in ${\mathrm{\nabla }}_x{u}_1$. In the enlarged regions, three different types of cracks are shown, namely

1. a crack due to decohesion of a particle from the surrounding matrix,

2. a crack in the aluminium matrix, and

3. a crack inside a particle.

Figure 6. Detection of different kind of cracks by the TGV model. (a) Microstructure images of the specimen under low and high load and (b) Results of the TGV model showing the displacement $u_1$ in $\text{μ}\mathrm{m}$ ($1\text{μ}\mathrm{m}=16\mathrm{px}$) and the strain ${\mathrm{\nabla }}_x{u}_1$.

Next, we are interested in crack propagation. Figure 7 shows the displacement $u_1$ and the strain ${\epsilon }_{11}={\mathrm{\nabla }}_x{u}_1$ for certain regions under different load. It is remarkable that even under low load, when the cracks are not or hardly visible in the images, the strain in the corresponding regions is large. Thus, it is a sensitive and useful tool to study crack initiation mechanisms.

Figure 7. Results on crack propagation by the TGV model. (a) Microstructure images without load (${\sigma }_0=0\mathrm{MPa}$) and stress–strain curve, (b) Displacement $u_1$ in $\text{μ}\mathrm{m}$ ($1\text{μ}\mathrm{m}=2\mathrm{px}$) and the strain ${\mathrm{\nabla }}_x{u}_1$. The colourbar for the strain is cut off at 0.1 to make smaller values under low load visible, and (c) Two enlarged regions and the strain ${\mathrm{\nabla }}_x{u}_1$ for increasing loads. The colour map for the strain is cut off at 0.1 to make smaller values under low load visible.

4.3. Comparison to correlation-based methods

Now, we draw our attention to a comparison of the proposed method to correlation-based methods. In the following, we use NCorr [3] for a comparison since it is a freely available software package. Note that other state-of-the-art software packages for strain analysis such as Veddac [1] or VIC [2] are based on similar methods. In particular in [36], it is shown that NCorr produces equally good results as VIC. The underlying idea of correlation- based methods is the comparison of windows around each pixel in certain search windows. The window sizes are parameters which must be appropriately chosen. The size of the search window directly influences the computation time and needs to be chosen according to the maximal displacement. The smoothness of the results can be steered by the size of the window around the pixels, similar to the regularization parameters in our model. Due to the extent of this windows, the local resolution of correlation-based methods is limited. Figure 8 shows three examples where the proposed TGV regularized method is compared to NCorr. The first example in Figure 8 shows the same specimen as Figure 5. Although the cracks are relatively large, it is not possible to resolve them with NCorr and there occur some artefacts around them. Thus, the low resolution is especially a drawback in our applications since we are interested in the local behaviour and in particular cracks. Besides the cracks, the smooth parts look very similar to the result using the proposed method. For the cracks, $a_1$ is a good result. A peak is visible at the crack tip, which is the position where the crack is expected to propagate. In contrast, the result of NCorr shows a peak everywhere around the crack.

Figure 8. Comparison of the proposed model to the correlation-based method NCorr for three different image pairs and enlarged regions.

The second and the third example show more complicated deformations. In addition to ${\mathrm{\nabla }}_x{u}_1$, we depict the smooth strain part $a_1$, that is split off by the TGV regularizer. It is visible that the result of NCorr looks very similar to $a_1$, in particular, the same peaks and structures in a degree of $\pm \pi /4$ can be observed. But in addition to this smooth part, which resembles the strain in the classical meaning in materials science, our method is able to resolve also the local damage.

In summary, the proposed method is on the one hand able to extract the information that also commercial software extracts. This warrants that the model in general gives correct results from the viewpoint of materials science. On the other hand, it is possible to visualize the local behaviour with a very high resolution where correlation-based methods fail.

4.4. Results for other materials and experimental settings

Finally, we leave the setting of tensile tests and show some results for other experimental settings and materials.

Figure 9 contains results for a compression test of a fibre reinforced aluminium matrix where load was applied in x-direction. Since the material is compressed in x-direction, cracks are expected to open up in y-direction. Hence, we depict also $u_2$ and ${\mathrm{\nabla }}_y{u}_2$ since these are more suitable for this setting. Also in this case, the proposed method provides reasonable results and resolves cracks.

Figure 9. Results for a compression test of a fibre reinforced aluminium matrix. (a) Microstructure images of the uncompressed and compressed state and (b) Displacement $u_1$, $u_2$ in $\text{μ}\mathrm{m}$ ($10\text{μ}\mathrm{m}=6\mathrm{px}$) and the strain ${\mathrm{\nabla }}_x{u}_1$, ${\mathrm{\nabla }}_y{u}_2$.

Figure 10 shows the results for a fatigue test. Here, we depict the images and results for one load cycle illustrated in Figure 10(a). The first image shows the specimen in state ‘A’ without load ${\sigma }_A=0\mathrm{MPa}$, then load is applied, i.e. state ‘B’ with ${\sigma }_B=632\mathrm{MPa}$. For the third image, i.e. state ‘C’, the load is released again, i.e. ${\sigma }_C=0\mathrm{MPa}$ and then the specimen is compressed to state ‘D’ with ${\sigma }_D=-630\mathrm{MPa}$. For the last image in state ‘E’, the load is released again, i.e. ${\sigma }_E=0\mathrm{MPa}$. We use the image of state ‘A’ as the reference image. The computed results are shown in Figure 10(b). Similar to the previous examples of tensile tests, we see that the strain is positive almost everywhere for the image under load and we see exactly where cracks open up. For state ‘C’, which corresponds to the state after releasing the tension, the strain is mainly positive as well since only the elastic deformation reverses. For the image under compression, the displacement $u_1$ slightly decreases from left to right and the strain is mainly negative. In the last state ‘E’ after completing the load cycle, the specimen is slightly more strained than in the initial state ‘A’. In particular, the crack in the particle opened up a bit. This is visible in the strain component ${\mathrm{\nabla }}_x{u}_1$. Also the stress–strain curve in Figure 10(a) shows a positive strain for state ‘E’.

Figure 10. Results for different states of a fatigue test. An exemplary hysteresis curve is shown in (a), the results for the states ‘A’–‘E’ are shown in (b). (a) Hysteresis curve of one cycle of a fatigue experiment and microstructure image of the initial state ‘A’ and (b) Results for the states shown in (a). Displacement $u_1$ in $\text{μ}\mathrm{m}$ ($1\text{μ}\mathrm{m}=32\mathrm{px}$) and the strain ${\mathrm{\nabla }}_x{u}_1$.

5. Conclusions

In this paper, we proposed a variational optical flow model for computing engineering strains on a microscale. Motivated by the setting of tensile tests, several first and the second order terms were evaluated for regularization. The method of our choice is a TGV regularized model with an optional constraint on the strain based on physical assumptions. The primal–dual method for finding a minimum of the corresponding energy function is especially suited since it directly computes the strain tensor within the iteration process. The results of the proposed model can be used for the detection of cracks and for an analysis of crack initiation and propagation. A comparison to state-of-the-art software used for strain analysis showed that the proposed method clearly outperforms this software, in particular when analysing the local strain. Further, it highlights microstructural damage as an additional benefit. We have shown that, besides tensile tests with AMCs, the method is also applicable for different materials as well as in other experiments such as compression and fatigue tests.

Acknowledgements

The authors thank K. Lichtenberg from the ‘Institute for Applied Materials (IAM)’ at the Karlsruhe Institute of Technology (KIT) for providing the data in Figure 9.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Funding by the German Research Foundation (DFG) within the Research Training Group 1932 ‘Stochastic Models for Innovations in the Engineering Sciences’, project area P3, is gratefully acknowledged.

Appendix 1. Soft and coupled shrinkage

The following two propositions state well-known facts (see e.g. [31]), which are required to compute the proximal maps in Algorithm 1.

Proposition 4

Let $x\in \mathbb{R}$. Then, $\hat{y}=\underset{y\in \mathbb{R}}{\arg \min }\{\lambda |y|+\frac{1}{2}(y-x{)}^2\}$ is given by the soft shrinkage (A1) $\hat{y}=S_{\lambda }\left(x\right)\text{≔}\left\{\begin{array}{ll}0& \mathrm{if}|x|\le \lambda ,\\ x\left(1-\frac{\lambda }{|x|}\right)& \mathrm{if}|x|>\lambda .\end{array}\right.$(A1)

Proposition 5

Let $x\in {\mathbb{R}}^d$. Then, $\hat{y}=\underset{y\in {\mathbb{R}}^d}{\arg \min }\{\lambda \parallel y{\parallel }_2+\frac{1}{2}\parallel y-x{\parallel }_2^2\}$ is given by the grouped or coupled shrinkage $\hat{y}={\mathbf{S}}_{\lambda }\left(x\right)\text{≔}\left\{\begin{array}{ll}0& \mathrm{if}\parallel x{\parallel }_2\le \lambda ,\\ x\left(1-\frac{\lambda }{\parallel x{\parallel }_2}\right)& \mathrm{if}\parallel x{\parallel }_2>\lambda .\end{array}\right.$

Computation of  (12(12) $\hat{y}=\underset{y\in {\mathbb{R}}^2}{\arg \min }\{|\alpha y_1+\beta y_2+\gamma |+\frac{1}{2}\parallel y-x{\parallel }_2^2\},$(12) ).Finally, we compute for given $x\in {\mathbb{R}}^2$ the value $\hat{y}=\underset{y\in {\mathbb{R}}^2}{\arg \min }\{|\alpha y_1+\beta y_2+\gamma |+\frac{1}{2}\parallel x-y{\parallel }_2^2\},$ where $\alpha ,\beta ,\gamma \in \mathbb{R}$. We need to distinguish the following cases:

1. $\alpha =\beta =0$: This gives obviously $\hat{y}=x$.

2. $\alpha =0,\beta \ne 0$: Here, we get ${\hat{y}}_1=x_1$ and $\begin{array}{rl}{\hat{y}}_2& =\underset{y_2\in \mathbb{R}}{\arg \min }\left\{|\beta y_2+\gamma |+\frac{1}{2}(x_2-y_2{)}^2\right\}\\ & =\underset{y_2\in \mathbb{R}}{\arg \min }\left\{|y_2+\frac{\gamma }{\beta }|+\frac{1}{2|\beta |}(x_2-y_2{)}^2\right\}\\ & =S_{|\beta |}\left(x_2+\frac{\gamma }{\beta }\right)-\frac{\gamma }{\beta }\end{array}$ with the soft shrinkage operator $S_{|\beta |}$

   For $\alpha \ne 0,\beta =0$, we get analogously $({\hat{y}}_1,{\hat{y}}_2)=\left(S_{|\alpha |}\left(x_1+\frac{\gamma }{\alpha }\right)-\frac{\gamma }{\alpha },x_2\right).$

3. $\alpha ,\beta \ne 0$: By substitution, we get the following equivalent minimization problem: $({\hat{z}}_1,{\hat{z}}_2)=\underset{z\in {\mathbb{R}}^2}{\arg \min }\left\{|z_1+z_2|+\frac{1}{2{\lambda }_1}(z_1-{\tilde{x}}_1{)}^2+\frac{1}{2{\lambda }_2}(z_2-{\tilde{x}}_2{)}^2\right\}$ with $\begin{array}{rl}{z}_1& \text{≔}\alpha y_1,{\tilde{x}}_1\text{≔}\alpha x_1,{\lambda }_1\text{≔}{\alpha }^2>0,\\ z_2& \text{≔}\beta y_2+\gamma ,{\tilde{x}}_2\text{≔}\beta x_2+\gamma ,{\lambda }_2\text{≔}{\beta }^2>0.\end{array}$

The following proposition provides a solution to the minimization problem.

Proposition 6

Let $(x_1,x_2{)}^{\mathrm{T}}\in {\mathbb{R}}^2$. Then, the unique minimizer of $F(z_1,z_2)\text{≔}|z_1+z_2|+\frac{1}{2{\lambda }_1}(z_1-{\tilde{x}}_1{)}^2+\frac{1}{2{\lambda }_2}(z_2-{\tilde{x}}_2{)}^2$ is given by $({\hat{z}}_1,{\hat{z}}_2)=\left\{\begin{array}{ll}({\tilde{x}}_1-{\lambda }_1,{\tilde{x}}_2-{\lambda }_2)& \mathrm{if}{\tilde{x}}_1+{\tilde{x}}_2>{\lambda }_1+{\lambda }_2,\\ ({\tilde{x}}_1+{\lambda }_1,{\tilde{x}}_2+{\lambda }_2)& \mathrm{if}{\tilde{x}}_1+{\tilde{x}}_2<-({\lambda }_1+{\lambda }_2),\\ \left(\frac{{\tilde{x}}_1{\lambda }_2-{\lambda }_1{\tilde{x}}_2}{{\lambda }_1+{\lambda }_2},\frac{{\tilde{x}}_2{\lambda }_1-{\lambda }_2{\tilde{x}}_1}{{\lambda }_1+{\lambda }_2}\right)& \mathrm{otherwise}.\end{array}\right.$

Proof.

Since F is not differentiable at $z_1=-z_2$, we distinguish whether ${\hat{z}}_1+{\hat{z}}_2=0$ or not.

1. ${\hat{z}}_1+{\hat{z}}_2\ne 0$: First, we assume that ${\hat{z}}_1+{\hat{z}}_2>0$. Since F is differentiable around $\hat{z}$, we can compute the gradient and set it to zero, i.e. $0=\left(\begin{array}{c}1\\ 1\end{array}\right)+\left(\begin{array}{c}\frac{1}{{\lambda }_1}({\hat{z}}_1-{\tilde{x}}_1)\\ \frac{1}{{\lambda }_2}({\hat{z}}_2-{\tilde{x}}_2)\end{array}\right).$ Hence, we obtain $({\hat{z}}_1,{\hat{z}}_2)=({\tilde{x}}_1-{\lambda }_1,{\tilde{x}}_2-{\lambda }_2)$ and ${\hat{z}}_1+{\hat{z}}_2>0$ if and only if ${\tilde{x}}_1+{\tilde{x}}_2>{\lambda }_1+{\lambda }_2$. Analogously, one can show $({\hat{z}}_1,{\hat{z}}_2)=({\tilde{x}}_1+{\lambda }_1,{\tilde{x}}_2+{\lambda }_2)$ for ${\tilde{x}}_1+{\tilde{x}}_2<-({\lambda }_1+{\lambda }_2)$.

2. ${\hat{z}}_1+{\hat{z}}_2=0$: Then ${\hat{z}}_1=-{\hat{z}}_2$ and we obtain ${\hat{z}}_1=\underset{z_1\in \mathbb{R}}{\arg \min }\left\{\frac{1}{2{\lambda }_1}(z_1-{\tilde{x}}_1{)}^2+\frac{1}{2{\lambda }_2}(z_1+{\tilde{x}}_2{)}^2\right\}.$ Setting the derivative to zero leads to $0=\frac{1}{{\lambda }_1}({\hat{z}}_1-{\tilde{x}}_1)+\frac{1}{{\lambda }_2}({\hat{z}}_1+{\tilde{x}}_2).$ Thus, the solution is given by ${\hat{z}}_1=\frac{{\tilde{x}}_1{\lambda }_2-{\lambda }_1{\tilde{x}}_2}{{\lambda }_1+{\lambda }_2}\mathrm{and}{\hat{z}}_2=\frac{{\tilde{x}}_2{\lambda }_1-{\lambda }_2{\tilde{x}}_1}{{\lambda }_1+{\lambda }_2}.\text{▪}$

Remark 7

The soft shrinkage (A1(H1) $\lambda {\mathrm{H}}^1\left(u\right)\text{≔}\lambda \parallel \mathbf{\nabla }u{\parallel }_2^2,$(H1) ) is obtained in the special case ${\lambda }_1={\lambda }_2$ and ${\tilde{x}}_1={\tilde{x}}_2$.