Inverse problems in elastography and displacement-flow MRI
Carrillo Lincopi, Hugo
DOI:
10.33612/diss.112422123
IMPORTANT NOTE: You are advised to consult the publisher's version (publisher's PDF) if you wish to cite from it. Please check the document version below.
Document Version
Publisher's PDF, also known as Version of record
Publication date: 2020
Link to publication in University of Groningen/UMCG research database
Citation for published version (APA):
Carrillo Lincopi, H. (2020). Inverse problems in elastography and displacement-flow MRI. Rijksuniversiteit Groningen. https://doi.org/10.33612/diss.112422123
Copyright
Other than for strictly personal use, it is not permitted to download or to forward/distribute the text or part of it without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license (like Creative Commons).
Take-down policy
If you believe that this document breaches copyright please contact us providing details, and we will remove access to the work immediately and investigate your claim.
Downloaded from the University of Groningen/UMCG research database (Pure): http://www.rug.nl/research/portal. For technical reasons the number of authors shown on this cover page is limited to 10 maximum.
Introduction
1.1
Motivation
In this thesis we study some inverse problems in MRI and in elasticity imag-ing. A problem is called inverse if its formulation involves the solution of other problem, called direct. In science, the direct problem comes from fundamental laws, such as energy conservation; the inverse problem consists in the determi-nation of the parameters of the equations governing the direct problem, from the knowledge of the solution entirely or a part of it. In other words, by solving an inverse problem, parameter values that cannot be obtained directly can be found.
Magnetic resonance imaging (MRI) is a noninvasive biomedical imaging technique. Its development comes from the study of the nuclear magnetic resonance by Felix Bloch and Edward Purcell [Blo46] and later the development of the tomographic technique by Paul Lauterbur and Peter Mansfield [Lau+73; Man82]. MRI has potential advantages over other techniques: lack of ionizing radiation, free choice of imaging planes, high resolution, capability for tissue characterization, etc. The characterization of tissues can be anatomical and qualitative or even quantitative giving a map of a specific parameter which can be obtained from the MRI technique, such as: proton density, magnetic relaxation times, water diffusivity, velocity and displacement. Each of the parameters mentioned above has the capability of characterizing a biological tissue, in order to detect any disease, see for example [JBP+05; CJ05; MBV06; GME12].
In this thesis we are interested in the recovery of the following parameters from magnetization measurements in MRI: the velocity of blood flowing in vessels and the displacements of tissues under a harmonic regime. Both velocity and displacement imaging is achieved by the so-called Phase-Contrast MRI (PC-MRI), that is, the derivation of the quantity of interest is performed by the substraction (modulo 2π) of two measured phases of the complex transverse
magnetization in MRI.
Imaging the velocity of blood flow finds its main application in the research into cardiovascular diseases, since heart failure is a major cause of death in Western societies and will become even more prevalent with increasing life expectancy [MBV06]. Blood velocity measurements can also lead to the recov-ery of further mechanical information like relative pressures [Yan+96; Ber+18; NB19; Nol+19] and wall shear stress [Sot+12; Sot+16].
Imaging the harmonic displacements in tissues finds its motivation in elas-tography, that is, two-step techniques consisting in measuring the displace-ments of a tissue and then recovering some elastic parameter of it. Examples of elastography techniques are ultrasound and magnetic resonance elastogra-phy (MRE). Several diseases involve changes in the mechanical properties of tissue and normal function of tissue, like skeletal muscle, heart, lungs, gut, etc. The changes in mechanical properties are not best seen precisely in the dis-placement map, but the best contrast is seen in the spatial distribution of any elastic parameter [HSB17]. These elastic parameters can be then reconstructed from the knowledge of the displacement map, using the elasticity equations.
Mathematically, elasticity imaging techniques have given rise to many prob-lems which are challenging in many branches of mathematics. Our focus will be on partial differential equations, studying nonlinear models for elasticity and the measurements. In the literature, inverse problems in elasticity are most studied for linear models, however those models do not describe accurately some phenomena, but nonlinear elasticity models do it.
1.2
Research objectives
The topics presented in this thesis are:
1. The proposal of a new method for the recovery of the velocity of blood from MRI measurements. The proposed method is more robust, meaning less affected by noise and able to measure a wider range of velocities than the usual methods. This method is motivated, proposed and tested in Chapter 2. In addition, in Chapter 3 the natural extension to MRE is discussed.
2. Analysis of nonlinear inverse problems in elasticity, in two cases: as-suming the measurements or the elasticity model being nonlinear. This analysis motivates techniques for the reconstruction of an elasticity pa-rameter. Examples will be presented. This topic is developed in Chapter 4.
In the following two sections we will provide the reader the fundamentals of MRI for the comprehension of velocity and displacement recovery, and also fundamentals of hybrid inverse problems and elasticity equations.
1.3
Magnetic Resonance Imaging
1.3.1
Generalities of MRI
The nuclear magnetic resonance phenomenon is modeled by the Bloch equation of the magnetization M = [M1, M2, M3]| of protons in position X(t) at time
t under the influence of a magnetic field B [Abr61; Sli13]: dM dt = γM × B + 1 T1 (M0− M3)e3− 1 T2 M⊥ (1.1)
where T1and T2are relaxation times, M0 is the magnetization at equilibrium,
M⊥ = M1e1+ M2e2 and (e1, e2, e3) are the basis vectors of the frame of
reference.
The Bloch equation represents the conservation of angular momentum of the spin of the protons, including quantum effects.
By components, this is: dM1 dt = γ(B3M2− B2M3) − M1 T2 , (1.2) dM2 dt = γ(−B3M1− B1M3) − M2 T2 , (1.3) dM3 dt = γ(B2M1− B1M2) − M3− M0 T1 . (1.4)
The main steps in MRI are defined by different stages of the “controlled” magnetic field B:
• Initially, B = B0e3, where B0 has a much larger magnitude compared
with all the other fields described below. The field B0e3 is constant
in time and space and it is present in all the experiment. Under this equilibrium condition, M1= M2= 0.
• The application of a radiofrequency pulse (rf pulse) causes the resonance in spins, that pulse is generated by adding a time-harmonic field B1(t)e1
for a very short period of time. This creates the transversal magnetization to tip towards the transversal plane, i.e. M16= 0, M26= 0.
• When the rf pulse is released, the problem is studied only using equations (1.2)-(1.3), and a non-zero initial condition, because the measurements that the scanner gives us are Mxy= M1+ iM2.
• In order to recover the desired physical quantity, we modify the compo-nent B3. The description for velocity and displacement recovery is given
For our purposes, we will also consider the relaxation time T2 large enough
for neglecting the term M⊥ T2
. Therefore, from equations (1.2)-(1.3), we see that Mxy satisfies the following equation:
dMxy dt = −iγB3Mxy. (1.5) Hence Mxy(X(t), t) = Mxy(t0) exp − iγ Z t t0 B3(X(t0), t0)dt0 .
We can fix Mxy(t0) ∈ R, and then we can write
Mxy(t) = Mxy(t0)eiϕ(X(t),t)
where the phase
ϕ(X(t), t) = −γ Z t
t0
B3(X(t0), t0)dt0 (1.6)
is the key quantity in this work.
1.3.2
Velocity Encoding
We now assume that the spins are moving “locally” following the equation
X(t) = x + tu (1.7)
where x, u are time-constant vectors. In particular, x corresponds to the orig-inal position of the spins at t = 0, and X(t) to the position of the spins at time t. In addition, let us consider the component B3 with the form
B3(X(t), t) = B0+ X(t) · G(t)
where G(t) is a constant-in-space vector. Then, equation (1.6) becomes, at the time of measurement TE, and passing the dependence from X(TE) to x and
u due (1.7): ϕ(X(TE), TE) = ϕ(x, u, TE) = − γB0TE+ γx · m0(G) + γu · m1(G) where m0= Z TE 0 G(t0)dt0, m1= Z TE 0 t0G(t0)dt0
are known as the zeroth and first moments, respectively [KJH18]. It is usual in MRI practice to consider a waveform G(t) such that m0(G) = 0 in order
to remove it from the velocity reconstruction problem. Then the model of the phase ϕ = ϕ(x, u, TE) reads:
ϕ = ϕ0+ γum1(G) , ϕ0:= γB0TE, (1.8)
with v the component of the velocity in direction of G, and G the magnitude of the vector G.
Since ϕ0is not fully known due to inhomogeneities in the magnetic field B0,
we cannot obtain v directly from equation (1.8). So, for each x we consider two measurements ϕG0 and ϕG1 satisfying (1.8) with waveforms G
0 and G1
respectively, so that we can obtain u via the formula
u = ϕ
G0− ϕG1
γ(m1(G0) − m1(G1))
. (1.9)
In practice, the user controls indirectly the strength of the magnetic fields by setting the velocity encoding parameter, VENC, defined by:
VENC(G0, G1) =
π
γ(m1(G0) − m1(G1))
,
which has the dimension of velocity. According to equation (1.9), the velocity can be determined by:
u = ϕ
G0− ϕG1
π VENC(G0, G1). (1.10) This is the standard phase contrast MRI technique for velocity recovery (PC-MRI). Note that u = u(x).
1.3.3
Magnetic Resonance Elastography (MRE)
Instead of assuming that the motion is modelled by (1.7), we adopt the follow-ing model [Man+01; HSB17]:
X(x, t) = x + U (x, t) (1.11) where the displacement U (x, t) is given by
U (x, t) = Reuc(x)eiωt
. (1.12)
This motion is caused by some periodic external force, where uc(x) ∈ Cd
describes a complex steady-state behavior.
The waveform G(t), called motion encoding gradient (MEG) in this context, is zeroth moment nulling. Typically, the MEG G(t) is set either as:
• symmetric, that is G(t) = G(−t), or • antisymmetric, that is G(t) = −G(−t),
where we set t = 0 at the middle of the time of the application of the MEG,that is the MEG is assumed to be applied for t ∈h−T
2, T 2
i
, where T = 2πω. Then the equation (1.6) for the phase of the spins moving under this model and under this gradient field is given by
ϕ(x) = ϕ0+ γ
Z T /2
−T /2
G(t) · X(x, t)dt
where G is the magnitude of the vector G and X is the component of X in direction of G. We will denote by uc the component of uc in the direction of
G. The phase ϕ0 is static and it can be measured by a null MEG, so we are
interested only in the dynamic phase: ϑ(x) = γ Z T /2 −T /2 G(t) · X(x, t)dt We have: ϑ(x) = γ Z T /2 −T /2
G(t)x + Re(uc(x)eiωt)
dt = γx Z T /2 −T /2 G(t)dt + γReuc(x) Z T /2 −T /2 G(t)eiωtdt = γReuc(x) Z T /2 −T /2 G(t)eiωtdt = γRe(uc) Z T2 −T 2 G(t) cos(ωt)dt − Im(uc) Z T2 −T 2 G(t) sin(ωt)dt
where we use the fact that G(t) is zeroth moment nulling, so the term γx
Z T /2
−T /2
G(t)dt in the second line is null. Therefore we obtain:
ϑ(x) ξ(ω, T ) =
Re(uc(x)) for G symmetric, ξ(ω, T ) = γ
Z T2
−T
2
G(t) cos(ωt)dt
Im(uc(x)) for G antisymmetric, ξ(ω, T ) = −γ
Z T2
−T
2
G(t) sin(ωt)dt
where ξ(ω, T ) is known as the encoding efficiency. This is a key controlled parameter in MRE.
Several measurements and the discrete Fourier transform in time To identify displacements at different states of the wave given by (1.12), we consider
Uc(x, τ ) = uc(x)eiωτ
where τ is the time corresponding to the displacement observed by the MRI experiment, which corresponds to the complex displacement at the time of the rf pulse, here uc corresponds to the complex displacement of maximum
amplitude. Then
U (x, τ, t) = ReUc(x, τ )eiωt
is the displacement occurring at time t after the rf pulse. Following the analysis done before, we can obtain by PC-MRI either Re(uc(x)eiωτ) or Im(uc(x)eiωτ)
in the direction of the MEG G(t).
For reconstructing any elastic parameter, it would be better to have a dis-placement field suitable in a time independent equation (see next section). For achieving that, let n = 0, . . . , N − 1, and consider the set of N displacements as follows:
• Considering the harmonic regime as above, let un(x) be the PC-MRI
measurement of the displacement of the spins having its position in x at a certain time τn. That is, we assume the model:
un(x) = Reuc(x)eiωτn if G is symmetric, Imuc(x)eiωτn if G is antisymmetric, (1.13)
where each un is obtained from the measurement ϑn:
un=
ϑn
ξ(ω, T ). (1.14)
• We apply the following discrete Fourier transform for {un}N −1n=0:
uF T = N −1 X n=0 une−i 2πn N and we choose τn = T n N .
• If G(t) is symmetric, then we have: un(x) = Re uc(x)ei 2πn N
and then usF T = N −1 X n=0 Reuc(x)ei 2πn N e−i2πnN = N −1 X n=0 h Re(uc)cos 2πn N −Im(uc) sin 2πn N ih cos2πn N −i sin2πn N i = N −1 X n=0 h Re(uc) cos2 2πn N − Im(uc) sin 2πn N cos2πn N i +ih− Re(uc) cos 2πn N sin2πn N + Im(uc) sin2 2πn N i = N 2uc since N −1 X n=0 cos22πn N = N −1 X n=0 sin22πn N =N 2 and N −1 X n=0 cos2πn N sin2πn N = 0.
• Similarly, we obtain for the case G antisymmetric: uaF T = iN 2 uc Therefore, uc = 2 Nu s F T if G is symmetric, 2i Nu a F T if G is antisymmetric. (1.15)
Now that we have obtain uc is we can recover completely the model (1.12) in
order to characterize the viscoelastic properties of a tissue [HSB17], and also we can use algorithms for recovering the shear modulus from a time independent elasticity equation, see for example [AWZ15].
1.4
Hybrid Inverse Problems in Elasticity
1.4.1
Hybrid Inverse Problems
Hybrid inverse problems are inverse problems that describe coupled-physics phenomena in order to reconstruct a parameter of interest. The main idea of hybrid inverse problems is given in two steps:
• First, a high-resolution inverse problem is solved to provide internal in-formation involving solutions and parameters of a differential equation. • Next, the obtained internal information, also called internal functional,
is used to reconstruct with high-contrast the parameter of interest of the inverse problem.
In some settings, a single modality gives a reconstruction with either high contrast or high resolution, but not both at the same time. Under convenient conditions, the physical coupling defining a hybrid inverse problem result in a reconstruction with a high contrast and resolution.
Hybrid inverse problems are useful in, for example, medical and geophys-ical imaging. Several examples of physgeophys-ical couplings defining hybrid inverse problems are described in [Bal13], for example: the photo-acoustic effect, the ultrasound modulation effect, transient elastography, current density imaging. The mathematical analysis of many hybrid inverse problems falls within the context of linear and nonlinear equations or systems of equations. In particular, in this thesis, the study of the second step is given by the theory of redundant systems of elliptic partial differential equations, coupling system
Lv = f in Ω, Bv = g on ∂Ω, with the functionals related with the measurements
Mv = H in Ω,
where v = (γ, u), γ is the set of parameters to reconstruct, u is the solution of the direct problem, L, B, M are operators acting on v, Ω is an open bounded subset of Rd with smooth boundary.
A main tool in our study of hybrid inverse problems is the symbol of a differential operator, since it will provide important properties of equations. We use the definition of symbol given in [Esk11] and we present it here.
Definition 1.1. If A(x, D) (with D = −i∂x∂ ) is a differential operator having
the form A(x, D)u = m X |k|=0 ak(x) − i ∂ ∂x k u,
then the symbol of A(x, D) is the polynomial A(x, ξ) defined by A(x, ξ) = m X |k|=0 ak(x)ξk
and its principal symbol is the polynomial A0(x, ξ) defined by
A0(x, ξ) =
X
|k|=m
ak(x)ξk.
Note that the definition is related with the following well-known Fourier transform property for partial derivatives
F [Dku](ξ) = ξkF [u](ξ)
where we can see that he symbol of the differential operator Dk is ξk.
1.4.2
Linear Elasticity Equations
We consider elastic deformations that are given by [Amm+15]:
ρUtt− ∇ · S[U ] = F (1.16)
where U (x, t) is the displacement field, F (x, t) is the excitation force density, and S is the Cauchy stress tensor which in case of linear elasticity is given by
S[U ] = 2µ[U ] + λ(tr([U ])Id
where
[U ] = 1
2(∇U + ∇U
|) = ∇SU
is the infinitesimal strain tensor.
As in section 1.3.3, if the force is periodic in time with frequency ω, the resulting displacements will be perodic in time as well. After applying a Fourier transform in time in equation (1.16), we obtain:
ω2u + 2∇ · µ∇Su + ∇(λ∇ · u) = 0 (1.17) where u = u(x, ω) is the Fourier transform of U (x, t) in the variable t, evaluated in the frequency ω.
1.5
Thesis Overview
The main contributions of this thesis are new methods for the reconstruction of the following parameters of interest: the velocity of blood in vessels via magnetic resonance imaging (MRI), the displacements of tissues under a har-monic regime via MRI, and the shear modulus from different models of hybrid imaging.
Inverse problems in MRI.
Imaging velocity fields by phase-contast MRI has limitations, some of them are:
• Since the phase is a number in the interval ] − π, π], a wrong choice of the VENC parameter causes aliasing, that is, there are pixels corresponding to jumps with respect to the true velocity. This is because the PC-MRI velocity satisfies |upc| ≤ |VENC|, so the true velocity has to be less than
|VENC| to avoid aliasing.
• The usual way to reduce the above limitation is to increase the VENC, but this adds noise to the image, since VENC∝VNR (veloticty-to-noise ratio).
In Chapter 2 a method (Optimal Dual-VENC, ODV) for tackling these limita-tions is proposed, in order to augment the range of velocities measured while keeping a low noise level. The approach uses dual-VENC measurements com-bined with appropriate cost function optimization having the form:
J (u) =1 2 2 X j=1 |ei ˆϑGj − eiϑGj(u)|2
where for j = 1, 2, ˆϑGj = ˆϕGj − ˆϕG0 is the “measured” velocity dependent
phase and ϑGj(u) its model depending on the velocity u. The VENCs involved
are defined from the pairs (G0, G1) and (G0, G2), respectively. Experiments
are performed with a phantom and with volunteers in order to verify the theoretical results.
Main Result: The global minimum of J in a “reasonable domain” (explained in Section 2.2.6) gives a new recovered velocity uODV which is not aliased
even if the VENC’s involved are less than the true velocity, keeping a low noise. On the other hand, the limitations described above are also present in imaging harmonic displacements. The adaptation of the method presented in Chapter 2 is developed to deal with harmonic displacements by MRI, and issues related with MRE are discussed. Experiments with a phantom are performed in order to confirm the theoretical results. This is presented in Chapter 3.
Inverse problems in elasticity.
Most of the inverse problems in elasticity have been mathematically studied in the frame of linear elasticity with linear measurements. In Chapter 4 the inverse problem of the recovery of the shear modulus µ is studied from the following models where Ω ⊂ Rd is a bounded domain with smooth boundary, and uj is the displacement field corresponding to the j-th measurement, where
j = 1, . . . , J and J is the number of measurements:
• Linear elasticity in harmonic regime with frequency ω, with elastic energy density measurements: that is, the model reads
ω2u j+ 2∇ · µ∇Suj = −∇pj in Ω, ∇ · uj = 0 in Ω, uj = gj on ∂Ω,
and is coupled to measurements given by: µ
2|∇
Su
j|2= Hj in Ω.
This problem is studied in dimension d = 2 assuming that the pressures pjare known. In addition, a forcing term is added in the model equation,
which is described by a differential operator of order at most 1.
• Saint-Venant model with internal measurements of the displacements, that is, the model becomes:
(Lµ,λ+ Nµ,λ)uj+ ω2uj = 0 in Ω,
uj = g on ∂Ω,
where
Lµ,λuj = 2∇ · µ∇Suj+ ∇(λ∇ · uj),
Nµ,λuj = 2cτ∇ · (µ∇u|j∇uj) + ∇(λ|∇uj|2)
and cτ is a constant in x. The measurements are given by:
uj= Hj in Ω.
This problem is studied in dimension d = 2, 3. In addition, we study the case of linear elasticity with the same model of measurements, which is useful for the analysis of Saint-Venant model.
In each case, the ellipticity of the respective linearized system is studied, ap-plying the curl operator in the models coupled to internal displacements surements. Then, the resulting systems augmented with the respective mea-surements can be written by:
Lv = f in Ω, Bv = g on ∂Ω,
where L, B are differential operators and v = (δµ, {δuj}Jj=1). Necessary
definitions for the development of the chapter are given in Sections 4.2 and 4.3.
Main Results: The following results are presented for each model:
• Ellipticity for the operator L(x, D), with J = 2 sets of measurements, that is, the principal symbol P(x, ξ) is a full rank matrix ∀x ∈ Ω and ∀ξ ∈ Rd\{0}. The ellipticity holds for certain conditions over the infinitesimal strain tensors ∇Suj, where j = 1, 2. This is presented for the different
models in Theorem 4.3, Corollary 4.1, Section 4.6.1 and Section 4.7.1. • Lopatinskii condition for the operator A = (L, B), and then the
exis-tence of a parametrix, wich gives a useful estimate in appropriate Sobolev spaces, showing injectivity of A except for a finite dimensional kernel. This is presented for the different models in Lemma 4.5, Theorem 4.4, Corollary 4.1, Section 4.6.2, and Section 4.7.2.
• A stability estimate which allows to show that the kernel mentioned above is trivial. As a consecuence, it is possible to define the inverse of A. In the case of elastic energy density, it is necessary a lower bound for the frequency ω. This is presented for the different models in Theorem 4.5, Lemma 4.2, Section 4.6.3 and Section 4.7.2.
• The existence of A−1 allows to define an iterative algorithm for the
recov-ery of µ, which is presented and its convergence is shown. See Sections 4.4.5, 4.5.3 and 4.7.3
