Medical imaging using elastography

Jérôme Fehrenbach, Davuth Ia, Mohamed Masmoudi, Rémi Souchon, Philippe

Trompette

► To cite this version:

Jérôme Fehrenbach, Davuth Ia, Mohamed Masmoudi, Rémi Souchon, Philippe Trompette. Medical imaging using elastography. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812991

HAL Id: hal-01812991

https://hal.science/hal-01812991

Submitted on 12 Jun 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.

Public Domain
Medical imaging using elastography

Jérôme Fehrenbach†, Davuth Ia†, Mohamed Masmoudi†,
Rémi Souchon‡, Philippe Trompette‡

† Laboratoire MIP UMR 5640, 118 route de Narbonne, 31062 Toulouse cedex 4

‡ Unité de recherche INSERM U556, 151 cours Albert Thomas, 69424 Lyon cedex 03

RÉSUMÉ. Le tenseur de Hooke d'un matériau élastique (linéaire isotrope) dépend de deux paramètres : le module d'Young E et le coefficient de Poisson ν. Un problème intéressant - du point de vue médical - est de détecter des inhomogénéités où le module d'Young prend une valeur différente du reste.

On propose ici une méthode pour reconstruire une approximation du module d'Young, en sup-
posant qu'il est constant sauf dans certaines régions où il prend une valeur différente. L'outil
principal est une méthode générale pour les problèmes inverses : il s'agit d'une implémentation
économne en mémoire et en calculs de la méthode de Gauss-Newton, basée sur l'utilisation de
la dérivation algorithmique en mode direct et en mode adjoint.

Cette méthode est validée par des résultats expérimentaux : la propriété de régularisation de l'algorithme de Gauss-Newton permet de localiser les plus grosses inhomogénéités.

ABSTRACT. The Hooke tensor of a linear isotropic elastic material depends on two parameters : the Young modulus E and the Poisson coefficient ν. An interesting problem - from the medical point of view - is to detect inhomogeneities where the Young modulus takes a different value from the background.

A method is proposed here to reconstruct an approximation of the Young modulus, assuming it is constant except in some regions where it takes a different value. The main tool is a general method for inverse problems : it is an implementation of Gauss-Newton's method that uses few memory and few computations, based on the use of direct and adjoint derivative.

This method is validated with experimental results : the regularization property of Gauss-
Newton's algorithm allows to locate the larger inhomogeneities.

MOTS-CLÉS : élastographie, imagerie médicale, optimisation topologique, problème inverse, élasticité, échographie.

KEYWORDS: elastography, medical imaging, topological optimization, inverse problems, elasticity, echography.
# 1. Introduction

Prostate and breast tumors can have a Young modulus much higher than the surrounding safe tissue. Their detection by clinical palpation requires that the hard nucleus have to be big or near enough from the skin. Elastography was introduced by J.Ophir et al. [OPH 91] to provide an elasticity image by comparison of two sonograms before and after a small external compression. The main drawback of this technique is that elastograms are noisy and difficult to read. A recent review article on elastography is [OPH 02].

The inverse problem framework [KIR 96] is used to improve elastograms. Gauss Newton algorithm is a regularization method for inverse problems, it chooses a solution of minimal norm. Indeed, it is shown here that the Gauss Newton method can be implemented thanks to both forward and reverse modes of algorithmic differentiation. This implementation of the Gauss Newton algorithm adds data from the physical model to go through the lack of measured data, it has the advantage of being no more expensive than the gradient method and gives a better convergence.

In this paper, the direct problem for elasticity is recalled (section 2). We explain then (section 3) the strategy of the inverse problem : reconstructing the Young modulus from the radial displacement under known boundary conditions. The implementation of Gauss-Newton algorithm using forward and reverse modes of algorithmic differentiation is then split in short algorithms. In the last section, experimental results are given.

# 2. Direct problem

Consider a smooth domain $\Omega$ in the plane. The boundary of $\Omega$ is divided into two parts : $\partial\Omega = \Gamma_N \cup \Gamma_D$. The domain $\Omega$ is filled with an elastic material, subject to a displacement $u \in L^2(\Omega, \mathbf{R}^2)$. There are no volume forces, there is a known displacement $d$ on $\Gamma_D$ and known forces $f$ on $\Gamma_N$. The Lamé coefficients of the material $\lambda = \lambda(x) = \lambda_0 E(x)$ and $\mu = \mu(x) = \mu_0 E(x)$ depend on the space variable (this is equivalent to assuming the Poisson ratio $\nu$ is constant, and the Young modulus $E(x)$ depends on the space variable).

The boundary values problem is to find a displacement field $u$ and a stress field $\sigma$ defined in $\Omega$ by :

$$ (BP) \left\{ \begin{array}{l} \sigma = \lambda \operatorname{tr} \epsilon(u) I + 2\mu \epsilon(u) \quad \Omega \\ -\operatorname{div} \sigma = 0 \quad \Omega \\ \sigma.n = f \quad \Gamma_N \\ u = d \quad \Gamma_D \end{array} \right. $$
Perform a degree one finite elements method : the finite elements version of the mixed boundary problem (BP) is equivalent to the linear system :

$$Aq = b, \qquad [1]$$

where $q \in \mathbb{R}^{2n}$ is the displacement vector, $b \in \mathbb{R}^{2n}$ is defined by the boundary conditions and $A$ is the stiffness matrix.

The stiffness matrix $A$ depends on the Young modulus distribution $E$, this matrix is denoted $A(E)$. Note that $E \mapsto A(E)$ is an affine map. Let $\tilde{A}$ be its linear part.

### 3. Inverse problem

We look for inhomogeneities having Young modulus $E_2$ in a material having Young modulus $E_1$. The known measure is the radial displacement under known boundary conditions, and the unknown is the location of the inhomogeneities.

The radial displacement at a point $M$ is given by $u_r = L.u$, where $L$ is a linear operator. The $n \times 2n$ matrix of the discretized version of this operator can be easily formed, this matrix is also denoted $L$.

We want to minimize the following quantity :

$$j(E) = \frac{1}{2} \|L.u_E - u_r^{\text{measur.}}\|^2 = \frac{1}{2} \|F(E)\|^2,$$

where $F(E) = L.u_E - u_r^{\text{measur.}}$, with $u_E$ solution of $A(E)u_E = b$.

The idea is to apply one step of the Gauss-Newton algorithm, starting from an homogeneous distribution $E \equiv E_1$ : we shall compute an iterate of Newton $E + d$ with $d$ solution of :

$$DF^T D F d = -DF^T F.$$

The vector $d$ will be searched by the conjugate gradient method. Indeed, the conjugate gradient method needs only to know the product of the given matrix by a vector, avoiding thus to compute the whole jacobian matrix, a matrix-matrix multiplication and the use of memory to store a (non sparse) matrix.

For the computation of $DF^T D F x$ we use the algorithmic differentiation rules, and proceed as follows :

1) **Computing $z = DFx$**

It is the directional derivative of a vectorial function : use forward differentiation.

2) **Computing $DF^T z$**

It is a scalar criterion : use reverse mode of algorithmic differentiation.

**Proposition 1 (Computing $z = DFx$)**

$$z = L\delta$$
where $\delta$ is solution of the linear system

$$A(E)\delta = -\hat{A}(x)u_E,$$

where $\hat{A}(x) = A(x) - A(0)$.

**Proposition 2 (Computing $DF^Tz$)** Let $v_E \in \mathbf{R}^{2n}$ be the solution of the linear system

$$A^T v = -L^T z,$$

Then for any vector $h \in \mathbf{R}^n$:

$$(DF^T)z.h = [D(F^Tz)].h = (\hat{A}(h)u_E | v_E)$$

**Corollary 1 (Computing $DF^T F$)** Let $w_E \in \mathbf{R}^{2n}$ be the solution of the linear system

$$A^T w = -L^T F, \qquad [2]$$

Then

$$\forall h \in \mathbf{R}^n, \quad (DF^T)F.h = (\hat{A}(h)u_E | w_E).$$

**Algorithm 1 : Finding the Young modulus distribution :**

**input** : mesh, $E_1$, $\nu$, boundary conditions, $u_r^{\text{measur.}}$.

* construct the matrix $L$ from mesh data,

* compute the stiffness matrix $A$ assuming homogeneous Young modulus $E \equiv E_1$, compute the direct state $u_E$,

* compute the adjoint state $w_E$,

* compute $DF^T F$ using equation [2]

* solve for $d$ the equation $DF^TDFd = -DF^T F$ using the conjugate gradient method, the multiplication of a vector by $DF^TDF$ is given by algorithm 2

* the estimated Young modulus distribution is $E + d$.

**Algorithm 2 : Computing $y = DF^TDFx$ for a given vector $x$**

**input** : mesh, $\nu$, stiffness matrix $A(E)$, direct state $u_E$, matrix $L$

* solve for $\delta$ the equation $A(E)\delta = -\hat{A}(x)u_E$,

* compute $z = L.\delta$ (see proposition 1),

* solve for $v_E$ the equation $A^T v_E = L^T z$,

* the components ($y_i$) of $y$ are given by $y_i = (A(e_i)u_E|v_E)$ (see proposition 2).
<figure><img src="image_1.png" /><figcaption>**Figure 1.** *The phantom experiment*</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** *The gradient*</figcaption></figure>

**4. Experimental results**

We have applied the algorithm above to experimental data obtained by R.Souchon [SOU 04] : inside a phantom of gelatin including six inhomogeneities (of different sizes), a ball is inflated (see fig.1) : the measures are the radial displacements between before and after inflation. The imaging system is based on an ultrasound scanner (Combison 311, Kretz, Austria) equipped with a transrectal probe (IRW 77AK, Kretz, Austria).

The data are treated as follows :

1) the boundary conditions are not exactly known, but reconstructed using an efficient implementation of Gauss-Newton algorithm similar to the one explained above (but simpler) : we look for Dirichlet conditions on the inner boundary, and Neumann conditions on the outer boundary.

In order to show that the gradient of the cost functional gives no relevant information, we show in fig. 2 the value of the gradient. The inhomogeneities are located inside the gelatin, not on the boundary.

2) once the boundary conditions are estimated, algorithm 1 is applied to estimate the Young modulus distribution. Indeed, a Tikhonov regularization term is added to the cost functional : $J(E) = \frac{1}{2}||F(E)||^2 + \alpha||E||^2$.

The results are shown in fig.3. We can distinguish four among the inclusions, the other two inclusions are smaller and do not appear on our results.

**5. References**

KIR 96 Kirsch A, *An introduction to the mathematical theory of inverse problems*,
(Springer, 1996)
<figure><img src="image_3.png" /><figcaption>**Figure 3.** *Estimation of the contrast*</figcaption></figure>

OPH 91 Ophir J, Cespedes I, Ponnekanti H, Yazdi Y and Li X, *Elastography : a quantitative method for imaging the elasticity of biological tissues*, Ultrasonic imaging, **13** (1991), 111–134

OPH 02 Ophir J, Alam S, Garra B, Kallel F, Konofagou E, Krouskop T, Merritt C, Righetti R, Souchon R, Srinivasan S and Varghese T, *Elastography : Imaging the Elastic Properties of Soft Tissues with Ultrasound*, J Med Ultrasonics, **29** (2002), 155–171

SOU 04 Souchon R, *Prostate cancer detection and HIFU therapy monitoring using elastography*, PhD thesis, Institut national des sciences appliquées de Lyon, 2004