Gabriel Peyré

Faouzi Triki

Clair Poignard

Liliana Borcea

Résumé

Les différentes modalités d'imagerie par ondes présentent chacune des limitations en termes de résolution ou de contraste. Dans ce travail, nous modélisons l'imagerie ultrasonore ultrarapide et présentons des méthodes de reconstruction qui améliorent la précision de l'imagerie ultrasonore. Nous introduisons deux méthodes qui permettent d'augmenter le contraste et de mesurer la position superrésolue et la vitesse dans les vaisseaux sanguins. Nous présentons aussi une méthode de reconstruction des paramètres microscopiques en tomographie d'impédance électrique en utilisant des mesures multifréquence et en s'aidant de la théorie de l'homogénisation.

Abstract

Different modalities in wave imaging each present limitations in terms of resolution or contrast. In this work, we present a mathematical model of the ultrafast ultrasound imaging modality and reconstruction methods which can improve contrast and resolution in ultrasonic imaging. We introduce two methods which allow to improve contrast and to locate blood vessels below the diffraction limit while simultaneously estimating the blood velocity. We also present a reconstruction method in electrical impedance tomography which allows reconstruction of microscopic parameters from multi-frequency measurements using the theory of homogenization.

Introduction

In medical imaging, inverse problems are often ill-posed, or limited in their resolution by the physics of the waves at play. There exist several techniques to overcome these difficulties, using additional information about the medium. For example, by reducing the set of admissible solutions and the number of unknowns, by looking for a inclusions with parameters significantly different from those of the surrounding medium [START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF].

Assuming different frequency responses for different tissue components, another approach uses signal separation techniques to reconstruct robust solutions using multi-frequency settings [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF], [START_REF] Giovanni | Disjoint sparsity for signal separation and applications to hybrid inverse problems in medical imaging[END_REF], [START_REF] Giovanni S Alberti | The linearized inverse problem in multifrequency electrical impedance tomography[END_REF].

A third promising technique for improving the robustness of wave-based imaging is to combine different physical types of waves. This allows to alleviate deficiencies of each separate type of waves and to combine their strengths. Example of multi-wave imaging modalities include photo-acoustic and thermo-acoustic imaging [53], magnetic resonance elastography [START_REF] Manduca | Magnetic resonance elastography: non-invasive mapping of tissue elasticity[END_REF], magneto-acousto-electrical tomography [START_REF] Bradley | The role of magnetic forces in biology and medicine[END_REF], magneto-acoustic tomography with magnetic induction [54], and impediography [START_REF] Jossinet | The phenomenology of acousto-electric interaction signals in aqueous solutions of electrolytes[END_REF].

Recently, nanoparticles have been proposed to be used as labels in molecular biology. Plasmon resonant nanoparticles have unique capabilities of enhancing the brightness and directivity of light and confining strong electromagnetic fields [START_REF] Link | Shape and size dependence of radiative, non-radiative and photothermal properties of gold nanocrystals[END_REF]. These nonlinear optical contrast mechanisms reveal new information from biological specimens and tissues.

Finally, one can use the specific dynamics of the imaged elements to improve the robustness of the imaging process. For example in blood flow imaging, the blood dynamics are used to locate blood vessels [START_REF] Bercoff | Ultrafast compound doppler imaging: Providing full blood flow characterization[END_REF]. Such techniques have been successfully used in Dynamic Optical Coherence Tomography [START_REF] Habib Ammari | A signal separation technique for sub-cellular imaging using dynamic optical coherence tomography[END_REF].

This thesis aims at investigating different methods for improving the resolution and contrast of wave-based imaging techniques for imaging of biological tissues. It is focused on two promising non-invasive imaging methods: ultrafast ultrasound, based on the propagation of sound waves, and electrical impedance tomography (EIT), based on the propagation of low frequency electro-magnetic waves. Both ultrasound imaging and EIT present the ad-10 CONTENTS vantage of being fast, relatively cheap and easy to operate. These methods are examples of the tradeoff between contrast and resolution that is encountered in many wave-based imaging techniques. Ultrasound imaging provides a high resolution of the order of less than a millimeter, but since acoustic impedance has a poor contrast in biological media, its contrast is very low. On the contrary, due to its low frequency nature, EIT provides a very low resolution, and due to its ill-posedness is very sensitive to measurements errors. Contrast is not a problem though, since different features in biological tissues present very distinct conductivities.

To overcome these difficulties, several approaches are possible. Ultrafast ultrasound is a promising imaging modality based on acoustic propagation. Instead of using focused waves as is the case in traditional echography, ultrafast ultrasound is based on plane waves produced by an array of piezzoelectrical elements. This allows for very high numbers of images per seconds, up to 10000Hz. This method induces a worse signal-to-noise ratio (SNR) than conventional echography, but by combining the information of a whole sequence of images, this allows for better imaging, for example in blood flow imaging. Signal processing techniques can then be used to improve reconstruction, using knowledge on dynamics of blood and tissue.

We make use of the dynamics of blood flow, and show that it can be used to improve the imaging of blood vessels.

In the case of a cell culture, it is impossible to directly image the microstructure. An idea developed in this thesis to improve the usefulness of EIT is to use the theory of homogenization and a very precise a priori model on the micro-structure. This model, combined with multi-frequency measurements, will allow us to reconstruct precise information about the micro-structure.

Overview of the thesis

The thesis is divided into four chapters.

First chapter

This chapter provides a mathematical analysis of ultrafast ultrasound imaging. This newly emerging modality for biomedical imaging uses plane waves instead of focused waves in order to achieve very high frame rates. We derive the point spread function of the system in the Born approximation for wave propagation and study its properties.

Second chapter

In this chapter, we consider dynamic data for blood flow imaging, and introduce a suitable random model for blood cells. We show that a singular value decomposition method can successfully remove the clutter signal by using the different spatial coherence of tissue and blood signals, thereby providing high-resolution images of blood vessels, even in cases when the clutter and blood speeds are comparable in magnitude. Several numerical simulations are presented to illustrate and validate the approach.

Third chapter

In this chapter, we introduce a signal processing method to produce simultaneous localization and velocity measurements of blood vessels, with superresolution. The method is based on L1 minimization and sparsity, with an added dynamic parameter. Numerical experiments show that this method allow for a reconstruction of both particle location and velocity.

Fourth chapter

In this chapter, we present a simplified electrical model for tissue culture. We derive a mathematical structure for overall electrical properties of the culture and study their dependence on the frequency of the current. We introduce a method for recovering the microscopic properties of the cell culture from the spectral measurements of the effective conductivity. Numerical examples are provided to illustrate the performance of our approach.

Chapter 1

Modeling of Ultrafast Ultrasound Imaging

Introduction

Conventional ultrasound imaging is performed with focused ultrasonic waves [START_REF] Thomas | Diagnostic Ultrasound Imaging: Inside Out[END_REF][START_REF] Kirk | Diagnostic ultrasound: Imaging and blood flow measurements[END_REF]. This yields relatively good spatial resolution, but clearly limits the acquisition time, since the entire specimen has to be scanned. Over the last decade, ultrafast imaging in biomedical ultrasound has been developed [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF][START_REF] Tanter | Ultrafast imaging in biomedical ultrasound[END_REF][START_REF] Demene | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity[END_REF]. Plane waves are used instead of focused waves, thereby limiting the resolution but increasing the frame rate considerably, up to 20,000 frames per second. Ultrafast imaging has been made possible by the recent technological advances in ultrasonic transducers, but the idea of ultrafast ultrasonography dates back to 1977 [START_REF] Bruneel | Ultrafast echotomographic system using optical processing of ultrasonic signals[END_REF]. The advantages given by the very high frame rate are many, and the applications of this new modality range from blood flow imaging [START_REF] Bercoff | Ultrafast compound doppler imaging: providing full blood flow characterization[END_REF][START_REF] Demene | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity[END_REF], deep super-resolution vascular imaging [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF] and functional imaging of the brain [START_REF] Mace | Functional ultrasound imaging of the brain[END_REF][START_REF] Mace | Functional ultrasound imaging of the brain: theory and basic principles[END_REF] to ultrasound elastography [START_REF] Gennisson | Ultrasound elastography: Principles and techniques[END_REF]. In this chapter we focus on blood flow imaging.

A single ultrafast ultrasonic image is obtained as follows [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF]. A pulsed plane wave (focused on the imaging plane -see Figure 1.1b) insonifies the medium, and the back-scattered echoes are measured at the receptor array, a linear array of piezoelectric transducers. These spatio-temporal measurements are then beamformed to obtain a two-dimensional spatial signal. This is what we call static inverse problem, as it involves only a single wave, and the dynamics of the medium is not captured. The above procedure yields very low lateral resolution, i.e. in the direction parallel to the wavefront, because of the absence of focusing. In order to solve this issue, it was proposed to use multiple waves with different angles: these improve the lateral 14 CHAPTER 1. ULTRAFAST ULTRASOUND IMAGING resolution, but has the drawback of reducing the frame rate.

In this work, we provide a detailed mathematical analysis of ultrasound ultrafast imaging. To our knowledge, this is the first mathematical study addressing the important challenges of this emerging and very promising modality. Even though in this work we limit ourselves to formalize the existing methods, the mathematical analysis provided gives important insights, which we expect will lead to improved reconstruction schemes.

First, we carefully study the forward and inverse static problems. In particular, we derive the point spread function (PSF) of the system, in the Born approximation for ultrasonic wave propagation. We investigate the behavior of the PSF, and analyze the advantages of angle compounding. In particular, we study the lateral and vertical resolutions. In addition, this analysis allows us to fully understand the roles of the key parameters of the system, such as the directivity of the array and the settings related to angle compounding.

This chapter is structured as follows. In Section 1.2 we describe the imaging system and the model for wave propagation. In Section 1.3 we discuss the static inverse problem. In particular, we describe the beamforming process, the PSF and the angle compounding technique.

The Forward Problem

The imaging system is composed of a medium contained in R 3 + := {(x, y, z) ∈ R 3 : z > 0} and of a fixed linear array of transducers located on the line z = 0, y = 0. This linear array of piezoelectric transducers (see [START_REF] Thomas | Diagnostic Ultrasound Imaging: Inside Out[END_REF]Chapter 7]) produces an acoustic illumination that is focused in elevation -in the y coordinates, near the plane y = 0 -and has the form of a plane wave in the direction k ∈ S 1 in the x, z coordinates (see Figure 1.1b). Typical sizes for the array length and for the penetration depth are about 10 -1 m.

We make the assumption that the acoustic incident field u i can be approximated as

u i (x, y, z, t) = A z (y) f t -c -1 0 k • (x, z)
, where c 0 is the background speed of sound in the medium. The function A z describes the beam waist in the elevation direction at depth z (between 4 • 10 -3 m and 10 -2 m). This is a simplified expression of the true incoming wave, which is focused by a cylindrical acoustic lens located near the receptor array (see [START_REF] Thomas | Diagnostic Ultrasound Imaging: Inside Out[END_REF]Chapters 6 and 7]). The function f is the waveform describing the shape of the input pulse:

f (t) = e 2πiν 0 t χ (ν 0 t) , χ (u) = e -u 2 τ 2 , (1.1)
Time (s) where ν 0 is the principal frequency and τ the width parameter of the pulse (see Figure 1.1a). Typically, ν 0 will be of the order of 10 6 s -1 . More precisely, realistic quantities are

×10 -6 -1.5 -1 -0.5 0 0.5 1
c 0 = 1.5 • 10 3 m•s -1 , ν 0 = 6 • 10 6 s -1 , τ = 1. (1.2)
Let c : R 3 → R + be the speed of sound and consider the perturbation n given by

n (x) = 1 c 2 (x) - 1 c 2 0 .
We assume that supp n ⊆ R 3 + . The acoustic pressure in the medium satisfies the wave equation

∆u (x, t) - 1 c 2 (x) ∂ 2 ∂t 2 u (x, t) = 0, x ∈ R 3 ,
with a suitable radiation condition on uu i . Let G denote the Green's function for the acoustic wave equation in R 3 [START_REF] Ammari | Mathematical and Statistical Methods for Multistatic Imaging[END_REF][START_REF] Watanabe | Integral transform techniques for Green's function[END_REF]:

G(x, t, x ′ , t ′ ) = - (4π) -1 |x -x ′ | δ (t -t ′ ) -c -1 0 |x -x ′ | .
In the following, we will assume that the Born approximation holds, i.e. we consider only first reflections on scatterers, and neglect subsequent reflections [START_REF] Ammari | Mathematical and Statistical Methods for Multistatic Imaging[END_REF][START_REF] Cho | Waves and fields in inhomogeneous media[END_REF] (in cases when the Born approximation is not valid, nonlinear methods

CHAPTER 1. ULTRAFAST ULTRASOUND IMAGING

have to be used). This is a very common approximation in medical imaging, and is justified by the fact that soft biological tissues are almost acoustic homogeneous, due to the high water concentration. In mathematical terms, it consists in the linearization around the constant sound speed c 0 . In this case, the scattered wave u s := uu i is given by

u s (x, t) = R R 3 n (x ′ ) ∂ 2 u i ∂t 2 (x ′ , t ′ ) G (x, t, x ′ , t ′ ) dx ′ dt ′ , x ∈ R 3 , t ∈ R + ,
since contributions from n∂ 2 t u s are negligible. Therefore, inserting the expressions for the Green's function and for the incident wave yields

u s (x, t) = - R 3 (4π) -1 |x -x ′ | n (x ′ ) A z ′ (y ′ ) f ′′ t -c -1 0 ((x ′ , z ′ ) • k + |x -x ′ |) dx ′ ,
where we set x = (x, y, z) and x ′ = (x ′ , y ′ , z ′ ). Since the waist of the beam in the y direction is small compared to the distance at which we image the medium, we can make the assumption

|x -(x ′ , y ′ , z ′ )| ≃ |x -(x ′ , 0, z ′ )| , x = (x, 0, 0) ∈ R 3 ,
so that the following expression for u s holds for x = (x, 0, 0) ∈ R 3 and t > 0:

u s (x, t) = R 2 -(4π) -1 |x -(x ′ , 0, z ′ )| f ′′ t -c -1 0 ((x ′ , z ′ ) • k + |x -(x ′ , 0, z ′ )|) ñ(x ′ , z ′ )dx ′ dz ′ ,
where ñ is given by

ñ(x ′ , z ′ ) := R n (x ′ ) A z ′ (y ′ ) dy ′ , x ′ = (x ′ , y ′ , z ′ ) ∈ R 3 . (1.3)
Since our measurements are only two-dimensional (one spatial dimension given by the linear array and one temporal dimension), we cannot aim to reconstruct the full three-dimensional refractive index n. However, the above identity provides a natural expression for what can be reconstructed: the vertical averages ñ of n. Since A z is supported near y = 0, ñ reflects the contribution of n only near the imaging plane. In physical terms, ñ contains all the scatterers in the support of A z ; these scatterers are in some sense projected onto y = 0, the imaging plane. For simplicity, with an abuse of notation from now on we shall simply denote ñ by n, since the original threedimensional n will not play any role, due to the dimensionality restriction

THE INVERSE PROBLEM

discussed above. Moreover, for the same reasons, all vectors x and x ′ will be two-dimensional, namely, x = (x, z) and similarly for x ′ . In view of these considerations, for x = (x, 0) ∈ R 2 and t > 0 the scattering wave takes the form

u s (x, t) = - R 2 (4π) -1 |x -x ′ | f ′′ t -c -1 0 (x ′ • k + |x -x ′ |) n (x ′ ) dx ′ . (1.4)
It is useful to parametrize the direction k ∈ S 1 of the incident wave by k = k θ = (sin θ, cos θ) for some θ ∈ R; in practice, |θ| ≤ 0.25 [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF].

The Inverse Problem

The static inverse problem consists in the reconstruction of n (up to a convolution kernel) from the measurements u s at the receptors, assuming that n does not depend on time. This process provides a single image, and will be repeated many times in order to obtain dynamic imaging, as it is discussed in the next sections.

Beamforming

The receptor array is a segment Γ = (-A, A) × {0} for some A > 0. The travel time from the receptor array to a point x = (x, z) and back to a receptor located in u 0 = (u, 0) is given by

τ θ x (u) = c -1 0 (x • k θ + |x -u 0 |) .
The beamforming process [START_REF] Thomas | Diagnostic Ultrasound Imaging: Inside Out[END_REF][START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] consists in averaging the measured signals on Γ at t = τ θ x (u), which results in the image

s θ (x, z) := x+F z x-F z u s u 0 , τ θ x (u) du, x = (x, z) ∈ R 2 + := {(x, z) ∈ R 2 : z > 0}.
The dimensionless aperture parameter F indicates which receptors are chosen to image the location x = (x, z), and depends on the directivity of the ultrasonic array (in practice, 0.25 ≤ F ≤ 0.5 [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF]). In general, F depends on the medium roughness and on θ, but this will not be considered this work. The above identity is the key of the static inverse problem: from the measurements u s ((u, 0), t) we reconstruct s θ (x, z). We now wish to understand how s θ is related to n. In order to do so, observe that by (1.4) we may write for

CHAPTER 1. ULTRAFAST ULTRASOUND IMAGING

x ∈ R 2 + s θ (x, z) = - x ′ ∈R 2 n (x ′ ) x+F z x-F z (4π) -1 |x ′ -u 0 | f ′′ τ θ x (u) -τ θ x ′ (u) du dx ′ = x ′ ∈R 2 g θ (x, x ′ ) n (x ′ ) dx ′ , (1.5)
where g θ is defined as

g θ (x, x ′ ) = - x+F z x-F z (4π) -1 |x ′ -u 0 | f ′′ τ θ x (u) -τ θ x ′ (u) du, (1.6)
(see Figure 1.2a for an illustration in the case when θ = 0). In other words, the reconstruction s θ is the result of an integral operator given by the kernel g θ applied to the refractive index n. Thus, the next step is the study of the point spread function (PSF) g θ (x, x ′ ), which should be thought of as the image corresponding to a delta scatterer in x ′ .

The point spread function

In its exact form, it does not seem possible to simplify the expression for g further: we will have to perform some approximations. First, observe that

setting h θ x,x ′ (u) = τ θ x (u) -τ θ x ′ (u) for x, x ′ ∈ R 2 + we readily derive (h θ x,x ′ ) ′ (u) = c -1 0 ( u -x |x -u 0 | - u -x ′ |x ′ -u 0 | ) ≈ c -1 0 ( u -x |x ′ -u 0 | - u -x ′ |x ′ -u 0 | ) = c -1 0 x ′ -x |x ′ -u 0 | ,
for x close to x ′ (note that, otherwise, the magnitude of the PSF would be substantially lower). As a consequence, by (1.6) we have

g θ (x, x ′ ) ≈ c 0 (4π) -1 x -x ′ x+F z x-F z (h θ x,x ′ ) ′ (u)f ′′ h θ x,x ′ (u) du = c 0 (4π) -1 x -x ′ f ′ (h θ x,x ′ (x + F z)) -f ′ (h θ x,x ′ (x -F z)) .
(1.7)

In order to simplify this expression even further, let us do a Taylor expansion of

w θ ± (x, z) := h θ x,x ′ (x ± F z) with respect to (x, z) around (x ′ , z ′ ). Direct calculations show that w θ ± (x ′ , z ′ ) = 0, ∇w θ ± (x ′ , z ′ ) = c -1 0 C F (C F sin θ ∓ F, 1 + C F cos θ),
where we define

C F := √ 1 + F 2 . Whence h θ x,x ′ (x ± F z) ≈ c -1 0 C F ((1 + C F cos θ)(z -z ′ ) +(C F sin θ ∓ F )(x -x ′ )).
Substituting this expression into (1.7) yields

g θ (x, x ′ ) ≈ gθ (x -x ′ ), (1.8)
where (see Figure 1.2b for an illustration in the case θ = 0), thereby allowing to write the image s θ given in (1.5) as a convolution of gθ and the refractive index n, namely

gθ (x) = c 0 4πx f ′ c -1 0 C F ((1 + C F cos θ)z + (C F sin θ -F )x) -f ′ c -1 0 C F ((1 + C F cos θ)z + (C F sin θ + F )x) , (1.9)
s θ (x) = x ′ ∈R 2 gθ (x -x ′ )n (x ′ ) dx ′ = (g θ * n)(x), x ∈ R 2 + .
The validity of this approximation, obtained by truncating the Taylor expansion of w θ ± at the first order, is by no means obvious. Indeed, by construction, the pulse f (t) is highly oscillating (ν 0 ≈ 6 • 10 6 s -1 ), and therefore even small variations in t may result in substantial changes in f (t). However, this does not happen, since if (x, z) is not very close to (x ′ , z ′ ) then the magnitude of the PSF is very small, if compared to the maximum value. The verification of this fact is quite technical, and thus is omitted: the details may be found in Appendix A.

Remark 1. From this expression, it is easy to understand the role of the aperture parameter F , which depends on the directivity of the array. Ignoring the second order effect in F and taking, for simplicity θ = 0, we can further simplify the above expression as

g0 (x) ≈ c 0 4πx f ′ c -1 0 (2z -F x) -f ′ c -1 0 (2z + F x) .
It is clear that F affects the resolution in the variable x: the higher F is, the higher the resolution is. Moreover, the aperture parameter affects also the orientation of the diagonal tails in the PSF. These two phenomena can be clearly seen in Figure 1.3. In general, the higher the aperture is the better for the reconstruction: as expected, the intrinsic properties of the array affects the reconstruction. Remark 2. It is also easy to understand the role of the angle θ. In view of

gθ (x) ≈ c 0 4πx f ′ c -1 0 ((1 + cos θ)z + (sin θ -F )x) -f ′ c -1 0 ((1 + cos θ)z + (sin θ + F )x) ,
an angle θ = 0 substantially gives a rotation of the PSF; see Figure 1.4.

We have now expressed g θ as a convolution kernel. In order to better understand the different roles of the variables x and z, it is instructive to use the actual expression for f given in (1.1). Since f ′ (t) = ν 0 e 2πiν 0 t χ(ν 0 t), with χ(t) = 2πiχ(t) + χ ′ (t), we can write

f ′ c -1 0 C F ((1 + C F cos θ)z + (C F sin θ ± F )x) = ν 0 e 2πiν 0 c -1 0 C F ((1+C F cos θ)z+(C F sin θ±F )x) χ ν 0 c -1 0 C F ((1 + C F cos θ)z + (C F sin θ ± F )x) ≈ ν 0 e 2πiν 0 c -1 0 (2z+(θ±F )x) χ 2ν 0 c -1 0 z ,
where we have approximated the dependence on F and θ at first order around F = 0 and θ = 0 in the complex exponential (recall that F and θ are small) and at zero-th order (F = 0 and θ = 0) inside χ: the difference in the orders is motivated by the fact that the variations of the complex exponentials have much higher frequencies than those of χ, since several oscillations are contained in the envelope defined by χ, as it can be easily seen in Figure 1.1a (and similarly for χ ′ ). This approximation may be justified by arguing as in Appendix A. Inserting this expression into (1.9) yields

gθ (x) ≈ c 0 ν 0 4πx e 2πiν 0 c -1 0 (2z+(θ-F )x) χ 2ν 0 c -1 0 z -e 2πiν 0 c -1 0 (2z+(θ+F )x) χ 2ν 0 c -1 0 z = - iν 0 c 0 2πx χ 2ν 0 c -1 0 z e 4πiν 0 c -1 0 z e 2πiν 0 c -1 0 θx sin(2πν 0 c -1 0 F x),
whence for every

x = (x, z) ∈ R 2 gθ (x) ≈ -iν 2 0 F χ 2ν 0 c -1 0 z e 4πiν 0 c -1 0 z e 2πiν 0 c -1 0 θx sinc(2πν 0 c -1 0 F x), (1.10)
where sinc(x) := sin(x)/x (see Figure 1.2c). This final expression allows us to analyze the PSF gθ , and in particular its different behaviors with respect to the variables x and z. Consider for simplicity the case θ = 0 (with τ = 1).

In view of the term χ 2ν 0 c -1 0 z , the vertical resolution is approximately 0.8 • ν -1 0 c 0 ; similarly, in view of the term sinc(2πν 0 c -1 0 F x), the horizontal resolution is approximately 1 2F ν -1 0 c 0 . Even though horizontal and vertical resolutions are comparable, in terms of focusing and frequencies of oscillations the PSF has very different behaviours in the two directions. Indeed, we can observe that the focusing in the variable z is sharper than that in the variable x: the decay of χ is much stronger than the decay of sinc. Moreover, in the variable z we have only high oscillations, while in the variable x the highest oscillations are at least four times slower (2 = 4 1 2 ≥ 4F ), and very low frequencies are present as well, due to the presence of the sinc. As it is clear from Figure 1.2, this approximation introduces evident distortions of the tails, as it is expected from the approximation F = 0 inside χ; however, the center of the PSF is well approximated. Similar considerations are valid for the case when θ = 0: as observed before, this simply gives a rotation.

The same analysis may be carried out by looking at the expression of the PSF in the frequency domain. For simplicity, consider the case θ = 0: the general case simply involves a translation in the frequency domain with respect to x. Thanks to the separable form of gθ given in (1.10), the Fourier transform may be directly calculated, and results in the product of the Fourier transform of χ and the Fourier transform of the sinc. More precisely, we readily derive 0 : the PSF is a low pass filter with cut-off frequency F ν 0 c -1 0 with respect to the variable x and a band pass filter around 2ν 0 c -1 0 with respect to z.

F gθ (ξ x ,ξ z ) = R 2 gθ (x, z)e -2πi(xξx+zξz) dxdz ≈ -iν 2 0 F R sinc(2πν 0 c -1 0 F x)e -2πixξx dx × R χ ν 0 c -1 0 z e -2πi(-2ν 0 c -1 0 +ξz)z dz.
Thus, since the Fourier transform of the sinc may be easily computed and is a suitable scaled version of the rectangle function, we have

F gθ (ξ x , ξ z ) ≈ -iν 2 0 F 1 2ν 0 c -1 0 F 1 [-F,F ] c 0 ν -1 0 ξ x × R χ ν 0 c -1 0 z e -2πi(-2ν 0 c -1 0 +ξz)z dz = - ic 0 ν 0 2 1 [-F,F ] c 0 ν -1 0 ξ x 1 ν 0 c -1 0 F χ -2ν 0 c -1 0 + ξ z ν 0 c -1 0 , whence F gθ (ξ x , ξ z ) ≈ -ic 2 0 1 [-F,F ] c 0 ν -1 0 ξ x F χ -2 + ν -1 0 c 0 ξ z /2. (1.11)
Therefore, up to a constant, the Fourier transform of the PSF is a lowpass filter in the variable x with cut-off frequency F ν 0 c -1 0 and a band pass filter in z around 2ν 0 c -1 0 (since χ is a low-pass filter). This explains, from another point of view, the different behaviors of gθ with respect to x and z. This difference is evident from Figure 1.5, where the absolute values of the Fourier transforms of the different approximations of the PSF are shown.

Angle compounding

We saw in the previous subsection that, while very focused in the direction z, the PSF is not very focused in the direction x due to the presence of the sinc function, see (1.10). In order to have a better focusing, it was proposed in [START_REF] Montaldo | Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[END_REF] to use multiple measurements corresponding to many angles in an interval θ ∈ [-Θ, Θ] for some 0 ≤ Θ ≤ 0.25. The reason why this technique is promising is evident from Figure 1.4: adding up several angles together will result in an enhancement of the center of the PSF, and in a substantial reduction of the artifacts caused by the tails in the direction x. Let us now analyze this phenomenon analytically.

In a continuous setting, angle compounding corresponds to setting

s ac Θ (x) = 1 2Θ Θ -Θ s θ (x) dθ, x ∈ R 2 + .
(1.12) Thus, by linearity, the corresponding PSF is given by

g ac Θ (x, x ′ ) = 1 2Θ Θ -Θ g θ (x, x ′ ) dθ, x, x ′ ∈ R 2 + . (1.13)
Let us find a simple expression for g ac Θ . By using (1.8), we may write

g ac Θ (x, x ′ ) ≈ gac Θ (x -x ′ ), where gac Θ is given by gac Θ (x) = 1 2Θ Θ -Θ gθ (x)
dθ, so that the image may be expressed as

s ac Θ (x) = (g ac Θ * n)(x), x ∈ R 2 + . (1.14)
Thus, in view of the approximation (1.10), we can write

gac Θ (x) = - iν 2 0 F 2Θ Θ -Θ χ 2ν 0 c -1 0 z e 4πiν 0 c -1 0 z e 2πiν 0 c -1 0 θx sinc(2πν 0 c -1 0 F x) dθ = -iν 2 0 F χ 2ν 0 c -1 0 z e 2iν 0 c -1 0 z sinc(2πν 0 c -1 0 F x)sinc(2πν 0 c -1 0 Θx). Therefore, we immediately obtain gac Θ (x) = g0 (x)sinc(2πν 0 c -1 0 Θx), x ∈ R 2 . (1.15)
This expression shows that the PSF related to angle compounding is nothing else than the PSF related to the single angle imaging with θ = 0 multiplied by sinc(2πν 0 c -1 0 Θx). Thus, for Θ = 0 we recover gθ for θ = 0, as expected. However, for Θ > 0, this PSF enjoys faster decay in the variable x. See Figure 1.6 for an illustration of g ac Θ and gac Θ and a comparison with g θ and Figure 1.5d for an illustration of the Fourier transform of g ac Θ . To sum up the main features of the static problem, we have shown that the recovered image may be written as s ac Θ = gac Θ * n, where gac Θ is the PSF of the imaging system with measurements taken at multiple angles. The ultrafast imaging technique is based on obtaining many of these images over time, as we discuss in the next section. 2), and F = 0.4). The better focusing in the variable x for g ac Θ is evident, as well as the good approximation given by gac Θ . The size of the square shown is 2 mm × 2 mm, and the horizontal and vertical axes are the x and z axes, respectively.

Conclusion

In this chapter, we have provided for the first time a detailed mathematical analysis of ultrafast ultrasound imaging. We have derived an approximate expression for the PSF of this imaging system, which allows for a detailed analysis of blood flow imaging in the next chapter.

Chapter 2 Blood Flow Imaging in Ultrafast Ultrasound

Introduction

In this chapter, we consider dynamic imaging, that is the ultrafast ultrasound imaging process described in the previous chapter is repeated many times, which gives several thousand images per second. In blood flow imaging, we are interested in locating blood vessels. One of the main issues lies in the removal of the clutter signal, typically the signal scattered from tissues, as it introduces major artifacts [START_REF] Bjaerum | Clutter filter design for ultrasound color flow imaging[END_REF]. Ultrafast ultrasonography allows to overcome this issue, thanks to the very high frame rate. Temporal filters [START_REF] Bercoff | Ultrafast compound doppler imaging: providing full blood flow characterization[END_REF][START_REF] Mace | Functional ultrasound imaging of the brain[END_REF][START_REF] Mace | Functional ultrasound imaging of the brain: theory and basic principles[END_REF], based on high-pass filtering the data to remove clutter signals, have shown limited success in cases when the clutter and blood velocities are close (typically of the order of 10 -2 m•s -1 ), or even if the blood velocity is smaller than the clutter velocity. A spatio-temporal method based on the singular value decomposition (SVD) of the data was proposed in [START_REF] Demene | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity[END_REF] to overcome this drawback, by exploiting the different spatial coherence of clutter and blood scatterers. Spatial coherence is understood as similar movement, in direction and speed, in large parts of the imaged zone. Tissue behaves with higher spatial coherence when compared to the blood flow, since large parts of the medium typically move in the same way, while blood flow is concentrated only in small vessels, which do not share necessarily the same movement direction and speed. This explains why spatial properties are crucial to perform the separation.

The analysis of the PSF provided in the previous chapter allows to study the Doppler effect, describing the dependence on the direction of the flow. Moreover, we consider a random model for the movement of blood cells,

CHAPTER 2. BLOOD FLOW IMAGING

which allows us to study and justify the SVD method for the separation of the blood signal from the clutter signal, leading to the reconstruction of the blood vessels' geometry. The analysis is based on the empirical study of the distribution of the singular values, which follows from the statistical properties of the relative data. We provide extensive numerical simulations, which illustrate and validate this approach.

This chapter is structured as follows. In Section 2.2 the dynamic forward problem is considered: we briefly discuss how the dynamic data are obtained and analyze the Doppler effect. In Section 2.3 we focus on the source separation to solve the dynamic inverse problem. We discuss the random model for the refractive index and the method based on the SVD decomposition of the data. In Section 2.4 numerical experiments are provided. Some concluding remarks and outlooks are presented in the final section.

The Forward Problem

The quasi-static approximation and the construction of the data

The dynamic imaging setup consists in the repetition of the static imaging method over time to acquire a collection of images of a medium in motion.

We consider a quasi-static model: the whole process of obtaining one image, using the image compounding technique discussed in Subsection 1.3.3, is fast enough to consider the medium static, but collecting several images over time gives us a movie of the movement over time. In other words, there are two time scales: the fast one related to the propagation of the wave is considered instantaneous with respect to the slow one, related to the sequence of the images.

In view of this quasi-static approximation, from now on we neglect the time of the propagation of a single wave to obtain static imaging. The time t considered here is related to the slow time scale. In other words, by (1.14) at fixed time t we obtain a static image s(x, t) of the medium n = n(x, t), namely

s(x, t) = (g ac Θ * n( • , t)) (x). (2.1)
Repeating the process for t ∈ [0, T ] we obtain the movie s(x, t), which represents the main data we now need to process. As mentioned in the introduction, our aim is locating the blood vessels within the imaged area, by using the fact that s(x, t) will be strongly influenced by movements in n.

THE FORWARD PROBLEM

The Doppler effect

Measuring the medium speed is an available criterion to separate different sources; thus, we want to see the influence on the image of a single particle in movement, as by linearity the obtained conclusions naturally extend to a group of particles. For a single particle, we are interested in observing the generated Doppler effect in the reconstructed image, namely peaks in the Fourier transform away from zero.

Intuitively, Figure 1.5d shows that there is a clear difference in the movements depending on their orientation. We want to explore this difference in a more precise way. Let us consider n(x, z, t) = δ (0,vt) (x, z), i.e. a single particle moving in the z direction with velocity v. The resulting image, as a function of time, is obtained via equations (1.15) and (2.1)

s(x, z, t) ≈ R 2 gac Θ (x -x ′ , z -z ′ )δ (0,vt) (x ′ , z ′ )dx ′ dz ′ = gac Θ (x, z -vt) = g0 (x, z -vt)sinc(2πν 0 c -1 0 Θx).
Therefore, arguing as in (1.11), we obtain that the Fourier transform with respect to the time variable t of the image is given by

F t (s)(x, z, ξ) ≈ R g0 (x, z -vt)e -2πiξt dt sinc(2πν 0 c -1 0 Θx) = 1 v e -2πi ξz v F 2 (g 0 )(x, - ξ v )sinc(2πν 0 c -1 0 Θx),
where F 2 is the Fourier transform with respect to the variable z. Adopting approximation (1.10), we obtain

F t (s)(x, z, ξ) ≈ - 1 v iν 2 0 F e -2πi ξz v sinc(2πν 0 c -1 0 Θx) × sinc(2πν 0 c -1 0 F x)F( χ) -ξ 2ν 0 c -1 0 v -1 .
Given the shape of χ, its Fourier transform has a maximum around 0, thus we can see a peak of |F t (s)(x, z, ξ)| when ξ is around -2ν 0 c -1 0 v, and so we have the Doppler effect.

In the case when the particle is moving parallel to the detector array, namely n(x, z, t) = δ (vt,0) (x, z), following an analogous procedure as before, we obtain s(x, z, t) ≈ g0 (xvt, z)sinc(2πν 0 c -1 0 Θ(xvt)),

CHAPTER 2. BLOOD FLOW IMAGING

and applying the Fourier transform in time yields

F t (s)(x, z, ξ) ≈ 1 v e -2πi ξx v F(g 0 (•, z)sinc(2πν 0 c -1 0 Θ•))(- ξ v ).
Using approximation (1.10), the convolution formula for the Fourier transform and the known transform of the sinc function, gives

F t (s)(x, z, ξ) ≈ -i e -2πi ξx v 4Θv ν 0 c 0 χ(2ν 0 c -1 0 z)e 4πiν 0 c -1 0 z × (1 [-F,F ] * 1 [-Θ,Θ] ) - ξ vν 0 c -1 0 .
The convolution of these characteristic functions evaluated at η is equal to the length of interval

[-F + η, F + η] ∩ [-Θ, Θ], because (1 [-F,F ] * 1 [-Θ,Θ] )(η) = R 1 [-F,F ] (η -s)1 [-Θ,Θ] (s)ds = R 1 [-F +η,F +η] (s)1 [-Θ,Θ] (s)ds.
Since both intervals are centered at 0, this value is maximized for η (and thus ξ) around 0, like in the static case, and so the observed Doppler effect is very small. These differences are fundamental to understand the capabilities of the method for blood flow imaging. This phenomenon will be experimentally verified in Section 2.4.

Multiple scatterer random model

We have seen the effect on the image s(x, z, t) of a single moving particle. We now consider the more realistic case of a medium (either blood vessels or tissue) with a large number of particles in motion. This will allow to study the statistical properties of the resulting measurements.

We consider a rectangular domain Ω = (-L x /2, L x /2) × (0, L z ), which consists in N point particles. Let us denote the location of particle k at time t by a k (t). In the most general case, each particle is subject to a dynamics

a k (t) = ϕ k (u k , t) , a k (0) = u k , (2.2)
where (u k ) k=1,...,N are independent uniform random variables on Ω and (ϕ k ) k=1,...,N are independent and identically distributed stochastic flows: for instance, they can be the flows of a stochastic differential equation or the deterministic flows of a partial differential equation. Thus, the a k s are independent and identically distributed stochastic processes. In view of these considerations, we consider the medium given by

n (x, t) = C √ N N k=1 δ a k (t) (x) , (2.3)
where C > 0 denotes the scattering intensity and 1

√

N is the natural normalization factor in view of the central limit theorem.

To avoid minor issues from boundary effects, which are of no interest to us in the analysis of this problem, we assume the periodicity of the medium. In other words, we consider the periodization

n p (x, t) = l∈Z 2 n(x + l • L, t), (2.4)
where L = (L x , L z ). Let g (x) := l∈Z 2 gac Θ (x + l • L) be the periodic PSF, which is more convenient than gac Θ (given by (1.15)) for a Ω-periodic medium. The dynamic image s is then given by

s(x, t) = (g ac Θ * n p (•, t)) (x) = (g * n( • , t))(x) = C √ N N k=1 g (x -a k (t)) .
Let us also assume for the sake of simplicity that, at every time t, a k (t) modulo Ω is a uniform random variable on Ω, namely

E l∈Z 2 w(a k (t) + l • L) = |Ω| -1 R 2 w(y) dy, w ∈ L 1 (R 2 ). (2.5)
As a simple but quite general example, it is worth noting that in the case when a k (t) = u k + F (t), where F (t) is any random process independent of u k , the above equality is satisfied, since

E l∈Z 2 w(u k + F (t) + l • L) = |Ω| -1 E l∈Z 2 Ω w(y + F (t) + l • L)dy = |Ω| -1 R 2 w(y) dy,
where the expectation in the first term is taken with respect to u k and F (t), while in the second term only with respect to F (t).

We now wish to compute the expectation of the random variables present in the expression for s(x, t). By (1.10)

E (g (x -a k (t))) = E l∈Z 2 gac Θ (x-a k (t)+l•L) = |Ω| -1 R 2 gac Θ (y)dy = 0. (2.6)
Let (x i ) i=1,...,mx and (t j ) j=1,...,mt be the sampling locations and times respectively. The data may be collected in the Casorati matrix S N ∈ C mx×mt defined by S N (i, j) = s(x i , t j ).

By (2.6), according to the multivariate central limit theorem, the matrix S N converges in distribution to a Gaussian complex matrix S ∈ C mx×mt , the distribution of which is entirely determined by the following correlations, for i, i ′ = 1, . . . , m x and j, j

′ = 1, . . . , m t E(S(i, j)) = 0, Cov (S(i, j), S(i ′ , j ′ )) = C 2 E (g (x i -a 1 (t j )) g (x i ′ -a 1 (t j ′ ))) , (2.7)
Cov S(i, j), S(i

′ , j ′ ) = C 2 E g (x i -a 1 (t j )) g (x i ′ -a 1 (t j ′ )) . (2.8)
More precisely, let w ∈ C mxmt be a column vector containing all the entries of S. Let v ∈ C 2mxmt and V ∈ C 2mxmt×2mxmt be defined by v = (w 1 , w 1 , w 2 , w 2 , ..., w mxmt , w mxmt ) T and V = E vv T .

The covariance matrix V can be easily computed from (2.7) and (2.8). Then the probability density function f of v can be expressed as [START_REF] Van Den Bos | The multivariate complex normal distribution-a generalization[END_REF]:

f (v) = 1 π mxmt det (V ) 1 2 exp - 1 2 v * V -1 v .
Moreover, it is possible to generate samples from this distribution: if X is a complex unit variance independent normal random vector, and if

√ V is a square root of V , then √ V X is distributed like v.
This allows for simulations of sample image sequences for a large number of particles with a complexity independent of the number of particles.

The analysis carried out here will allow us to study the distribution of the singular value of the matrix S, depending on the properties of the flows ϕ k . This will be the key ingredient to justify the correct separation of blood and clutter signals by means of the singular value decomposition of the measurements. As explained in the introduction, the aim of the dynamic inverse problem is blood flow imaging. In other words, we are interested in locating blood vessels, possibly very small, within the medium. The main issue is that the signal s(x, t) is highly corrupted by clutter signal, namely the signal scattered from tissues. In the linearized regime we consider, we may write the refractive index n as the sum of a clutter component n c and a blood component n b , namely n = n c + n b . Blood is located only in small vessels in the medium, whereas clutter signal comes from everywhere: by (1.3), since blood vessels are smaller than the focusing height, even pixels located in blood vessels contain reflections coming from the tissue. Let us denote the location of blood vessels by Ω b ⊂ Ω. The inverse problem is the following: can we recover Ω b from the data s(x, t) = s c (x, t) + s b (x, t)? Here, s c and s b are given by (2.1), with n replaced by n c and n b , respectively. In this section, we provide a quantitative analysis of the method described in [START_REF] Demene | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity[END_REF] based on the singular value decomposition (SVD) of s.

THE INVERSE PROBLEM: SOURCE SEPARATION

The SVD algorithm

We now review the SVD algorithm presented in [START_REF] Demene | Spatiotemporal clutter filtering of ultrafast ultrasound data highly increases doppler and fultrasound sensitivity[END_REF]. The Casorati matrix S ∈ C mx×mt is defined as in previous section by

S(i, j) = s (x i , t j ) , i ∈ {1, ..., m x } , j ∈ {1, ..., m t } .
Without loss of generality, we further assume that m t ≤ m x . We remind the reader that the SVD of S is given by

S = mt k=1 σ k u k v k T ,
where (u 1 , ..., u mx ) and (v 1 , ..., v mt ) are orthonormal bases of C mx and C mt , and

σ 1 ≥ σ 2 ≥ ... ≥ σ mt ≥ 0. For any K ≥ 1, S K = K k=1 σ k u k v k T
is the best rank K approximation of S in the Frobenius norm. The SVD is a wellknown tool for denoising sequences of images, see for example [START_REF] Iizuka | Data adaptive signal estimation by singular value decomposition of data matrix[END_REF]. The idea is that since singular values for the clutter signal are quickly decaying after a certain threshold, the best rank K approximation of S will contain most of the signal coming from the clutter, provided that K is large enough. This could be used to recover clutter data, by applying a "denoising" algorithm, and keeping only S K . But it can also be used to recover the blood location, by considering the "power Doppler"

Ŝb,K (i) := mt k=K+1 σ 2 k |u k | 2 (i) = mt j=1 |(S -S K ) (i, j)| 2 , i ∈ {1, ..., m x } .
As we will show in the following subsection, clutter signal can be well approximated by a low-rank matrix. Therefore, S K will contain most of the clutter signal for K large enough. In this case, even if the intensity of total blood reflection is small, S -S K will contain more signal coming from the blood than from the clutter and therefore high values of Ŝb,K (i) should be located in blood vessels. Before presenting the justification of this method, let us briefly provide a heuristic motivation by considering the SVD of the continuous data given by

s(x, t) = ∞ k=1 σ k u k (x)v k (t).
In other words, the dynamic data s is expressed as a sum of spatial components u k moving with time profiles v k , with weights σ k . Therefore, since the tissue movement has higher spatial coherence than the blood flow, we expect the first factors to contain the clutter signal, and the remainder to provide information about the blood location via the quantity Ŝb,K .

Justification of the SVD in 1D

We will assume that the particles of the blood and of the clutter have independent dynamics described by (2.2)-(2.4). We add the subscripts b and c to indicate the dynamics of blood and clutter, respectively.

In this subsection, using the limit Gaussian model presented in §2.2.3, we present the statistics of the singular values in a simple 1D model. These are useful to understand the behavior of SVD filtering. The results of §2.2.3 allow to simulate large number of sample signals s, given that we can compute the covariance matrices (2.7) and (2.8). Since these matrices are very large, we restrict ourselves to the 1D case, so that all sampling locations x i are located at x = 0, and are thus characterized by their depth z i . We will therefore drop all references to x in the following. We also consider very simplified dynamics, which can be thought of as local descriptions of the global dynamics at work in the medium. Let a b = a 1,b and a c = a 1,c be the random variables for the dynamics of blood and clutter particles, respectively, as introduced in (2.2). The dynamics is modeled by a Brownian motion with drift, namely

a α (t) = u α + v α t + σ α B t , α ∈ {b, c} .

THE INVERSE PROBLEM: SOURCE SEPARATION

Here, u α represents the position of the particle at time t = 0, and is uniformly distributed in (0, L α ), where L b ≪ L c . The deterministic quantity v α is the mean velocity of the particles. In order to take into account the random fluctuations of the particles in movement, we added a diffusion term σ α B t , where B t is a Brownian motion and σ 2 α is a diffusion coefficient quantifying the variance of the fluctuations of the particle position relative to the mean trajectory. We also make the simplifying assumption that the diffusion terms are independent over different particles. More precisely, we have the following conditional expectation and variance:

E ( a α (t)| u α ) = u α + v α t, Var ( a α (t)| u α ) = tσ 2 α .
The difference between clutter and blood dynamics is in the diffusion coefficient: in the case of clutter, since it is an elastic displacement, σ 2 c ≈ 0. For simplicity, from now on we set σ c = 0. In the case of blood, which is modeled as a suspension of cells in a fluid, we have σ 2 b = σ 2 > 0. This coefficient is expressed in m 2 s -1 , and models the random diffusion in a fluid transporting red blood cells due to turbulence in the fluid dynamics and collisions between cells. In practice, σ 2 is much larger than the diffusion coefficient of microscopic particles in a static fluid, and depends on the velocity v b [START_REF] Gerald | The mechanics of the circulation[END_REF]. As for the mean velocities, in the most extreme cases, v b and v c can be of the same order, even though most of the time v b > v c .

Let S b and S c denote the data matrix constructed in §2.2.3, related to blood and clutter signal, respectively. We now compute the covariance matrix V of S α :

Cov(S α (i, j), S α (i ′ , j ′ )) = C 2 α E (g (z i -a α (t j )) g (z i ′ -a α (t j ′ ))) = C 2 α L E L 0 g z i -y -v α t j -σ α v α B t j × g z i ′ -y -v α t j ′ -σ α v α B t j ′ dy = C 2 α E C gg z i -z i ′ + v α t j ′ -t j + σ α (B t j ′ -B t j ) ,
where C gg (z) = 1 L L 0 g(y)g(z + y)dy and C b and C c denote the intensity of the blood and clutter signals, respectively. The expectation operator is taken over all possible positions u α and all possible drifts B t j and B t j ′ in the first line, and only over all drifts in the second and third lines. By standard properties of the Brownian motion, B t j ′ -B t j is Gaussian distributed, of expected value 0 and variance |t jt j ′ | and so it has the same distribution as B t j ′ -t j . Thus, in the case of the blood, we can write

Cov (S b (i, j), S b (i ′ , j ′ )) = C 2 b E C gg z i -z i ′ + v b (t j ′ -t j + σ b B t j ′ -t j ) . Likewise, Cov S b (i, j), S b (i ′ , j ′ ) = C 2 b E C gḡ z i -z i ′ + v b (t j ′ -t j + σ b B t j ′ -t j ) ,
where C gḡ (z) = 1 L L 0 g(y)ḡ(z + y)dy. The tissue model is then given by σ c = 0, and is therefore deterministic given the initial position. Thus

Cov (S c (i, j), S c (i ′ , j ′ )) = C 2 c C gg (z i -z i ′ + v c (t j ′ -t j )) , Cov S c (i, j), S c (i ′ , j ′ ) = C 2 c C gḡ (z i -z i ′ + v c (t j ′ -t j )) .
On one hand, in the case of blood, since C gḡ and C gg are oscillating and with very small support (see Figures 2.1a and 2.1b), the integration done when taking the expectation in the blood case should yield small correlations as long as |t j ′t j | is large enough. On the other hand, in the case of clutter, correlations will be high between the two signals as long as z iz i ′ and v c (t jt j ′ ) are of the same order and almost cancel out. This heuristic is confirmed by numerical experiments. In Figure 2.1c, we compare the clutter model and the blood model in one dimension: velocities are in the z direction, and we only consider points aligned on the z axis. As we can see, correlations are quickly decaying as we move away from (0, 0) in the case of blood. In the case of clutter, there are correlations at any times at the corresponding displaced locations.

Once the correlation matrix is computed, we can generate a large number of samples to study the distribution of the singular values in different cases. In Figure 2.2a, we compare the distribution in the two models (blood and clutter), using the Gaussian limit approximation for the simulations, with the same intensity for both models. A comparison with a white noise model with the same variance shows that blood and noise have approximately the same singular value distribution. On the contrary, the distribution of the singular values of clutter presents a much larger tail. A comparison of the distribution of the singular values for the clutter model at different velocities shows no real difference in the tail of the distribution (Figure 2.2b).

As a consequence, the clutter signal s c is well approximated by a low rank matrix, and the blood signal can be thought of as if it were only noise. Therefore, the SVD method act as a denoising algorithm and extracts the clutter signal, according to the discussion in the previous subsection.

Numerical Experiments

In this section, we consider again a more realistic 2D model, given by (2.2). This framework will allow us to simulate generic blood flow imaging sequences from particles. The dynamics of blood and clutter are modeled as follows.

Let us assume that clutter is subject to a deterministic and computable flow ϕ c . The randomness of the motion of red blood cells in vessels is modeled by a stochastic differential equation, given by

dy = v b (t, y) dt + σ(y) dB t , (2.9)
where B t is a two dimensional Brownian motion and σ is determined by the effective diffusion coefficient K = 1 2 σ 2 . In blood vessels, this diffusion coefficient is proportional to the product γr 2 where γ is the shear stress in the vessel, and r is the radius of red blood cells. As in the previous section, let a c = a 1,c and a b = a 1,b . Let ϕ b be the flow associated to (2.9). We assume that ϕ b represents the dynamics of blood particles, relative to overall clutter movement, so that

a c (t) = ϕ c (u c , t) , ϕ c (u c , 0) = u c , (2.10)
and

a b (t) = ϕ c (ϕ b (u b , t) , t) , ϕ b (u b , 0) = u b . (2.11)
The dynamics of all the other particles are then taken to be independent realizations of the same dynamics. The velocity field v b and the clutter dynamics ϕ c are computed beforehand and correspond to the general blood flow velocity and to an elastic displacement, respectively. In our experiments, we let ϕ c be an affine displacement of the medium, changing over time: a global affine transformation, with slowly varying translation and shearing applied to the medium at each frame, namely

ϕ c (u, t) = 1 w 1 (t) 0 1 u + w 2 (t) w 3 (t) ,
where w i are smooth and slowly varying (compared to ϕ b ) functions such that w i (0) = 0. As for the blood velocity flow v b , it is parallel to the blood vessels, with its intensity decreasing away from the center of the blood vessel [START_REF] Thomas | Diagnostic Ultrasound Imaging: Inside Out[END_REF]Section 11.3]. More precisely, v b is a Poiseuille laminar flow, namely the mean blood flow velocity is half of the maximum velocity, which is the fluid velocity in the center of the vessel.

The relative blood displacements b k,j = ϕ k,b (u b,k , t j ) are computed according to the following discretization of the stochastic differential equation (2.9):

b k,j+1 = b k,j + δtv b (t j , b k,j ) + √ δtσ (b k,j ) X k,j + o (δt) ,
where (X k,j ) are centered independent Gaussian random variables and δt = t j+1t j is taken to be constant. The blood particle positions a k,b (t j ) are then computed simply by applying the precomputed flow ϕ c . In order to validate the SVD approach, we explore the effects of the blood velocity and of the direction of the blood vessels on the behavior of the singular values and on the quality of the reconstruction. In each case, the clutter displacement is the same composition of time-varying shearing and translation, and the mean clutter velocity is 1 cm•s -1 . We choose C c = 5 and C b = 1, for the same density of scatterers from clutter and blood: per unit of area, the clutter intensity is therefore five times higher than the blood intensity. A single frame of ultrafast ultrasound imaging is presented in Figure 2.3: it is clear that without further processing, it is impossible to locate the blood vessels. In Figure 2.4, the results for various velocities and orientations are presented. The reconstruction intensities are expressed in decibels, relatively to the smallest value in the image. The SVD method allows for reconstruction of blood vessels, even if the maximum blood velocity is close to, or oven lower than, the mean velocity of clutter. We always use the threshold K = 20. As we can see, due to the better resolution in the z direction discussed in Section 1.3, vessels oriented parallel to the receptor array have a reconstruction with a better resolution. But due to the oscillating behavior of the PSF in the z direction, and the low-pass filter behavior of the PSF in the x direction, the sensitivity is better for vessels oriented perpendicularly to the receptor array, and the SVD method is able to reconstruct smaller vessels with lower velocities. This follows from the discussion in Subsection 2.2.2. In order to visualize this phenomenon even better, Figure 2.5 presents the time behavior of a single pixel from the data of Figure 2.4c. We can clearly see the Doppler effect in the case when the flow is perpendicular to the receptor array, and the low frequency behavior of the signal in the case when it is parallel to the receptor array.

In Figure 2.6, results of an investigation on the effect of the threshold K on the reconstruction are presented. Except for K, the parameters of Figure 2.4b are used. If the threshold is too low, the reconstruction is not satisfactory and artifacts appear everywhere in the reconstructed image. If the threshold is too high, the reconstruction still works but the contrast becomes lower. With our parameters, K = 20 seems to produce the best results.

In order to further validate the method, we consider the impact of mea- surement noise on the recovery. To this end, we add independent white Gaussian noise to the data, and consider the quality of the reconstruction as a function of the noise intensity. Let us define the contrast of the reconstruction as the ratio between the mean intensity of the reconstructed image inside and outside the blood domain. The parameters of Figure 2.4b are used. Blood intensity is five times lower than clutter intensity, and therefore a noise intensity of 10% corresponds to half the intensity of blood. In Figure 2.7, sample reconstructions at different noise levels are provided. We can conclude that contrast is robust to moderate levels of noise, since blood vessels can still be identified up to 7.5% of noise if they are oriented along the z axis, and up to 2.5% of noise if they are oriented along the x axis. Figure 2.7 also clearly quantifies the better contrast for vessels oriented along the z axis.

Conclusion

In this chapter, by using a random model for the movement of the blood cells, we have shown that a SVD approach can separate the blood signal from the clutter signal. Our model and results open a door for a mathematical and numerical framework for realizing super-resolution in dynamic optical coherence tomography [START_REF] Joo | Diffusive and directional intracellular dynamics measured by field-based dynamic light scattering[END_REF], in ultrafast ultrasound imaging by tracking micro-bubbles [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF], as well as in acousto-optic imaging based on the use of ultrasound plane waves instead of focused ones, which allows to increase the imaging rate drastically [START_REF] Laudereau | Ultrafast acousto-optic imaging with ultrasonic plane waves[END_REF].

Chapter 3

Dynamic super-resolution

Introduction

The super-resolution problem is the following: given possibly noisy low frequency measurements of a medium -Fourier coefficients below a certain value, convolution by a low pass filter -is it possible to reconstruct the original medium with a precision which exceeds the diffraction limit? This problem is impossible in the general case, but there can be situations where it is indeed possible, with conditions on sparsity for example. One possible mathematical formulation for the super-resolution problem can be the so called sparse spike reconstruction problem: let

µ 0 = N i=1 α i δ x i be a complex measure with finite support defined on Ω ⊂ R d . Let F : M(Ω) → R n ,
where M(Ω) is the set of complex measures on Ω, be the measurement operator, so that the observed vector y is

y = Fµ 0 .
The super-resolution is then that of recovering the measure µ 0 given the measurements y. Since M(Ω) is infinite dimensional, F is not injective, and therefore one has to use regularization to invert it. A common method for solving such sparse problems is to use an infinite dimensional variant of the Lasso program: μ ∈ argmin ν ν TV such that Fν = y.

(3.1)

Mathematical theory on the sparse spike recovery given low frequency measurements has been flourishing in the past years. It includes stable reconstruction of spikes with separation in one and multiple dimensions [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF], robust recovery of positive spikes in the case of a Gaussian point spread function, with no condition of separation [START_REF] Bendory | Robust recovery of positive stream of pulses[END_REF], exact reconstruction for positive spikes in a general settings in one dimension [START_REF] De | Exact reconstruction using beurling minimal extrapolation[END_REF], with estimations on the stability [START_REF] Denoyelle | Support recovery for sparse super-resolution of positive measures[END_REF]. The resolution of ultrafast ultrasound is determined by the wavelength of the incident wave, and by other factors such as the length of the receptor array and the range of angles used in angle compounding [START_REF] Giovanni S Alberti | Mathematical analysis of ultrafast ultrasound imaging[END_REF]. Due to diffraction theory, the minimum resolution one can obtain is of the order of half a wavelength, which is of the order of 300 nm. This implies that in imaging of small blood vessels, blood vessels separated by less than 300 nm cannot be distinguished.

Localization microscopy has been shown to surpass the diffraction limit in optical imaging [START_REF] Michael A Thompson | Extending microscopic resolution with single-molecule imaging and active control[END_REF], [START_REF] Samuel T Hess | Ultrahigh resolution imaging by fluorescence photoactivation localization microscopy[END_REF], and a similar technique has been proposed in ultrasound imaging [START_REF] Errico | Ultrafast ultrasound localization microscopy for deep super-resolution vascular imaging[END_REF], [START_REF] Desailly | Sono-activated ultrasound localization microscopy[END_REF]. In ultrasound imaging, this method consists in using a contrast agent -micro-bubbles for example -which are randomly activated in blood vessels and produce very localized spikes in the observations. If separated by at least several wavelengths, using sparse recovery methods, it is possible to achieve sub-wavelength recovery of the position of particles.

This motivates the introduction of our dynamic model for super-resolution. Instead of considering a single measure µ 0 , we consider a time-varying measure µ t , where t ∈ [-δ, δ]. To simplify the problem, we assume that each point is moving with a constant velocity:

µ t = N i α i δ x i +v i t , t ∈ [-δ, δ],
where v i ∈ R d . The measurements vector is then composed of the measurements at discrete times

t k = kτ, k ∈ [-K, K],
where Kτ = δ:

y k = Fµ t k , k ∈ [-K, K].
In this work, we show that under certain conditions, we are able to recover simultaneously the positions x i and the velocities v i with infinite precision, using a sparse spike recovery based method. Figure 3.1 illustrates the idea of the method in the case of 2d particles and a convolution operator.

The chapter is organized as follows: first, we present theoretical results on the dynamic super-resolution problem. Second, we present numerical

CHAPTER 3. DYNAMIC SUPER-RESOLUTION and Gν = ( ν, ϕ l,k ) l,k , ν ∈ M(Ω x × Ω v ).
In this case, the measurements y = (y k ) k where y k = Fµ t k are given by

y = Gµ.
This allows us to propose the following program to recover positions and velocities:

μ ∈ argmin ν ν TV such that Gν = y. (3.2)
We will call (3.2) the dynamic recovery, whereas (3.1) will be called the static recovery.

The perfect low-pass case

From now on to the end of this section, we consider the specific case of Fourier coefficients:

Ω x = [0, 1], Ω v = [-v max , v max ], ϕ l (x) = e -2π x,l , l ∈ {-f c , . . . , f c } d .
Analog to [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF], a sufficient condition for the solution to (3.2) to be µ is that for any η ∈ C N such that |η| = 1, ∀j ∈ 1, . . . , N , there exists a dual polynomial, called a dual certificate,

q(x, v) = K k=-K ||l||∞≤fc c k,l e i2πl•(x+kτ v) . (3.3) obeying: q(x i , v i ) = η j , j ∈ {1, . . . , N } q(x, v) < 1, everywhere else. (3.4)
Compared to the dual polynomials in the two dimensional case presented in [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF], the only difference with our case, when d = 1, is in the allowed frequencies: whereas in the two dimensional case, all 2 dimensional low frequencies are allowed, here the only allowed frequencies are

(l t , kτ l t ), k = -K, . . . , K, ||l|| ∞ ≤ f c .
The following proposition presents the uniqueness of the reconstruction in the case where we have uniqueness in more than 3 frames:

Proposition 1. Let ν i = α i |α i | .
Assume that there exist m distinct timesteps, m ≥ 3, k 1 , k 2 , . . . , k m such that for all p ∈ {1, . . . , m}0, there exist a trigonometric polynomial q j defined on R d with frequencies in {-f c , . . . , f c } d such that

q p (x i + t kp τ v i ) = ν j , i ∈ {1, . . . , N }, p ∈ {1, . . . , m} |q j (x)| < 1, everywhere else. (3.5)
Moreover, assume that for any m distinct indices i 1 , i 2 , . . . , i m ∈ {1, . . . , N } such that η ip = η iq , for all p, q,

m p=i (x, v) : (x -x ip ) + k p τ (v -v ip ) = 0 = ∅. (3.6)
Then the polynomial q(x, v) = 1 m m p=1 q p (x + t kp τ v) satisfies (3.3) and (3.4). Proof. It is immediate that q defined as such verifies:

q(x i , v i ) = ν j , i ∈ {1, . . . , N } |q(x, v)| ≤ 1, everywhere. (3.7)
Moreover, let (x, v) be such that |q(x, v)| = 1. By the properties of q j , this means that each of the terms in the sum must be equal and have modulus 1/m. Therefore, there must exist i 1 , i 2 , . . . , i m such that η ip = η iq for all p, q and: x + k p τ v = x ip + k p τ v ip , for all p ∈ {1, . . . , m}.

By (3.6), i 1 , i 2 , . . . , i m cannot be distinct. Therefore, at least two are equal to some j ∈ {1, . . . , N } and we can conclude that x = x j and v = v j .

Geometrically, condition (3.6) ensures that there cannot be a virtual point hiding between a true point at each of the frames where there is separation. Numerical experiments show indeed that such situations produce instability in the reconstruction. Figure 3.2 illustrates the simplest of these situations: three points are static (their velocity is 0) and are equally spaced by ∆. In this case, two ghost points arise. One point at the position of the center point at t = 0 and with a velocity v = ∆ /τ, this is the point illustrated in the figure. The other one is the symmetric with the same location but negative velocity.

There are also simple conditions that imply that the condition is always verified: if m > n or if all α i s have a different phase. This simple result ensures that our proposed method works at least as well as what we call the static method, which consists in identifying the position of particles in each frame and then estimating their velocity. In the following proposition, we list a number of concrete cases where we can apply this result.

|x i + t kp v i -x j -t kp v j | > C f c , where C = 2 if d = 1 and C = 2.38 if d = 2.
ii) d = 1, α i > 0, ∀i and the x i + t kp v i are all distinct.

Then µ is the only solution to (3.2).

Proof. Condition i) comes from [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF], condition ii) comes from [START_REF] De | Exact reconstruction using beurling minimal extrapolation[END_REF].

Numerical simulations

Methods

Solving minimization problem (3.2) in all its generality is not an easy task. Indeed, it is not linear and infinite dimensional. It is possible to use an analog discrete problem, where the locations and velocities are fixed on a grid, whose size determines the resolution we want to obtain, and replacing the TV norm by a l 1 norm. However, this methods becomes intractable for a fine resolution. In [START_REF] Emmanuel | Towards a mathematical theory of super-resolution[END_REF], a continuous solution is obtained in the one dimensional case, however, its not trivial to adapt their algorithm to our method.

NUMERICAL SIMULATIONS 51

In this work, we propose to use a continuous reconstruction method. In [START_REF] Boyd | The alternating descent conditional gradient method for sparse inverse problems[END_REF], an algorithm is presented to solve the following problem:

min μ F μ -Y 2 2 subject to μ T V ≤ M. (3.8)
The proposed algorithm is limited to positive weights, but this is a realistic expectation in the case of physical signals. Let us assume that we can guess the value of µ T V beforehand. Then we can set M = µ T V . In this setting, we have the following proposition:

Proposition 3. Assume that µ is the unique solution of (3.2). Then µ is the unique solution of (3.8).

Proof. Since (3.2) admits a unique minimizer, every μ = µ such that F μ = Y verifies μ T V > K. Therefore, μ is the unique minimizer of (3.8).

In order to evaluate the numerical experiments, we introduce the superresolution factor in space as follow:

SRF x = x max f c ∆ x , (3.9)
where x max is such that Ω = [0, x max ]. In the case of velocity reconstruction, we introduce a super-resolution factor in velocity defined as:

SRF v = x max f c ∆ v δ .
(3.10)

1D Fourier examples.

We consider the perfect low-pass filter described in section 3.2. We consider a generic example where a number N of points are scattered in the medium with uniform and independent positions and velocities. The parameters for these simulations are the following:

f c = 20 x max = 1cm K = 2 τ = 1 /60s. (3.11)
The number of points in taken between 4 and 10 and their velocities are taken between + andcm.s -1 . We then realize this setup 1000 times and try to reconstruct using the minimization scheme described above.

APPLICATIONS TO ULTRAFAST ULTRASOUND

55

of static frames where the reconstruction is correct, whereas numerical simulations show that dynamic reconstruction works in cases where the static reconstruction works in no single frame.

Applications to ultrafast ultrasound

This section is devoted to apply the described method to the case of blood flow imaging in ultrafast ultrasound (UFUS). The goal is to locate blood vessels and compute the velocity of blood inside them.

In the case of 2D UFUS with angle compounding, the point spread function is given up to a constant by

g(x, z) = χ 2ν 0 c -1 0 z exp 4iπν 0 c -1 0 z sinc 2πν 0 c -1 0 F x sinc 2πν 0 c -1 0 Θx
, where F, Θ are constants, ν 0 is the base pulsation, and χ is a smooth lowfrequency function [START_REF] Giovanni S Alberti | Mathematical analysis of ultrafast ultrasound imaging[END_REF]. To remove high frequencies in measurements, we demodulate the signal by multiplying our measurements by exp -4iπν 0 c -1 0 z . Therefore, if we sample our image at points (x l , z l ) l , the associated function ϕ l are given by

ϕ l (x, z) =χ 2ν 0 c -1 0 z -z l exp 4iπν 0 c -1 0 z × sinc 2πν 0 c -1 0 F (x -x l ) sinc 2πν 0 c -1 0 Θ(x -x l ) , (3.13)
Thanks to the Shannon sampling theorem, since the demodulated signal (3.13) has almost all its energy in low frequencies, the solution to the lasso minimization will be near the correct measure µ, as long as there is enough stability. This is the case when there is enough separation in the static images.

Figure 3.6 illustrates reconstruction in the case of an UFUS sequence on a small 2D patch.

Conclusion

In this chapter, we have shown that perfect reconstruction of positions and velocity is possible with low frequency measurements. We have proven that in the case where we have separation of the points in the static images, the minimization procedure exactly recovers positions and velocity. We have also shown that this procedure can be applied to ultrafast ultrasound imaging.

What remains to be done is to determine a tighter criterion to have perfect recovery, since numerical simulations show that reconstruction works even in the case where we have no separation in still images.

Chapter 4

Spectrography of cell cultures

Introduction

Cell culture production processes, such as those from stem cell therapy, must be monitored and controlled to meet strict functional requirements. For example, a cell culture of cartilage, designed to replace that in the knee, must be organized in a specific way.

Hyaline cartilage is located on the joint surface and play an important role in body movement. In normal articular cartilage, there is a depthdependent stratified structure known as zonal organization. As a simplified model, cartilage comprises three different layers [START_REF] Joseph M Mansour | Biomechanics of cartilage. Kinesiology: the mechanics and pathomechanics of human movement[END_REF]: a superficial zone in outer 10%, a middle zone that is 50% of the height, and a deep zone consisting in the inner 40%. At the microscopic level, cartilage tissue is composed of cells, collagen fibers, and glycosaminoglycans (GAGs). The concentration and organization of each micro-structure differs among the three layers. In the superficial zone, cells are anisotropic and horizontally aligned, collagen orientation is also horizontal and GAGs have a lower concentration than in the other layers. In the middle zone, there are fewer cells and they are isotropic, collagen is randomly oriented and there is a medium concentration of GAGs. In the deep zone, cells are isotropic, cell density is higher than in the middle zone, collagen is vertically aligned and there is a high GAG density. As these parameters all contribute to the function of collagen in the knee, and must be replicated in the cell culture.

It is important that the method for monitoring cell cultures is nondestructive. Destructive methods require hundreds of samples to be cultured for a single functional tissue, and for the samples to be monitored multiple times during maturation. Here, we propose a microscopic electrical impedance tomography (micro-EIT) method for monitoring cell cultures 58 CHAPTER 4. SPECTROGRAPHY OF CELL CULTURES that exploits the distinctive dielectric properties of cells and other microstructures. In this method, electrodes inject a current into the medium at different frequencies and the corresponding dielectric potentials are recorded, thus enabling reconstruction of the microscopic parameters of the medium. The parameters of interest are cell density, collagen orientation, and GAG density, as well as the orientation and shape of cells.

EIT uses a low-frequency current (below 500 kHz) to visualize the internal impedance distribution of a conducting domain such as a tissue sample or the human body. Recent studies measured electrical conductivity values and anisotropy ratios of engineered cartilage to distinguish extracellular matrix samples containing differing amounts of collagen and GAGs. During chondrogenesis over a six-week period, these measurements could distinguish the stages of the process and provide information regarding the internal depthdependent structure.

In this work, we provide a mathematical framework for determining the microscopic properties of the cell culture from spectral measurements of the effective conductivity. For simplicity, we consider a micro-structure comprising two components in a background medium. One of the components has a frequency dependent on the material parameters arising from the cell membrane structure, while the other has constant conductivity and permittivity over the frequency range. First, we derive in Theorem 4 the overall electrical properties of the culture, which depend on the volume fraction of each component and associated membrane polarization tensors defined by (4.10) and (4.11). Then, we show that the spectral measurements of the overall electrical properties of the culture can be used to determine the volume fraction of each component and the anisotropy ratio of the first component. For doing so, we study the dependence of the membrane polarization tensors on the operating frequency and use the spectral theorem to recover in Proposition 7 from the measurement of the effective conductivity on a range of frequencies the coefficients of its expansion with respect to the frequency. Proposition 7 also provides the anisotropy ratio of the cell culture. This chapter is organized as follows: In Section 4.2, we present a simplified model of the tissue culture. In Section 4.3, we derive an equivalent effective conductivity for the solution at the macroscopic scale. In Section 4.4, we present a method based on spectral measurements, in which microscopic properties are measured from the effective conductivity. This process is known as inverse homogenization or dehomogenization. Finally, we provide some numerical examples to illustrate our main findings.

THE DIRECT PROBLEM

The direct problem

In this section, we propose a simple electrical model for the tissue and derive an effective conductivity using periodic homogenization.

Problem setting

We consider the domain of interest -the cell culture -to be described by a domain Ω ⊂ R 3 . We assume that Ω = D × (0, 1) where D ⊂ R 2 denotes a floor of the culture medium. Following [START_REF] Milton | The theory of composites[END_REF], we describe the conductivity of the medium by a scalar field

σ ω,ε (x) = σ ω x, x ε ,
where ω denotes the angular frequency of the injected current, and ε > 0 is a small parameter representing the microscopic scale of the medium; σ is 1-periodic in every direction in the second variable. Let us consider the following unit domain:

Y = - 1 2 ; 1 2 d .
For a fixed x, σ x, x ε describes the conductivity in a single cartilage tissue with cell size ε at a location x ∈ Ω. To have a complete model of the tissue, σ must describe the conductivity of both cells and of the other inclusions, i.e., collagen and GAGs. The biological fluid conductivity is noted k 0 and is assumed to be frequency independent. The cells are made of biological fluid enclosed in a very thin and very resistive membrane [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF] of thickness εδ for some small parameter δ > 0. The conductivity of the membrane is frequency dependent and is noted k m (ω). The cell shape varies slowly with the parameter x ∈ Ω compared to the micro-scale ε. The other inclusions are described by some frequency independent conductivity function k

k m k 0 k 0 + k i
i x, x ε . Let ψ : Ω × R d :→ R be a C 1 Ω × R d function, 1
-periodic in every direction with respect to the second variable. We assume that the function ψ is the level set function for the membrane boundary given by

Ω + ε = x : ψ x, x ε > δ (resp. Ω - ε = x : ψ x, x ε < -δ ).
We also assume that the support of k i (x, y) is strictly included in {(x, y) : ψ (x, y) > δ}. We can now describe the conductivity σ ω , which is schematically represented at a fixed x in Figure 4.2:

σ ω (x, y) =      k 0 + k i (x, y) if ψ (x, y) > δ, k 0 if ψ (x, y) < -δ, k m (ω) else. (4.1)
Now that we have an expression for the conductivity in the medium, as commonly accepted in EIT, we use the quasi-static approximation for the electrical potential. For an input current g(x) sin (ωt) on the boundary ∂Y, with ∂Ω g = 0, the real part of the corresponding time-harmonic potential, denoted by u ω,ε , satisfies the following problem approximately:

∇ • σ ω,ε ∇u ω,ε = 0 in Ω, σ ω,ε ∇u ω,ε • ν = g on ∂Ω. (4.2)
where ν is the outer normal vector on ∂Ω. Here, we impose the normalization Ωε u ω,ε = 0.

Remark 3. Let us briefly explain how the expression of σ ω in (4.1) is derived. We should note that the frequency dependent behaviors of σ ω,ε in (4.2) are attributed to thin cell membranes. Imagine that we inject an oscillating current at the angular frequency ω into the cube Y. Then, the resulting time-harmonic potential w = u + iv in Y is governed by

∇ • ((σ ′ (y) + iωσ ′′ (y))∇w(y)) = 0 for y ∈ Y,
where σ ′ denotes the conductivity distribution and σ ′′ is the permittivity distribution in Y. In [START_REF] Kim | Asymptotic analysis of the membrane structure to sensitivity of frequency-difference electrical impedance tomography[END_REF], it was shown that, under some conditions on the membrane, the real part u approximately satisfies

∇ • |σ ′ + iωσ ′′ | 2 σ ′ ∇u = 0 in Y. (4.3)
Since σ ′′ ≪ σ ′ outside the membrane, we have

|σ ′ + iωσ ′′ | 2 σ ′ ≈ σ ′ outside the membrane.
Hence, it is reasonable to assume that the conductivity outside the membrane, as a coefficient of the elliptic PDE (4.3), does not change with frequency. On the other hand, since σ ′ on the membrane is very small, the effect of σ ′′ is not negligible. Hence, the conductivity, k m , on the membrane changes with frequency as follows:

|σ ′ + iωσ ′′ | 2 σ ′ = σ ′ + ω 2 σ ′′ σ ′
on the membrane.

Homogenization of the tissue

We are now interested in getting rid of the micro-scale oscillations of σ ω,ε , since boundary measurements will only allow us to image macro-scale variations of the conductivity. To this end, we proceed to the homogenization of equation (4.2). Assume that k 0 + k i is bounded from below and from above:

0 <σ ≤ k 0 + k i ≤ σ.
From [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF], we have two-scale convergence [START_REF] Allaire | Homogenization and two-scale convergence[END_REF][START_REF] Milton | The theory of composites[END_REF][START_REF] Nguestseng | A general convergence result for a functional related to the theory of homogenization[END_REF] of u ω,ε to u ω , which is a solution to

     ∇ • σ * ω ∇u ω = 0 in Ω, σ * ω ∇u ω • ν = g on ∂Ω, Ω u ω = 0, (4.4)

CHAPTER 4. SPECTROGRAPHY OF CELL CULTURES

for an input current g(x) sin (ωt) on the boundary ∂Ω. Here, σ * ω is called the effective conductivity which can be represented by [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF] σ * ω (x)e p • e q = Y σ ω (x, y) ∇ (y p + v p (y)) • e q dy, ∀p, q ∈ {1, ..., d} = k 0 δ p,q + ∂Y ∂v p ∂ν y q ds(y) , As δ → 0, v p can be approximated [START_REF] Khelifi | Asymptotic expansions for the voltage potentials with two-dimensional and three-dimensional thin interfaces[END_REF] by the solution of the following equation, where β (ω) = δ km(ω) :

               ∇ • (σ ω (x, y) ∇(v p (y) + y p )) = 0 for y ∈ Y\∂C, k 0 ∂ ∂ν (v + p (y) + y p ) = k 0 ∂ ∂ν (v - p (y) + y p ) for y ∈ ∂C, v + p (y) -v - p (y) = β (ω) k 0 ∂ ∂ν (v + p (y) + y p ) for y ∈ ∂C, v p 1-periodic, Y v p (y) + ydy = 0. (4.7)
Here, ∂C denotes the membrane of the cell C and βk 0 is the effective thickness of the membrane.

Imaging the micro-structure from effective conductivity measurements

In this section, we do not care about the space dependence of σ * ω , and will therefore drop it. We will thus assume that σ * ω is constant equal to some matrix in M d (C) := {m ∈ C d×d : m i,j = m j,i for i, j = 1, 2, • • • , d}. We will show what kind of information on the micro-structure we can recover from the knowledge of σ * ω in a range of frequencies ω ∈ (ω 1 , ω 2 ). First, in section 4.3.1, we will obtain a simple representation of the effective conductivity in the dilute case, where the volume fraction of both cells and other inclusions is small compared to the volume of biological fluid. Then, in the following sections we will use this representation and will show how to recover information about the micro-structure using the spectral measure.

IMAGING THE MICRO-STRUCTURE

Effective conductivity in the dilute case

Here, we consider some reference cell C 0 and some reference inclusion B 0 with there C 2 boundaries ∂C 0 and ∂B 0 . We assume that C = x C + ρ C C 0 and β (ω) = ρ C β 0 (ω) for some reference β 0 (ω) and let B = x B + ρ B B 0 , where x C and x B respectively indicate the locations of the cell and inclusion and ρ C and ρ B their characteristic sizes. We assume that the conductivity k i of the inclusion is given by

k i (y) = (k 0 -k 1 ) χ B (y) ,
where χ B denotes the characteristic function of B.

The effective conductivity is therefore expressed as

σ * ω e p • e q = Y σ (y) ∇ (y p + v p (y)) • e q dy, ∀p, q ∈ {1, • • • , d} , where, for p ∈ {1, • • • , d},                              ∇ • (k 0 ∇(v p (y) + y p )) = 0 in Y\ (B ∪ ∂C) , ∇ • (k 1 ∇(v p (y) + y p )) = 0 in B, k 0 ∂ ∂ν (v + p (y) + y p ) = k 0 ∂ ∂ν (v - p (y) + y p ) on ∂C, v + p -v - p = β (ω) k 0 ∂ ∂ν (v + p (y) + y p ) on ∂C, v + p -v - p = 0 on ∂B, k 0 ∂ ∂ν (v + p (y) + y p ) = k 1 ∂ ∂ν (v - p (y) + y p ) on ∂B, v p periodic, Y (v p (y) + y) dy = 0. (4.8)
From now on, I denotes the inclusion map H 1/2 (∂C) → H -1/2 (∂C), where H 1/2 and H -1/2 are the Sobolev spaces of order 1/2 and -1/2 on ∂C. We will now proceed to prove the following result.

Theorem 4. Let f k = ρ d k , k ∈ {B, C} and f = max (f B , f C ).
Then we have the following expansion:

σ * ω = k 0 [I + f B M B 0 + f C M C 0 (ω)] + o (f ) , (4.9)
where

M C 0 (ω)e p • e q = ∂C 0 ν q (y) 1 β 0 (ω) k 0 I + L #,C 0 -1 [ν p ](
with λ = k 1 + k 0 2 (k 1 -k 0 )
.

We begin be reviewing properties of periodic layer potentials. Let us define the periodic Green's function

G # (x) = - n∈Z d \{0} e 2iπn•x 4π 2 |n| 2 .
Thanks to Poisson's summation formula, in the sense of distribution,

G # satisfies ∆G # (x) = n∈Z d δ(x -n) -1. (4.12)
We write G # (x, y) := G # (xy). Let us introduce the periodic single layer potential, for a Lipschitz domain D ⊂ Y:

S #,D : H -1/2 (∂D) → H 1 loc R d \∂D ϕ → x → ∂D G # (x, y) ϕ(y)ds(y),
the periodic double layer potential

D #,D : H 1/2 (∂D) → H 1 loc R d \∂D ϕ → x → ∂D ∂G # ∂ν(y) (x, y) ϕ(y)ds(y),
and the periodic Neumann-Poincaré operator

K #,D : H 1/2 (∂D) → H 1/2 (∂D) ϕ → x → ∂D ∂G # ∂ν(y) (x, y) ϕ(y)ds(y),
and its adjoint given by

K * #,D : H -1/2 (∂D) → H -1/2 (∂D) ϕ → x → ∂D ∂G # ∂ν(x) (x, y) ϕ(y)ds(y).
We review the jump properties of the layer potentials [START_REF] Ammari | Polarization and moment tensors: with applications to inverse problems and effective medium theory[END_REF].

IMAGING THE MICRO-STRUCTURE

65

Lemma 1. We have the following jump relations along the boundary ∂D:

S #,D [ϕ](x)| + = S #,D [ϕ](x)| -, ∂ ∂ν S #,D [ϕ](x) ± = ± 1 2 I + K * #,D [ϕ](x), D #,D [ϕ](x)| ± = ∓ 1 2 I + K #,D [ϕ](x), ∂ ∂ν D #,D [ϕ](x) + = ∂ ∂ν D #,D [ϕ](x) - .
where the subscript ± means f D (x)| ± = lim t→0 + f D (x ± tν(x)) for x ∈ ∂D.

We denote by L #,D the operator ϕ → ∂ ∂ν D #,D [ϕ]. We write ν p = ν • e p on ∂B and ∂C. Using these jump relations, we have the following representation theorem for v p , p ∈ {1, ..., d}.

Theorem 5. We have the following representation for v p :

v p = C p + S #,B [ϕ 1,p ] -D #,C [ϕ 2,p ] , (4.13)
where C p is a constant and (ϕ 1 , ϕ 2 ) satisfies the following system:

   λI -K * #,B [ϕ 1,p ] + ∂ ∂ν D #,C [ϕ 2,p ] = ν p on ∂B, 1
βk 0 I + L #,C [ϕ 2,p ] -∂ ∂ν S #,B [ϕ 1,p ] + = ν p on ∂C. (4.14)
Lemma 2. For any (F, G) ∈ H -1/2 (∂B) × H -1/2 (∂C), the system

   λI -K * #,B [ϕ 1 ] + ∂ ∂ν D #,C [ϕ 2 ] = F on ∂B, 1 βk 0 I + L #,C [ϕ 2 ] -∂ ∂ν S #,B [ϕ 1 ] + = G on ∂C, admits a unique solution (ϕ 1 , ϕ 2 ) ∈ H -1/2 (∂B) × H 1/2 (∂C).
Proof. As shown in Appendix B,

1 β I + L #,C and λI -K * #,B are invertible for λ / ∈ (-1/2, 1/2]. Moreover, since ∂ ∂ν D #,C : H 1/2 (∂C) → H -1/2 (∂B) and ∂ ∂ν S #,B : H -1/2 (∂B) → H -1/2 (∂C)
are compact, the operator

H -1/2 (∂Ω) × H 1/2 (∂Ω) → H -1/2 (∂Ω) × H -1/2 (∂Ω) (ϕ 1 , ϕ 2 ) → λI -K * #,B [ϕ 1 ] - ∂ ∂ν D #,C [ϕ 2 ], 1 βk 0 I + L #,C [ϕ 2 ] - ∂ ∂ν S #,B [ϕ 1 ]
+ is a Fredholm operator. It is therefore sufficient to show that it is injective. Let (ϕ 1 , ϕ 2 ) be such that

λI -K * #,B [ϕ 1 ] + ∂ ∂ν D #,C [ϕ 2 ] = 0 on ∂B, 1 βk 0 I + L #,C [ϕ 2 ] -∂ ∂ν S #,B [ϕ 1 ] = 0 on ∂C. Let v = S #,B [ϕ 1 ] -D #,C [ϕ 2 ].
Then v is 1-periodic in every direction, and v is a solution by construction to the following problem: By the uniqueness of the solution to (4.15) up to a constant, v(x) = c, ∀x ∈ Y. Then, we have ϕ 1 = 0 on ∂C and ϕ 2 = 0 on ∂B because they are equal to the jumps of v (resp. ∂v ∂ν ) across ∂B (resp. ∂C). This concludes the proof.

                   ∇ • (k 0 ∇(v p (y) + y)) = 0 for y ∈ Y\ (B ∪ ∂C) , ∇ • (k 1 ∇(v p (y) + y)) = 0 for y ∈ B, k 0 ∂ ∂ν (v + p (y) + y) = k 0 ∂ ∂ν (v - p (y) + y) for y ∈ ∂C, v + p (y) -v - p (
We can now proceed to prove Theorem 5.

Proof. Let (ϕ 1 , ϕ 2 ) be a solution of (4.14), and let

v p = S #,B [ϕ 1 ] -D #,C [ϕ 2 ] .
Then using the jump relations of the layer potentials, we have that v p is a solution of (4.8), except that we have not necessarily ∂Y v p = 0. We just have to adjust C p accordingly.

We now proceed to compute the representation of the effective conductivity.

Theorem 6. We have the following representation for σ * ω :

σ * ω = k 0 (I + M * ) ,
where M * = (M * pq ) d p,q=1 is defined by

(M * ) pq = ∂B
x p ϕ 1,q ds-∂C ν p ϕ 2,q ds, ∀p, q ∈ {1, ..., d} .

Proof. We recall the expression of σ * ω in (4.5):

σ * ω e p • e q = k 0 δ p,q + ∂Y ∂v p ∂ν (y)y q ds(y) .

Using representation (4.13), we obtain The same reasoning applies to the second part of the equation:

∂Y ∂D #,C [ϕ 2,p ] ∂ν (y)y q ds(y) = ∂C D #,C [ϕ 2,p ]| + (y)ν q (y)ds(y) - ∂C D #,C [ϕ 2,p ]| -(y)ν q (y)ds(y) = ∂C ϕ 2,p (
y)ν q (y)ds(y).

Therefore, σ * ω e p • e q = k 0 δ p,q + ∂Y ∂v p ∂ν (y)y q ds(y)

= k 0 δ p,q + ∂B y q ϕ 1,p (y)ds(y)-∂C ϕ 2,p ν q (y)(y)ds(y) .

We turn to the proof of Theorem 4. We first review asymptotic properties of the periodic Green's function G # . The following result from [9, Chapter 2] holds. Lemma 3. We have the following expansion for G # :

G # (x) = G(x) + R d (x),
where G is the Green function and R d is a smooth function on R d and its Taylor expansion at 0 is given by

R d (x) = R d (0) - 1 2d |x| 2 + O |x| 4 . (4.16)
Using this expansion, we obtain by exactly the same arguments as those in [9, Chapter 8] the following expansion, which is uniform in z ∈ ∂B 0 ,

λI -K * B 0 [ψ B,p ](z) = ν B 0 ,p (z) + o(1) 1 β 0 k 0 I + L C 0 [ψ C,p ](z) = ν C 0 ,p (z) + o(1),
where K * B 0 is the standard Neumann-Poincaré operator and L C 0 denotes the operator ∂ ∂ν D C 0 associated with the standard double layer potential D C 0 :

K * B 0 [φ](x) := ∂B 0 ∂G ∂ν(x) (x, y)φ(y)ds(y), L C 0 [φ](x) := ∂ ∂ν ∂C 0 ∂G ∂ν(y)
(x, y)φ(y)ds(y).

Therefore, we arrive at the result stated in Theorem 4.

Spectral measure of the tissue

Expansion (4.9) yields

σ * ω = k 0 I + ρ d B M B 0 + ρ d C M C 0 (ω) + O ρ d with M C 0 (ω)e p • e q = ∂C 0 ν q (y) 1 β 0 (ω) k 0 I + L C 0 -1 [ν p ](y) ds(y).
In order to use the spectral theorem in a Hilbert space, we have to modify the expression of

M C 0 . Let L -1 C 0 be the inverse of L C 0 : H 1/2 0 (∂C 0 ) → H -1/2 0 (∂C 0 ). Then we write 1 β 0 (ω) k 0 I + L C 0 -1 [ν p ] = 1 β 0 (ω) k 0 L -1 C 0 • I + I H 1/2 -1 L -1 C 0 [ν p ].
The following result holds.

IMAGING THE MICRO-STRUCTURE

69 Lemma 4. L -1 C 0 •I can be extended to a self-adjoint operator L † : L 2 (∂C 0 ) → L 2 (∂C 0 ) , whose image is a subset of H 1/2 (∂C 0 ). Proof. Let J 1 : L 2 (∂C 0 ) ֒→ H -1/2 (∂C 0 ) and J 2 : H 1/2 (∂C 0 ) ֒→ L 2 (∂C 0 ). Let L † = J 2 • L -1 C 0 • J 1 . Then obviously L † extends L -1 C 0 • I and its image is a subset of H 1/2 (∂C 0 ). Let us show that it is self-adjoint. Let (ϕ, ψ) ∈ L 2 (∂C 0 ) × L 2 (∂C 0 ). Let , L 2
and , H 1/2 ,H -1/2 respectively denote the L 2scalar product and the duality pairing between H 1/2 (∂C 0 ) and H -1/2 (∂C 0 ). We have

L † [ϕ], ψ L 2 = L -1 C 0 [ϕ], ψ L 2 = L -1 C 0 [ϕ], ψ H 1/2 ,H -1/2 = L -1 C 0 [ψ], ϕ H 1/2 ,H -1/2 = L -1 C 0 [ψ], ϕ L 2 = L † [ψ], ϕ L 2 , since L C 0 is self-adjoint from H 1/2 (∂C 0 ) onto H -1/2 (∂C 0 ).
From this result, we can now proceed. From the spectral theorem, there exists a spectral measure E such that for any z ∈ C \ Λ L † and for any

(ϕ, ψ) ∈ (L 2 (∂C 0 )) 2 , L † z + I -1 [ϕ], ψ L 2 = Λ(L † ) 1 x z + 1 ϕ(x)ψ(x)dE (x) . (4.17)
where Λ(L † ) denotes the spectrum of L † . Let

F p,q (z) = δ p,q +ρ d B M B 0 e p •e q +ρ d C Λ(L † ) 1 x z + 1 L -1 C 0 [ν p ](x)•ν q (x)dE (x) . (4.18)
where δ p,q = 1 if p = q and δ p,q = 0 if p = q. Therefore, we have

σ * ω e p • e q ≃ k 0 [F p,q (β 0 (ω) k 0 )] . Since lim z→0 F (z) = I + ρ d B M B 0 ,
there is no singularity of F in 0. Since 0 / ∈ Λ L † , (4.17) is valid on a neighborhood of 0. Proposition 7. Let F = (F p,q ) d p,q=1 be defined by (4.18). Then the following expansion of F in a neighborhood of 0 holds:

F p,q (z) = ∞ k=0 a k,p,q z k , (4.19)
where a 0,p,q = I + ρ d B M B 0 e p • e q , and a 1,p,q = ρ d C ν p • ν q .

Proof. Identity (4.19) holds using the analyticity of F in a neighborhood of 0. We also have

a 0,p,q = lim z→0 F p,q (z) = δ p,q + ρ d B M B 0 e p • e q .
In order to obtain the next coefficients, we begin by establishing the following limit:

lim z→0 L † + zI -1 [ν p ] = L C 0 [ν p ], p = 1, 2. Indeed, let ϕ (z) = L † + zI -1 [ν p ]. Then ϕ (z) = 1 z ν p -L † ϕ p . Since the range of L † is a subset of H 1/2 (∂C 0 ), ϕ (z) ∈ H 1/2 (∂C 0 ). Therefore, ϕ (z) = L C 0 [ν p ] -zL C 0 [ϕ](z) → z→0 L C 0 [ν p ].
This yields

lim z→0 1 z (F p,q (z) -F p,q (0)) = ρ d C ν p • ν q .
In the following, we write

F (z) = (F p,q (z)) p,q∈{1,...,d} , z ∈ C\Λ (L C 0 ) , and
A k = (a k,p,q ) p,q∈{1,...,d} , k ∈ N. (4.20)
Since F p,q is analytic on C\Λ (L C 0 ), the values of a k can be recovered from the values of F p,q on a subset of C with a limiting point. Therefore, we can reconstruct the values a k,p,q from the measurements of the effective conductivity σ * ω in a band of frequencies ω ∈ (ω 1 , ω 2 ). Further details on this will be provided in the following section.

Inverse homogenization

Imaging of the anisotropy ratio

The anisotropy ratio (the ratio between the largest and the lowest eigenvalue of the effective conductivity tensor) depends on the frequency [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF]. Furthermore, in the general case, the anisotropy orientation (the direction of the effective conductivity tensor eigenvectors) depends also on the frequency. However, in the special case where we have an axis of symmetry of a single inclusion or a cell, the anisotropy orientation is independent of the frequency.

We denote by

O d (R) := {R ∈ R d×d | R T R = 1, det(R) = 1}
the set of rotational matrices. Here, the superscript T denotes the transpose. For convenience, we write R(x) := Rx for x ∈ Y and R(D) := {Rx : x ∈ D}. We will need the following covariance result :

Lemma 5. Let R ∈ O d (R) and f ∈ L 2 (∂C 0 ). Then L C 0 [f • R] • R = L C 0 [f ].
Proof. We have, for any x ∈ ∂C 0 ,

L C 0 [f • R] (R (x)) = lim h→0 ∇D C 0 [f • R] (R(x) + hν (R(x))) • ν (R(x)) . Moreover, D C 0 [f • R] (R(x)) = ∂C 0 ∇G (R (x) -y) • ν (y) f (R (y)) ds (y) = ∂C 0 ∇G (R (x) -R (y)) • ν (R (y)) f (y)ds (y) .
Since G is isotropically symmetric,∇G (R (xy)) = R (∇G (xy)), therefore for any x, y ∈ ∂C 0 ,

∇G (R (x) -R (y))•ν (R (y)) = R (∇G (x -y))•R (ν (y)) = ∇G (x -y)•ν (y) so that D C 0 [f • R] (R(x)) = D C 0 [f ] (x) , ∀x ∈ ∂C 0 .
This in turn implies that

L C 0 [f • R] (R (x)) = lim h→0 ∇D C 0 [f • R] (R(x) + hν (R(x))) • ν (R(x)) = lim h→0 ∇D C 0 [f ] (x + ν (x)) • ν (x) = L C 0 [f ] (x).
The following corollary holds immediately.

Corollary 8. Let R ∈ O d (R). Then, M R(C 0 ) = RM C 0 R T .
Let us begin with the two-dimensional case.

Proposition 9. Let d = 2, and (ε 1 , ε 2 ) be an orthonormal basis of R 2 . Let ξ be the orthogonal symmetry of axis ε

1 . If ξ (C 0 ) = C 0 , then F (z)ε 1 • ε 2 = 0, ∀z ∈ C \ Λ(L † ).
Proof.

F (z)ε 1 • ε 2 = ρ d C ∂C 0 L † z + I -1 [ν • ε 1 ] (x) ν(x) • ε 2 ds(x) = ρ d C ∂C 0 L † z + I -1 [ν • ε 1 ] (ξ(x)) ν(ξ(x)) • ε 2 ds(x) = -ρ d C ∂C 0 L † z + I -1 [ν • ε 1 ] (x) ν(x) • ε 2 ds(x) because ν(ξ(x)) • ε 1 = ν(x) • ε 1 and ν(ξ(x)) • ε 2 = -ν(x) • ε 2 . Therefore, F (z)ε 1 • ε 2 = 0, ∀z ∈ C \ Λ(L † ).
We have a similar result in three dimensions. The following proposition holds.

Proposition 10. Let d = 3, and (ε 1 , ε 2 , ε 3 ) be an orthonormal basis of R 3 . Let ξ 1 (resp. ξ 2 ) be the orthogonal symmetry of axis ε

1 (resp. ε 2 ). If ξ 1 (C 0 ) = ξ 2 (C 0 ) = C 0 , then F (z)ε j • ε k = 0, ∀z ∈ C, ∀k = j ∈ {1, 2, 3} .
Proof. The proof is exactly the same as in the d = 2 case and is therefore omitted.

Remark 4. It is also true that the symmetry axes of B 0 correspond to the eigenvectors of the polarization tensor M B 0 . Therefore, the anisotropy direction of the frequency-independent background can also be recovered as the principal directions of M B 0 .

Remark 5. Even if each of inclusion and cell has an axis of symmetry, the direction of eigenvectors of the effective conductivity tensor can be frequency dependent. The following numerical test is conducted to show an example of frequency dependency. There are an ellipsoidal inclusion with major axis e 1 and minor axis e 2 and an ellipsoidal cell with major axis e 2 and minor axis e 1 in the unit square as shown in Figure 4.4 (a). For the square domain Y = (-1 2 , 1 2 ) 2 , each axis length of cell and inclusion is 1/8, and 1/24. The center of ellipsoidal cell and inclusion are (1/3, 1/6) and (0, -1/3) respectively. The ratio between membrane thickness and size of a cell is 5 × 10 -3 . The conductivity value of medium, membrane, inclusion are 0.5 S/m, 10 -5 S/m, and 10 -12 S/m respectively. We use (4.5) to compute the effective conductivity tensor. For the numerical computation, we take advantage of using u j satisfying ∇ • (σ∇u j ) = 0 in Ω with boundary condition u j (y)| ∂Ω = y j | ∂Ω for y = (y 1 , y 2 ). Then, v j can be replaced with v j = u jy j . Hence, the eigenvectors of the effective conductivity can be computed and the main direction of anisotropy changes in terms of the frequency as shown in Figure 4.4 (b). At each frequency, in order to compute the true effective conductivity given by (4.5), we perform a finite element computation using FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF]. Comparison between the true effective conductivity and the expansion from Theorem 4.9 can be seen in Figures 4. [START_REF] Habib Ammari | Spectroscopic imaging of a dilute cell suspension[END_REF] To recover the moments from the effective conductivity, we approximate as a rational function, F p,q (z) ≃ p 0 + p 1 z + ... + p N z N q 0 + q 1 z + ... + q N z N . for some N ∈ N. Such an approximation of F is called a Padé approximation of F . Then we approximate the moments by the following values: ã0,p,q = p 0 q 0 , ã1,p,q = p 1 q 0 -q 1 p 0 q 2 0 .

Implementation of the inverse homogenization

Numerically, this is done as a simple least square inversion: the coefficients of the polynomials P (z) = p 0 +p 1 z+...+p N z N and Q(z) = q 0 +q 1 z+...+q N z N are computed to minimize the quantity

K k=1 F p,q (z k ) - P (z k ) Q(z k ) 2 ,
where z 1 , ..., z K are the frequency values where F is measured.

We now consider a toy example where C is an ellipse in R 2 . In this case, if λ 1 and λ 2 are the eigenvalues of A 1 defined by (4.20) for k = 1, the ratio r := λ 2 /λ 1 is independent of the volume fraction and is given by r = Since the right-hand side of (4.21) can be regarded as a function of a/b, the anisotropy ratio a/b can be easily obtained by solving (4.21) with the known value r. In Figure 4.8 (resp. in Figure 4.9), we illustrate the reconstruction of the ratio r using the Padé approximation of F as a function of the anisotropy ratio a/b compared to its theoretical value given by the preceding formula in the case where there is no inclusion B (resp. with an inclusion B with ρ B = 0.1). As we can see, the reconstruction is almost perfect in the case where there is no inclusion, and there is a slight bias induced by the inclusion B.

After recovering the anisotropy ratio a/b, we can recover the volume fraction ρ C from the product of λ 1 , λ 2 of the eigenvalues of A 1 . Indeed, we have Value of anisotropy 4.1 presents numerical reconstruction of the volume fraction ρ C using the preceding formula, with an anisotropy ratio equal to 2.

λ 1 λ 2 =
To reconstruct the angle of the inclusions, we simply use the orientation of the eigenvalues of the moments of A 0 for B and A 1 for C. This is illustrated by results in Figure 4.10 when both B and C are ellipses of anisotropy ratio 2 and with ρ B = ρ C = 0.1.

Therefore, the absolute error E(x) due to the truncation of the Taylor series of w ± at first order is given by

E(x) = c 0 (4π) -1 [E + (x) -E -(x)] ,
where

E ± (x) = 1 x f ′ ( ∓c -1 0 F √ 1 + F 2 x + c x x 2 2 ) -f ′ ( ∓c -1 0 F √ 1 + F 2 x) .
We now consider two cases, depending on x. First, consider the case when x > 5 • 10 -3 m. From the above calculations we immediately have Next, consider the case when x ≤ 5 • 10 -3 m. By using again the mean value theorem we obtain

|E(x)| ≤ c 0 (4π) -1 4 x f ′ ∞ ≤ 2
E ± (x) = c x x 2 f ′′ (θ x ), θ x = ∓c -1 0 F √ 1 + F 2 x + δ x c x x 2 2
for some δ x ∈ [0, 1]. Since |f ′′ (t)| is even and decreasing for t > 0, we have that

|E ± (x)| ≤ c -1 0 x 2z |f ′′ ( c -1 0 F √ 1 + F 2 x -c -1 0 x 2 2z )|,
since the inequality x ≤ 5 • 10 -3 m guarantees that

c -1 0 F √ 1+F 2 x -c -1 0 x 2
2z > 0. Therefore we have

|E(x)| ≤ (4π) -1 xz -1 |f ′′ ( c -1 0 F √ 1 + F 2 x -c -1 0 x 2 2z )|.
Let us look at the right hand side of this inequality. As x → 0 the error tends to 0: this is expected, because of the Taylor expansion around 0. On the other hand, for big x, the value of |f ′′ (

c -1 0 F √ 1+F 2 x -c -1 0 x 2
2z )| is very small, since |f ′′ (t)| decays very rapidly for large t. Therefore, the maximum of the right hand side is attained in a point x * ∈ (0, 0.005). The value in this point may be explicitly calculated, and we have We now wish to estimate the relative error E ∞ / g ∞ . In order to do this, let us compute g(0). Since the Taylor expansion becomes exact as x → 0, we may very well compute g(0) by using the approximated version. Thus, setting G = F/ √ 1 + F 2 we have

g(0) = lim x→0 - c 0 4πx f ′ (c -1 0 Gx) -f ′ (-c -1 0 Gx) = lim x→0 -G(4π) -1 f ′ (c -1 0 Gx) -f ′ (0) c -1 0 Gx + f ′ (-c -1 0 Gx) -f ′ (0) -c -1 0 Gx = -2G(4π) -1 f ′′ (0),
whence |g(0)| ≥ 8.8•10 13 by a direct calculation of |f ′′ (0)|. Finally, combining this inequality with (A.1) allows to bound the relative error by

E ∞ g ∞ ≤ 5%.
We have proven that the relative error of the approximation obtained by truncating the Taylor expansions of w ± at the first order is less than 5%. This has been proven only in the particular case when z = z ′ : the general case may be done by extending the above argument to two dimensions. Therefore, S #,C is negative semi-definite, which concludes the proof.

Theorem 12. For λ / ∈ -1 2 , 1 2 , the operator λI -K * #,C is invertible on H -1/2 (∂C).

Proof. Since λI -K * C is invertible, K * #,C -K * C is a compact operator [START_REF] Ammari | Polarization and moment tensors: with applications to inverse problems and effective medium theory[END_REF], λI -K * #,C is a Fredholm operator and it is enough to show that it is one-toone. The proof goes exactly as in [START_REF] Chang | Spectral properties of the layer potentials on lipschitz domains[END_REF]. Let us assume that λI -K * #,C is not one-to-one. Then there exists some f ∈ H Then A = 0 or B = 0 since f is not identically zero. Then by Green's formula together with the jump formulas, we have .

A = - 1 2 I + K * #,C [f ], S #,C [f ]

  The real part of the input pulse f . The imaging system. The incident wave is supported near the imaging plane {y = 0}, within the focusing region bounded by the two curved surfaces.

Figure 1 . 1 :

 11 Figure 1.1: The pulse f of the incident wave u i and the focusing region.

  The approximation of the PSF given in (1.9).

  The approximation of the PSF given in (1.10).

Figure 1 . 2 :

 12 Figure 1.2: The real part of the point spread function g 0 and its approximations are shown in these figures (with parameters as in (1.1) and (1.2), and F = 0.4). The size of the square shown is 2 mm × 2 mm, and the horizontal and vertical axes are the x and z axes, respectively. The relative error in the L ∞ norm is about 7% for the approximation shown in panel (b) and about 9% for the approximation shown in panel (c).

20 CHAPTER 1

 201

Figure 1 . 3 :

 13 Figure 1.3: The exact PSF with different values of the aperture parameter F (with parameters as in (1.1) and (1.2), and θ = 0). The size of the square shown is 2 mm×2 mm, and the horizontal and vertical axes are the x and z axes, respectively.

Figure 1 . 4 :

 14 Figure 1.4: The exact PSF with different values of the angle θ (with parameters as in (1.1) and (1.2), and F = 0.4). The size of the square shown is 2 mm × 2 mm, and the horizontal and vertical axes are the x and z axes, respectively.

Figure 1 . 5 :

 15 Figure 1.5: The absolute values of the Fourier transforms of the point spread functions and its approximations (with parameters as in (1.1) and (1.2), and F = θ = 0). The frequency axes are normalized by ν 0 c -10 : the PSF is a low pass filter with cut-off frequency F ν 0 c -1 0 with respect to the variable x and a band pass filter around 2ν 0 c -1 0 with respect to z.

  The PSF g ac Θ with Θ = 0.25.

Figure 1 . 6 :

 16 Figure 1.6: A comparison of the PSF related to the single illumination with the PSF associated to multiple angles (with parameters as in (1.1) and (1.2), and F = 0.4). The better focusing in the variable x for g ac Θ is evident, as well as the good approximation given by gac Θ . The size of the square shown is 2 mm × 2 mm, and the horizontal and vertical axes are the x and z axes, respectively.

33 2. 3

 333 The Inverse Problem: Source Separation 2.3.1 Formulation of the dynamic inverse problem

  The real part of C gg . The real part of C gg . Clutter: jCov(S(0; 0); S(z; t)Clutter: jCov(S(0; 0); S(z; t)Absolute values of the correlations in the clutter model (σ = 0, v c = 10 -2 m•s -1 ) and in the blood model (σ 2 = 10 -6 m 2 s -1 , v b = 10 -2 m•s -1 ).

Figure 2 . 1 :Clutter

 21 Figure 2.1: Correlations of the Casorati matrix.

Figure 2 . 2 :

 22 Figure 2.2: The distribution of the singular values of the Casorati matrix S in different cases.

Figure 2 . 3 :

 23 Figure 2.3: Single frame of ultrafast ultrasound (real part).

  Maximum blood velocity: 2 cm•s -1 ; mean clutter velocity: 1 cm•s -1 . Maximum blood velocity: 1 cm•s -1 ; mean clutter velocity: 1 cm•s -1 . Maximum blood velocity: 0.5 cm•s -1 ; mean clutter velocity:1 cm•s -1 .

Figure 2 . 4 :

 24 Figure 2.4: The SVD method for different velocities and orientations. In each case, we have from left to right: the blood velocity and location, the reconstructed blood location, the decay of the singular values. The squares are 5 mm × 5 mm, and the horizontal and vertical axes are the x and z axes, respectively. The parameters used are those given in (1.1) and (1.2), F = 0.4 and Θ = 7. The density of particles for both blood and clutter is 2,000 per mm 2 , and σ = 2.5 • 10 -5 .

Figure 2 . 5 :

 25 Figure 2.5: Time behavior of a single pixel (real part), located in a constant velocity flow.

  Flow parallel to the receptor array.

  Flow perpendicular to the receptor array.

Figure 2 . 6 :

 26 Figure 2.6: Effect of the threshold K on the reconstruction. From left to right: K = 10, 20, 30, 40.

Figure 2 . 7 :

 27 Figure 2.7: Effect of noise on the reconstruction. The parameters are the same used in Figure 2.4b.

50 CHAPTER 3 Figure 3 . 2 :

 50332 Figure 3.2: Geometry of ghost points

Figure 4 . 1 :

 41 Figure 4.1: Organization of the cells in the cartilage tissue.

Figure 4 . 2 :

 42 Figure 4.2: Typical values of σ ω on Y.

(4. 5 )

 5 where e p := (0,• • • , 1, • • • , 0) with 0 components except p-th component 1,and v p is the solution to the following equation on Y for p = 1, ..., d: σ ω (x, y) ∇(v p (y) + y p )) = 0 for y ∈ Y, v p 1-periodic, Y (v p (y) + y p )dy = 0.

  63

  y)ds(y), (4.10) and M B 0 e p • e q = ∂B 0 λI -K * #,B 0 -1 [ν p ](y)y q ds (y) (4.11)

  y) = β (ω) k 0 ∂ ∂ν (v + p (y) + y) for y ∈ ∂C, v + p (y)v - p (y) = 0 for y ∈ ∂B, k 0 ∂ ∂ν (v + p (y) + y) = k 1 ∂ ∂ν (v - p (y) + y)for y ∈ ∂B.

  ∂Y∂v p ∂ν (y)y q ds(y)= ∂Y ∂S #,B [ϕ 1,p ] ∂ν (y)y q ds(y) -∂Y ∂D #,C [ϕ 2,p ] ∂ν (y)y q ds(y)and∂Y ∂S #,B [ϕ 1,p ] ∂ν (y)y q ds(y) = ∂B ∂S #,B [ϕ 1,p ] ∂ν + (y)y q ds(y) -∂B ∂S #,B [ϕ 1,p ] ∂ν -(y)y q ds(y) = ∂B y q ϕ 1,p (y)ds(y).

2 Figure 4 . 3 :

 243 Figure 4.3: A domain presenting a symmetry. In this case, the anisotropy direction is frequency independent.

Following [ 6 ] 2 blue arrows for eigenvectors at frequency 10 4

 624 , we use the following values: • The size of cells: 50µm; • Ratio between membrane thickness and size of a cell: 0.7 × 10 -3 ; • Medium conductivity: 0.5 S/m; • Membrane conductivity: 10 -8 S/m; 74 CHAPTER 4. SPECTROGRAPHY OF CELL CULTURES Hz red arrows for eigenvectors at frequency 10 9 Hz (a) (b)

Figure 4 . 4 :

 44 Figure 4.4: (a) shows voltage map with current flows for each y 1 -and y 2 -direction current at 10 4 and 10 9 Hz. (b) shows eigenvectors of the effective conductivity. Blue arrows represent eigenvectors at frequency ω/2π = 10 4 Hz while red arrows are representing eigenvectors at frequency ω/2π = 10 9 Hz.

  and 4.7, in the case θ B = 0 and θ C = 0, and ρ B = ρ C = 0.1.

  real( β(ω)) imag( β(ω))

Figure 4 . 5 :

 45 Figure 4.5: Values of β(ω) for ω/2π ∈ 10 4 ; 10 9 .

Figure 4 . 6 :

 46 Figure 4.6: Real part of the effective conductivity.

Figure 4 . 7 :

 47 Figure 4.7: Imaginary part of the effective conductivity.

2π 0 b 2 0 a 2 2 a 2

 20222 cos 2 (t) √ b 2 cos 2 (t)+a 2 sin 2 (t) dt 2π sin 2 (t) √ b 2 cos 2 (t)+a 2 sin 2 (t) cos 2 (t)+sin 2 (t)

Figure 4 . 8 :

 48 Figure 4.8: Reconstruction of r when there is no inclusion B.

Figure 4 . 9 :

 49 Figure 4.9: Reconstruction of r when there is an inclusion B with ρ B = 0.1.

Figure 4 . 10 :

 410 Figure 4.10: Reconstruction of the orientation of the inclusions B and C.

Table

5 c 0

 0 10 3 ν 0 ≤ 3.7 • 10 12 .

  |E(x)| ≤ 4 • 10 12 , 0 ≤ x ≤ 5 • 10 -3 m.To summarize the above derivation, we have shown that the absolute error E(x) is bounded by|E(x)| ≤ 4 • 10 12 ,x ≥ 0. (A.1)

  86 APPENDIX B. SPECTRUM OF INTEGRAL OPERATORSL 2 (∂C), S #,C [ϕ], ϕ L 2 = -n∈Z d \{0}∂C ∂C e 2iπn•(x-y) 4π 2 |n| 2 ϕ (x) ϕ (y) ds (x) dS (y)

2 I

 2 -1/2 (∂C) such that λI -K * #,C [f ] = 0.Let us writeλI -K * #,C [f ] = λ --K * #,C [f ], 1 L 2 = 0, we have f, 1 L 2 = 0. Let u = S #,C [f ] ∈ H 1 (Y\∂C). Let A = C |∇u(x)| 2 dx and B = Y\C |∇u(x)| 2 dx.

2 I

 2 + K * #,C [f ], S #,C [f ]

L 2

 2

  Schwartz function in the variable z, we have R 2 gac Θ (y)dy = 0. Thus, by (2.5) the expected value may be easily computed as

and (1.15), since gac Θ is a derivative of 32 CHAPTER 2. BLOOD FLOW IMAGING a

Table 4 .

 4 1: Reconstructed values of ρ C with anisotropy ratio of 2.

	Values of ρ C	0.01	0.02	0.03	0.05	0.1	0.2	0.3
	Reconstructed 0.0098 0.0196 0.0294 0.0491 0.0981 0.1963 0.2945

Remerciements APPENDIX B. SPECTRUM OF INTEGRAL OPERATORS

Bibliography

Appendix A

The Justification of the Approximation of the PSF This appendix is devoted to the formal justification of the PSF approximation (1.8) which was obtained by truncating the Taylor expansion of w θ ± at the first order: we shall show here that the error caused by this truncation is small. For simplicity, we shall consider only the case when z = z ′ and θ = 0: the general case may be tackled in a similar way. Without loss of generality, we may set x ′ = 0 and suppose x ≥ 0. We also suppose that we are not too close to the detectors, namely z ≥ 10 -2 m. Moreover, in order to be able to be quantitative, we consider the particular case when F = 0.4 and τ = 1.

The expression of the PSF that we want to approximate is (see (1.7))

where w ± (x) is given by

(Note that, for simplicity of notation, we have removed the dependence of w on θ and z.) An immediate calculation shows that

Hence, there exists

Appendix B

Spectrum of some periodic integral operators

We extend these results to the case of periodic Green's function.

Theorem 11. For any λ > 0, the operator λI + L #,C :

Proof. We first show that the operator L #,C is a Fredholm operator. Note that, L #,C = L C + R where R is an integral operator with a smooth kernel and is therefore compact. Moreover, since L C has a dimension 1 kernel and image, it is a Fredholm operator. Therefore, L #,C is Fredholm. Now we show that L #,C is positive semi-definite, and the result will follow from the Fredholm alternative. Since

for any ϕ, ψ ∈ H 1/2 (∂C), we just have to show that S #,C is negative semidefinite. From the expression (4.12) for G # , we compute, for any ϕ ∈

We have therefore a contradiction : we have |β| ≤ 1 2 since A, B ≥ 0. Therefore, β = -1 2 which implies that B = 0. Therefore, u is constant in R d \ ∪ n∈Z d {C + n}. Since u is continuous across ∂C, u is harmonic on C and is constant on ∂C, and by uniqueness of the Dirichlet problem on C, u is constant on C. Therefore,

which is a contradiction.

Résumé

Les différentes modalités d'imagerie par ondes présentent chacune des limitations en termes de résolution ou de contraste. Dans ce travail, nous modélisons l'imagerie ultrasonore ultrarapide et présentons des méthodes de reconstruction qui améliorent la précision de l'imagerie ultrasonore