# PGD et mesure de champs par corrélation d'images

Jean-Charles Passieux, Luis Perini, Emmanuel de Luycker, Jean-Noël Périé

► To cite this version:

Jean-Charles Passieux, Luis Perini, Emmanuel de Luycker, Jean-Noël Périé. PGD et mesure de champs par corrélation d'images. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717111

HAL Id: hal-01717111

https://hal.science/hal-01717111v1

Submitted on 25 Feb 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_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# PGD et mesure de champs par corrélation d'images

Jean-Charles PASSIEUX, Luis A. GOMES PERINI, Emmanuel DE LUYCKER, Jean-Noël PÉRIÉ

Université de Toulouse, Institut Clément Ader, EA 814 - Toulouse, France, passieux@insa-toulouse.fr
gomesper@insa-toulouse.fr,
emmanuel.de-luycker@iut-tarbes.fr,
jean-noel.perie@iut-tlse3.fr

**Résumé** — Un algorithme alternatif est proposé pour la mesure *globale* de champs de déplacements continus par corrélation d'images numériques. L'idée consiste à utiliser la *Proper Generalized decomposition* (PGD) pour résoudre le problème non-linéaire de la corrélation. En effet, l'incrément de déplacement est recherché sous la forme d'un produit de fonctions de variables séparées. Les problèmes dans les 3 différentes dimensions de l'espace sont découplées par un algorithme de point fixe basique utilisé classiquement en PGD. La méthode s'applique tant à des images 2D qu'a des images 3D volumiques (issues de tomographie, d'IRM...). L'incertitude de mesure est comparable à celle des approches éléments finis, tandis que le temps de calcul est réduit particulièrement quand on descend en résolution. Enfin, dans le cadre de mesure de champs de déformations, cet algorithme est bien adapté à la nature tensorielle des fonctions d'interpolation B-Splines.

**Mots clés** — CIN globale, séparation de variables, analyse isogéométrique

## 1 Introduction

La mesure de champs par corrélation d'images numériques (CIN) est aujourd'hui très largement utilisée en mécanique des solides [1]. L'idée consiste à imager une éprouvette avant et après une transformation imposée par un chargement :

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Images en niveaux de gris : (gauche) image de référence $f(x)$, (droite) image sous charge $g(x)$.</figcaption></figure>

Même si la plupart du temps, la CIN est basée sur une approche dite *locale*, une alternative récente consiste à rechercher le champ déplacement d'une manière *globale*, par une méthode de type Galerkin [2]. Cette dernière est basée sur une formulation faible de l'équation de conservation des niveaux de gris
entre une image de référence $f(\mathbf{x})$ et une image déformée $g(\mathbf{x})$ [3] :

$$ \mathbf{u}(\mathbf{x}) = \arg\min_{\Omega} \int_{\Omega} [f(\mathbf{x}) - g(\mathbf{x} + \mathbf{u}(\mathbf{x}))]^2 d\mathbf{x} \quad (1) $$

L'inconnue **u** étant dans un espace de dimension infinie, une approximation est recherchée dans un sous espace engendré par une base de fonctions d'interpolation $\varphi_i(\mathbf{x})$

$$ \mathbf{u}(\mathbf{x}) = \sum_{i} \varphi_{i}(\mathbf{x})x_{i} \quad (2) $$

Grâce à ce formalisme, le choix du sous-espace d'approximation peut être très varié. Elle permet par exemple de généraliser les approches locales en considérant des fonctions constantes ou linéaires par morceaux sur des sous ensembles de pixels avoisinants un ensemble de points [1]. Une interpolation élément fini a été proposée plus récement [4, 2, 5] qui consiste à imposer la continuité du champs de déplacement ce qui a un effet régularisant sur la mesure ; ainsi les incertitudes de mesure sont réduites. Cette interpolation a également l'intérêt de créer une passerelle entre la simulation et la mesure [2]. Des fonctions d'interpolation plus régulières (B-Splines, NURBS) ont été utilisées notamment pour la mesure de gradients du champs de déplacement [6, 7]. D'autres types de fonctions (analytiques, mécaniques, et numériques) peuvent également être utilisées.

Une fois l'interpolation choisi, la résolution du problème non-linéaire (1) conduit à l'inversion de plusieurs systèmes fournissent la solution sur l'ensemble de la région d'intérêt.

$$ \mathbf{M}\mathbf{x}^n = \mathbf{b}^n \quad (3) $$

avec

$$ \mathbf{M}_{ij} = \int \varphi_i(\mathbf{x})^T \nabla f \nabla^T f \varphi_j(\mathbf{x}) d\mathbf{x} \quad (4) $$

$$ \mathbf{b}_i = \int \varphi_i(\mathbf{x})^T \nabla f (f(\mathbf{x}) - g(\mathbf{x} + \mathbf{u}^{n-1}(\mathbf{x}))) d\mathbf{x} \quad (5) $$

Dans le cas des approches dites locales de la corrélation d'image, l'opérateur **M** est diagonal par bloc, ce qui conduit à un coût de résolution infime. C'est sans doute ce qui a participé à la popularité de ces approches. A l'inverse, dans le cas d'une interpolation élément fini, l'opérateur **M** est certes creux, mais c'est une matrice bande. Si le nombre d'inconnues augmente (en particulier pour un maillage haute résolution), le coût d'inversion de ces systèmes peut devenir significatif. Dans la suite, nous proposons d'associer à la corrélation d'image un solveur PGD [8, 9] pour réduire le coût de cette étape.

# 2 Corrélation d'images par PGD

## 2.1 Principe

Dans ce cadre global, une méthode de CIN basée sur la séparation des variables spatiales par PGD a été proposée récemment [10]. En effet, contrairement à une méthode CIN-EF standard, le champ de déplacement est recherché sous la forme du produit de fonctions de variables séparées [11] :

$$ \mathbf{u}(x, y, z) = \begin{pmatrix} u(x, y, z) \\ v(x, y, z) \\ w(x, y, z) \end{pmatrix} = \sum_{i=1}^{m} \begin{pmatrix} u_i^x(x) \cdot u_i^y(y) \cdot u_i^z(z) \\ v_i^x(x) \cdot v_i^y(y) \cdot v_i^z(z) \\ w_i^x(x) \cdot w_i^y(y) \cdot w_i^z(z) \end{pmatrix} \quad (6) $$

dans laquelle les fonctions $u_i^x(x), u_i^y(y), u_i^z(z), v_i^x(x), v_i^y(y), v_i^z(z), w_i^x(x), w_i^y(y), w_i^z(z)$ sont les inconnues du problème et $m$ le rang de l'approximation. On choisit de construire une meilleure approximation de rang 1 à chaque prédiction linéaire du problème de corrélation. Pour se faire, elle est introduite dans l'équation (1). Le problème est résolu par un algorithme de point fixe standard ce qui permet de séparer les problèmes selon les différentes dimensions de l'espace.

Des exemples synthétiques et réalistes à base de fonctions de forme élément fini seront proposées afin de démontrer ses performances par rapport à une méthode CIN-Q4 prise comme référence. Avec la méthode proposée, le nombre d'inconnues dans chaque dimension a beaucoup moins d'impact sur le temps de calcul. On montrera que cette méthode est donc un bon candidat pour faire des mesures à haute résolution.
## 2.2 Quantification des incertitudes de mesure

Dans cette section, une analyse de performances a priori est menée sur un cas test 2D synthetisé à partir de bruit de perlin et d'un champs de déplacement constant. En suivant la démarche proposée par [2], une série de p images $g_j$ sont générées par l'advection de f par une translation de corps rigide dans la direction x. L'amplitude $u_j^{ref}$ de ce déplacement imposé prends différentes valeurs entre 0 et 1 pixel. Une technique basée sur un shift dans l'espace de Fourier es utilisée pour générer les image transportées. Par suite, l'analyse de corrélation est réalisée sur toutes les paires $(f, g_j)$. La précision de la mesure du déplacement $u_j$ est quantifiée par deux moments statistiques. L'erreur moyenne sur le déplacement $\langle \delta_u \rangle$ et l'écart type du déplacement sensé être constant $\langle \sigma_u \rangle$ et définis par :

$$ \langle \delta_u \rangle = \frac{1}{p} \sum_{j=1}^{p} |\langle u_j \rangle - u_j^{ref}| $$

$$ \langle \sigma_u \rangle = \frac{1}{p} \sum_{j=1}^{p} \langle (u_j - \langle u_j \rangle)^2 \rangle^{1/2} $$

L'évolution de ces indicateurs en fonction du nombre de pixel dans chaque éléments sont représentés sur la Figure 2a et 2b. Il sont comparés aux même quantités, mais pour une approche éléments fini traditionnelle Q4, avec le même maillage, les même images et les même schémas d'interpolation des niveaux de gris.

<figure><img src="image_3.png" /><figcaption>(a) Standard displacement uncertainty $\langle \sigma_u \rangle$ as a function of the element size $h$ in Q4-DIC (○) and PGD-DIC(+).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Mean displacement error $\langle \delta_u \rangle$ as a function of the element size $h$ in Q4-DIC (○) and PGD-DIC(+).</figcaption></figure>

<caption>Fig. 2 – Quantification d'incertitudes a priori.</caption>

On peut observer que la méthode proposée possède, à peu de choses près, le même performances que Q4-DIC, lorsque l'on observe l'incertitude de mesure et le biais. Ce résultat n'est pas supprenant dans la mesure ou les interpolations du champs de déplacement sont les mêmes, seul l'algorithme diffère. Cet algorithme représente donc une alternative à Q4-dic qui nous sert de référence. Une comparaison entre les méthodes de corrélation d'images éléments finis et les méthodes dites locales à été proposée récement dans [12].

## 2.3 Application à un essai biaxial sur une éprouvette composite Carbone/Carbone

Une éprouvette cruciforme faite en composite Carbone/Carbone et soumise à un chargement biaxial plan est analysée avec la méthode. La méthode est comparée à une méthode éléments finis Q4 qui nous servira de référence. Cet example est étudiant dan [13] et plus récement dans [14] pour l'identification d'un loi d'endommagement. Le déplacement est mesurée dans une région d'intérêt de 850×780 pixels. Des éléments de 24 pixels sont utilisés dans chaque directions. L'image de référence et les maillage utilisés sont présentés sur la Figure 3. Le problème est résolu avec la méthode proposée PGD-DIC. L'amplitude du déplacement est affichée sur le domaine dans sa configuration déformée sur la Figure 3(b) et le résidu de corrélation associé est illustré sur la figure 3(c). Les solutions obtenues en utilisant la PGD est reconstruite et compatée à une solution de référence (Fig. 4) obtenue avec Q4-DIC basé sur la même ROI avec le même nombre de pixels par éléments. Il apparaît que les solutions Q4 et PGD
<figure><img src="image_5.png" /><figcaption>Fig. 3 – Application de la PGD-DIC à un essai réel : éprouvette cruciforme en composite C/C soumise à une chargement biaxial plan. Le mouchetis peint en noir et blanc (a). Le champs de déplacement mesuré (b) et la carte des résidus de corrélation en pourventage de la dynamique (c)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 – Comparaison de la composante y du déplacement mesuré avec PGD-DIC et Q4-DIC. La même précision est obtenue sur la composante $u_x$.</figcaption></figure>

sont quasiment les mêmes (Fig. 4). L'écart entre des deux solutions est de l'ordre de grandeur du critère d'arrêt de l'algorithme itératif de la corrélation (Fig. 4).

## 2.4 PGD multigrille

Par ailleurs, une version multi-échelle du solveur PGD est proposée ici. La PGD approxime le champs de déplacement inconnu par une série. Ici, on s'autorise d'utiliser des discrétisations différentes pour chacun des termes de la série. Ainsi, la stratégie consiste à utiliser des discrétisations relativement grossières pendant les première itérations et progressivement, les maillages dans chaque dimensions s'affinent. Cet algorithme PGD multiéchelle est associé à un dérafinement des images, dont il est montré qu'il permet de régulariser le problème en évitant les bassins d'attraction associés à des minimum locaux [2, 5, 15]. La convergence en est donc améliorée, particulièrement dans le cas de déplacement d'aplitude grande devant la taille des éléments.

La méthode est appliquée ici au cas de mesures de champs dans le volume (CVN) à partir d'images 3D d'Imagerie par Résonance Magnétique (IRM). L'analyse d'une paire d'images IRM d'un essai sur un os bovin in situ réalisé par [16]. Il s'agit d'un essai de compression sur un os dont la structure trabéculaire sert de texture naturelle pour la corrélation d'image (cf. figure ci-dessous).

Le déplacement mesuré et la carte de résidu de corrélation est présentée en Figure 5(a) et 5(b), avec 20 voxels par éléments dans chaque dimensions. Par ailleurs, sur la Figure 5(c) on peut remarquer que l'algorithme PGD monoéchelle ne converge pas avec une taille de 15 voxels par élément, alors que lorsque 3 niveaux de grilles sont utilisés, le solveur atteint le niveau de précision requis.
<figure><img src="image_7.png" /><figcaption>Fig. 5 – Exemple de mesure de déplacement par PGD-CVN à partir d'images IRM d'un os trabéculaire : composante selon x du champ de déplacement (gauche) et carte seuillée des résidus de correlation (milieu) et convergence de l'algorithme multigrille droite)</figcaption></figure>

## 2.5 Quelle interpolation pour la mesure de déformations

Dans le formalisme global de la CIN, une interpolation basée sur l'analyse isogéométrique [17, 6, 7] s'est montrée très efficace pour minimiser la sensibilité au bruit et les incertitudes de mesure [7], comparativement aux éléments finis, notamment pour la mesure de champs de déformations. Un élément NURBS en 2D ou 3D est construit par le produit tensoriel de fonctions B-Splines unidimensionnelles [17, 18]. La PGD semble donc être un algorithme naturel pour ce type d'interpolation.

Pour illustrer ceci, on choisi de caractériser la qualité de la mesure du déplacement et de la déformation à partir d'un cas synthétique. Une image de référence est transportée par le champs de déplacement analytique suivant [18] :

$$u_s(x) = \sin\left(\frac{5\pi}{2}\right)\left(1 - \left|2\frac{x}{x_{max}} - 1\right|\right)$$

Sur la figure 6, le déplacement mesuré est comparé au déplacement analytique de référence, pour différents types d'interpolation : des éléments finis linéaires, des éléments finis quadratiques, des B-Splines $C^{p-1}$ de degré 3, et des polynômes de Hermite. On peut voir que quelque soit la base d'interpolation du champs de déplacement, le déplacement est continu et relativement bien mesuré.

En revanche, lorsque l'on compare l'évolution de la dérivée du déplacement, les méthodes éléments finis p1 et p2 concuisent à une déformation discontinu au bord des éléments. Les fonctions B-Spline se montrent plus régularisantes, comme montré dans [18]. On se propose également de comparer ces trois premières approches à une interpolation de hermite, afin d'assurer la continuité des dérivées du déplacement. L'intérêt de ces fonctions est que le support des fonctions d'interpolation sont constantes, ce qui n'est pas le cas avec les B-Splines.

Enfin nous appliquons cette méthode à la mesure de champs de déformation sur l'éprouvette en composite stratifié de la figure 1. Le champs de glissement $ε_{12}$ est représenté sur la figure 8 pour des éléments finis linéaires et pour des éléments B-Splines. Dans [18], il est montré que ces fonctions plus régulières améliorent la mesure de déformation. Cependant, ces fonctions en 2D ou 3D sont bien souvent construites par produit tensoriel des fonction 1D. En ceci l'algorithme PGD-DIC, parfaitement adapté à la nature tensorielle de ces fonctions, représente une alternative intéressante.

## 3 Conclusion

Une algorithme de corrélation d'image numérique (CIN) alternatif est proposé. Il est basé sur un sol-
veur proper generalized decomposition (PGD). La méthode conserve les avantage des méthodes éléments
finis standar (modularité, continuité, lien avec la simulation...) en limitant son principal inconvénient, le
coût de calcul. En effet, quelque soit la dimension du problème(2D or 3D), la complexité des problèmes
mis en jeux par la PGD sont linéaires. Avec une telle propriété la méthode est une algorithme alternatif
crédible, particulièrement en 3D dans le cas de mesure de champs volumique (DVC).

Une version multigrille ce solveur associée à un dérafinement des images permet d'améliorer la
convergence de corrélation d'image, à l'image des travaux antérieurs [?, 2]. Des illustrations 3D per-
<figure><img src="image_8.png" /><figcaption>Fig. 6 – Mesure du champ de déplacement $u_x$ avec la méthode proposée : avec, de gauche à droite : des éléments finis linéaires, des éléments finis quadratiques, des B-Splines $C^{p-1}$ de degré 3, et des polynômes de Hermite.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 – Mesure du champ de déformation $\epsilon_{xx}$ avec la méthode proposée : avec, de gauche à droite : des éléments finis linéaires, des éléments finis quadratiques, des B-Splines $C^{p-1}$ de degré 3, et des polynômes de Hermite.</figcaption></figure>

mettent de mettre en évidence cet effet sur des images issues d'un essai réel.

Comme montré dans [18], l'utilisation de fonction plus régulières peuvent améliorer la mesure de déformation. Ces fonctions de haut degré de régularité en 2D ou 2D sont construites par produit tensoriel des fonction 1D. En ceci l'algorithme proposé est donc naturellement bien adapté à ce type d'interpolations.
<figure><img src="image_10.png" /><figcaption>Fig. 8 – Mesure du champ de déformation ε<sub>xy</sub> avec la méthode proposée : avec des éléments finis, 18 px par éléments (gauche) et des B-Splines, 30 px par éléments (droite)</figcaption></figure>

tion.

Références

[1] M.A. Sutton, J.-J. Orteu, and H. Schreier. Image correlation for shape, motion and deformation measure-
ments : Basic Concepts, Theory and Applications. Springer, New York, NY (USA), 2009.

[2] G. Besnard, F. Hild, and S. Roux. “finite-element” displacement fields analysis from digital images : Appli-
cation to portevin–le châtelier bands. *Experimental Mechanics*, 46(6) :789–803, 2006.

[3] B.K.P. Horn and G. Schunck. Determining optical flow. *Artificial Intelligence*, 17 :185–203, 1981.

[4] Y. Sun, J.H.L. Pang, C. Khuen Wong, and F. Su. Finite element formulation for a digital image correlation method. *Applied Optics*, 44(34):7357–7363, 2005.

[5] J. Fehrenbach and M. Masmoudi. A fast algorithm for image registration. *Comptes Rendus Mathématique*, 346(9-10):593–598, 2008.

[6] P. Cheng, M. Sutton, H.W. Schreier, and S. R. McNeill. Full-field speckle pattern image correlation with b-spline deformation function. *Experimental mechanics*, 42(3):344–352, 2002.

[7] J. Réthoré, S. Roux, and F. Hild. Hybrid analytical and extended finite element method (HAX-FEM) : A new enrichment procedure for cracked solids. *International Journal for Numerical Methods in Engineering*, 81(3):269–285, 2010.

[8] E. Cueto F. Chinesta, P. Ladevèze. A short review on model order reduction based on proper generalized decomposition. *Arch Comput Methods Eng*, 18 :395–404, 2011.

[9] P. Ladevèze, J.-C. Passieux, and D. Néron. The LATIN multiscale computational method and the proper generalized decomposition. *Comput. Methods Appl. Mech. Engn.*, 199(21):1287–1296, 2010.

[10] J.-C. Passieux and J.-N. Périé. High resolution digital image correlation using proper generalized decompo-
sition : PGD-DIC. *International Journal for Numerical Methods in Engineering*, 92(6):531–550, 2012.

[11] B. Bognet, F. Bordeu, F. Chinesta, A. Leygue, and A. Poitou. Advanced simulation of models defined in plate geometries : 3D solutions with 2D computational complexity. *Computer Methods in Applied Mechanics and Engineering*, 201-204(1):1–12, 2012.

[12] F. Hild and S. Roux. Comparison of local and global approaches to digital image correlation. *Experimental Mechanics*, 52(9):1503–1519, 2012.

[13] J.-N. Périé, S. Calloch, C. Cluzel, and F. Hild. Analysis of a multiaxial test on a c/c composite by using digital image correlation and a damage model. *Experimental Mechanics*, 42:318–328, 2002.

[14] J.-N. Périé, H. Leclerc, S. Roux, and F. Hild. Digital image correlation and biaxial test on composite material for anisotropic damage law identification. *International Journal of Solids and Structures*, 46(11-12):2388–2396, 2009.

[15] J. Réthoré, F. Hild, and S. Roux. Shear-band capturing using a multiscale extended digital image correlation technique. *Comput Methods Appl Mech Eng*, 196(49–52):5016–5030, 2007.

[16] A. Benoit, S. Guérard, B. Gillet, G. Guillot, F. Hild, D. Mitton, J.-N. Périé, and S. Roux. 3d analysis from micro-mri during in situ compression on cancellous bone. *Journal of Biomechanics*, 42(14):2381–2386, 2009.

[17] T. J. R. Hughes, J.A. Cottrell, and Y. Bazilevs. Isogeometric analysis : Cad, finite elements, nurbs, exact geometry and mesh refinement. *Computer Methods in Applied Mechanics and Engineering*, 194 :4135–4195, 2005.

[18] T. Elguedj, J. Réthoré, and A. Buteri. Isogeometric analysis for strain field measurements. *Computer Methods in Applied Mechanics and Engineering*, 200(1-4):40–56, 2010.