Agathe Marmin

Stefan Catheline

Amir Nahas

Simon Chatelin

Gabrielle Laloy-Borgna

Sybille Facca

Sylvain Gioux

Manuel Flury

Nahas Amir

Poster

Emmanuel Martins

Emmanuel Martins Seromenho

Optique Nice

Nadia Bahlouli

Stephane Perrin

Poster -Digital Holography based Full Field Quantitative elastography

Keywords: élastographie, holographie numérique, corrélation de bruit, tissus biologiques elastography, digital holography, noise correlation, biological tissues Contents

en français

Cette thèse de doctorat porte sur le développement d'un système d'élastographie quantitative plein-champ pour l'imagerie de tissus biologiques. Un montage basé sur l'holographie numérique a été conçu afin d'accéder à des images de phase plein-champ à chaque image caméra. A partir de ces images de phase, deux approches principales ont été utilisées pour générer des images de rigidité : l'élastographie impulsionnelle et l'élastographie par corrélation de bruit. Ces approches ont été d'abord testées sur des images de simulation en différences finies mais aussi sur des échantillons polymères aux propriétés optiques et mécaniques contrôlées ainsi que sur des tissus biologiques.

-Chapter 1 -

General Introduction

After brief definitions of the terms in the title of this thesis, its context and objective are introduced in the following sections as well as the organisation of the manuscript and the scientific communications presented throughout this thesis. § 1.1 Definitions This thesis presents a method and experimental set-up based on a digital holography imaging system to perform quantitative elastography. Digital holography is a type of imaging based on light interferences that can be utilized to retrieve displacement images with high-sensitivity. Elastography designates the imaging of mechanical properties, typically the stiffness. Elastography is based on the recording of a medium's response to a mechanical constraint as the displacements within the sample after it was exposed to a force are characteristic of its mechanical properties. The aim of this manuscript is to demonstrate that the combination of an elastography approach with a full-field interferometric system yields quantitative stiffness images with promising results for biological tissues imaging. § 1.

Context and motivations

Since the beginning of medicine, palpation has been used to evaluate the mechanical properties of biological tissues in order to diagnose various pathologies and diseases. Structural abnormalities in tissues such as tumors or fibrosis usually present a higher stiffness than healthy tissues, which can be detected by touch. Palpation is however qualitative, limited to areas accessible by touch and lesions of significant size. Over the past 40 years, elastography technologies have been developed to overcome these limitations, providing physicians with quantitative stiffness images. These techniques were initially based on ultrasound or magnetic resonance imaging and it was not until the late 1990s that elastography was first introduced in optics. Optical systems allow very sensitive and non-contact measurements and recent works showcase the potential of these systems for clinical translation. In particular, the penetration of light in ocular tissue makes optical elastography an excellent candidate for applications in ophthalmology, and the noninvasive, non-contact nature of optical imaging is particularly well suited for intraoperative use.

1 Optical imaging of biological tissues has been widely developed for biology notably through the field of microscopy. Yet, optical imaging has only been recently introduced to medical settings compared to common imaging modalities such as scanners, MRI and ultrasounds. Optics does indeed present significant benefits for biomedical tissues imaging: contactless imaging, high resolution, compact systems, the absence of ionizing radiations and new contrasts.

Most approaches developed in optical elastography are based on the generation of a mechanical impulse in the tissues. The principle is then to film progressive mechanical waves in the field of view in order to calculate locally the stiffness. This classical approach usually requires a controlled source of shear waves in the medium and high imaging rates to track the rapid propagation of the impulse. These considerations reduce the possible clinical applications of this approach. The use of diffuse-field and noise-correlation methods inspired from seismology recently overcame this limitation. Noise-correlation is particularly promising for applications in biological tissue imaging since noise is naturally present in vivo. Noise-correlation could therefore be performed without any external source of mechanical waves, making elastography more easily adaptable to the operating room. In this thesis, an adaptation of this passive method to a full-field interferometric setup is presented. § 1.3 Objectives of the thesis Elastography contrast holds a strong interest for a wide range of applications from diagnostics to surgical imaging. Commercial systems have notably already shown promising contrasts in cornea imaging and liver fibrosis assessments but the elastography systems and techniques are still being developed to improve the image quality and widen the range of applications. This thesis falls within this development and aims at developing elastography contrast for biomedical applications based on an optical imaging system to provide contactless images in real-time.

In this thesis, different elastography approaches are explored, mainly classical transient elastography and noise-correlation elastography. The first objective of this thesis is the design of the optical imaging system and the proof of concept of the noise-correlation approach on the experimental images. The proposed experimental set-up aims at providing contactless quantitative full-field elastography images without the limitations of classical approaches. After establishing the proof of concept on manufactured samples and biological tissues, the method is then refined to also provide in-plane anisotropy images.

Another objective of this work is the design of an experimental set-up adapted to perform different elastography methods, notably classical approaches for control measurements. This versatility notably allows for the exploration of different types of mechanical excitation sources. During this thesis, the noise-correlation algorithm is also generalized to spatially coherent mechanical sources. Classical transient elastography movies are processed successfully with improved image quality compared to conventional approaches. With the evolution of technologies and science, new imaging modalities have emerged for biomedical applications. Physicians and biologists are now provided with a number of systems providing images with a range of contrasts over varied field of views. This Chapter introduces the main modalities for imaging biological tissues and the contrasts they possibly provide, whether it is functional or anatomical imaging. This Chapter aims at providing a basic introduction to existing imaging systems in order to place the work of this thesis in a broader context.

The main modalities used in medical diagnosis are presented in Fig. 2.1 and compared in Fig. 2.2. The presented medical imaging systems often combine a type of probing wave (magnetic resonance, ultrasounds, visible light, x-rays or gamma-rays) with an imaging contrast (phase, amplitude, tomographic, fluorescence, Doppler or spectral imaging). In the first section, the general behavior of waves in complex media is first introduced to lay the foundation for understanding most imaging systems. A brief description of the basic physics behind the contrast of the most common medical imaging systems such as X-rays, ultrasounds, Magnetic Resonance Imaging (MRI) and nuclear imaging is then presented. The different medical imaging modalities are introduced here depending on the type of wave used for tissue imaging. The main characteristics, biomedical applications and limitations are discussed for each imaging modality.

The following sections focus on optical imaging. Optics are mainly used in medicine for specific applications such as microscopes for anatomopathology, optical coherence tomography (OCT) for ophthalmology or cameras for endoscopy. However, the use of optics has been recently emerging as a promising tool for biomedical imaging, providing new contrasts as well as fast high-resolution images. The different optical imaging modalities for biological tissues are introduced. Optical imaging techniques based on polarization, absorption or scattering contrast are the subject of a first section. Since the focus of this thesis is on holographic imaging, interferometric imaging techniques are specifically presented in another section.

At last, the final section presents imaging techniques based on the combination of acoustic and optical imaging, whether it is combining both waves or adapting techniques from one domain to the other.

Waves in complex media

Waves are defined here as a disturbance propagating in a regular and organized way. Two types of waves are often used for imaging: electromagnetic waves and mechanical waves. Electromagnetic waves classically designate oscillating magnetic and electric fields produced by the vibration of charged particles whereas mechanical waves designate the propagation of energy through particle-to-particle interaction, often referred to as a pressure or elastic wave. Mechanical waves (such as sound, water or seismic waves) consequently require a medium to propagate but electromagnetic waves (such as light, radio-frequencies or x-rays) can transfer energy through vacuum. As a reference, the propagating speed of electromagnetic radiations in vacuum is c = 2.998 × 10 8 m/s and the speed of sound in water is around c = 1.480m/s × 10 3 . Yet, although mechanical and electromagnetic waves are very different in nature, they share common characteristics and behaviors that are described by classical wave theory and defined in this section.

General descriptions and characteristics of waves

General definitions

Waves characteristics vary tremendously with their type and frequency. They are often grouped under subcategories as shown for the acoustic and electromagnetic spectra in Fig. 2.1. Other adjectives can be used to describe a wave's temporal profile (pulsed, gated, periodic...) or its propagation direction in respect to the oscillations (transverse, longitudinal, surface...).

The wavelength (λ), wave-number (k), temporal period (T), frequency (f or ν), speed (c) and angular frequency (ω) are used here to characterize the waves. For electromagnetic waves, the energy (E) can also be used. The following formula present the links between these quantities :

           k = 2π λ f = 1 T ω = 2πf c = f λ E = hc λ = f h (2.1)
where h denotes Plank's constant (≈ 6, 626 × 10 -34 J).

The phase and amplitude of a wave characterize the local phase and amplitude of the oscillations. The intensity of a wave is the energy arriving per unit time per unit area averaged over a period of oscillation, usually considered equal to the square of the amplitude.

Wave equation

In classical wave theory, the equation of propagation is a second-order partial differential equation and can be written:

∆ψ - 1 c 2 ∂ 2 ψ ∂t 2 = 0 (2.2)
with ψ the wave function and ∆ = ∇ 2 the Laplace operator.

This basic formula describes wave propagation in the absence of a source in the case of a lossless propagation medium. The one-dimensional version of this wave equation was first derived by d'Alembert in 1747 for the case of the vibrating string. For an homogeneous medium, the general solution in one dimension (x) is a propagative monochromatic plane wave: ψ(x, t) = f (x -ct) + g(x + ct) (2.3) with f and g arbitrary twice-differentiable functions, corresponding to forward and backward propagation along the (0x) axis.

Traditionally in wave-propagation studies, it is advantageous for the description of traveling waves to take the first function f corresponding to a forward propagating wave. The general monochromatic wave-function ψ(r, t) is commonly written:

ψ(r, t) = A(r)cos(ωt -k • r + ϕ) (2.4)
with ω the wave frequency, k the wave vector and A the wave's amplitude and ϕ the phase angle.

A complex-number representation is often convenient in theoretical studies. It replaces the amplitude and phase by a single complex number and condenses the writing of mathematical relations. The complex representation rely on Euler's theorem :

e iϕ ≡ cos ϕ + i sin ϕ, (2.5)
Consequently, the general wave function ψ(r, t) can be represented by the following complex expression: ψ(r, t) = A(r)e i(ωt-k•r+ϕ) (2.6) where the real part of this complex function corresponds to the expression in Eq. 2.4.

Wave and matter interactions

Waves interactions with matter can be described through four main phenomena: absorption, reflection, refraction and scattering, see Fig. 2.3. A general description of these behaviors is suggested in the following subsections. In simple cases such as a single interface between two semi-infinite medium, geometric definitions or the wave equation are sufficient to describe wave and matter interactions. However, in complex media, such as biological tissues, averages quantities are defined, notably for describing absorption and scattering in a medium. Both approaches are used here to describe the interactions of waves with matter.

Reflection and refraction

Refraction is described by a change in direction of an incident wave after being transmitted through an interface between media with different properties. Reflection is when an incident beam is partially reflected at an interface. These behaviors emerge from a mismatch between the wave speeds of the media.

The refractive index in a homogeneous medium can be defined as a wavelengthdependent complex quantity ñ(λ): ñ(λ) = n(λ) -iα(λ)

(2.7)

The imaginary part α encompasses attenuation through absorption, which is introduced in the following subsection. The real refractive index (conventionally called refractive index) of a medium n characterizes the phase velocity c medium of the wave in the medium and its wavelength λ medium :

c medium (λ) = c n(λ) λ medium = λ n(λ) (2.8)
with c and λ the speed and wavelength in vacuum for electromagnetic radiations or the speed and wavelength in water for acoustic waves.

From Eq. 2.8, it can be concluded that wave propagation in tissues affects the wavelength and speed of the incident wave but has no impact on its frequency.

For acoustic waves, the acoustic impedance is most often used to describe behaviors when there is a mismatch between the properties of two media. Acoustic impedance Z is defined by the medium's density ρ and the acoustic wave velocity c us [1]:

Z = ρ c us (2.9)
Considering two homogeneous semi-infinite media with different refractive indexes and an incident wave beam with its wavelength small compared to the interface or its curvature, reflection and refraction occur (see Fig. 2.4). A description of the phenomena of refraction and reflection at an interface can be retrieved by using the wave equation 2.2 with a planar boundary:

c 2 sin(θ 1 ) = c 1 sin(θ 2 ) θ 3 = θ 1 (2.10)
with c 1 and c 2 the wave speed of the two media, θ 1 , θ 2 and θ 3 , the angles between the normal to the interface and the incident, refracted and reflected waves, respectively.

From this equation, in the case where the wave speed is smaller in the incident medium, a critical angle can be observed at θ c = arcsin(c 1 /c 2 ) from which the wave is fully reflected. This equation if often formulated differently for electromagnetic waves and acoustic waves. The real refractive index is used for electromagnetic waves and the following equation is known as Snell-Descartes' law:

n 1 sin(θ 1 ) = n 2 sin(θ 2 )
(2.11)

For a normal incidence on a plane boundary, the fraction of energy transmitted (T) and reflected (R) by the interface can be expressed using Fresnel's equations:

T = 4n 1 n 2 (n 1 +n 2 ) 2 R = 1 -T = (n 1 -n 2 ) 2 (n 1 +n 2 ) 2
(2.12) As for acoustic waves, the acoustic impedance is preferred:

Z 2 sin(θ 1 ) = Z 1 sin(θ 2 )
(2.13) Also using the acoustic impedance, the reflected and transmitted energy ratio for an ultrasound beam with normal incidence at an acoustic boundary between two mediums (with respective impedances Z 1 and Z 2 ) can be quantified similarly to Eq. 2.12:

T = 4Z 1 Z 2 (Z 1 +Z 2 ) 2 R = (Z 1 -Z 2 ) 2 (Z 1 +Z 2 ) 2 (2.14)
When the incidence of the plane wave is not normal, the reflection and transmission coefficients are different from Eq. 2.12 and 2.14 and their expressions involve the wave's polarization state and can be found in literature for electromagnetic waves [2] and for acoustic waves [3].

Absorption

Absorption of the wave by a medium can be explained by the conversion of the incident energy from the wave to the medium, for example through thermal diffusion. Absorption can be described by the presence of absorbing particles. Yet, at a bigger scale, a statistical approach is preferred to describe this phenomenon. The absorption coefficient µ a corresponds physically to the intensity variation (dI) along an elementary distance of the propagation axis (dz):

dI = -µ a Idz (2.15)
µ a is medium and wavelength dependent and is expressed in cm -1 . When integrating Eq. 2.15 on a width z of the medium, Beer-Lambert's attenuation's law is obtained:

I = I 0 exp(-µ a z) (2.16)
where I 0 is the initial intensity of the beam and I the transmitted intensity through the medium.

The absorption length l a represents the average distance traveled by the wave before an absorption event. It can be expressed as:

l a = 1 µ a (2.17)
For a medium of a defined thickness, the transmittance T and absorbance A are also often used to characterize the absorption of a collimated beam. These coefficients are defined as:

T = I 0 I A = log 10 ( I 0 I ) (2.18)
with I the intensity after propagation through the medium and I 0 the initial intensity.

Scattering

Scattering can occur either when a medium contains particles with a different refractive index or when the medium has a continuous but fluctuating refractive index. The particle or fluctuation then acts as a secondary wave source with a specific intensity distribution. This intensity distribution depends on the incident wavelength λ, the local index gradient in the case of a fluctuation or the particle's size a, its refractive index n s as well as its surrounding medium's n m . Scattering occurs in electromagnetic and mechanical waves [4].

In the case of spherical scattering particles, scattering regimes of electromagnetic and acoustic [3] waves can be described by three main processes:

• When the particle's size is small compared to the wavelength (usually <λ/10), this scattering regime can be characterized using Mie theory (theory of scattering by spherical objects) under Rayleigh's limit. This scattering is often described as Rayleigh elastic scattering. It can be characterized by its effective cross-section σ s :

dσ s dΩ = 8π 2 n 2 m ( n 2 s -n 2 m n 2 s + 2n 2 m
) 2 a 2 λ 4 (1 + cos 2 (θ))

(2. 19) with θ the angle between incident wave's direction and the direction of diffusion considered and dΩ the elementary solid angle.

• Mie scattering designates wave scattering described by Mie theory in the cases where neither Rayleigh nor geometrical approximations are valid. In this case, the scattering particle can be described as a group of dipoles excited by the incident wave with different phases. The intensity scattered in different directions consequently reflects the interferences between the different secondary waves emitted by the oscillating dipoles.

• Geometrical scattering occurs for particles that are bigger than the incident wavelength. In this case, geometrical approaches are sufficient to describe scattering.

Contrary to absorption and as described earlier, scattering is a directional and often anisotropic process. To define anisotropic scattering, a scattering phase function p(s, s ′ ) is necessary to describe the fraction of energy from incident light on the scatterer in the direction s that is scattered in the direction s ′ . However, under the assumption an uniform distribution of density of the scatterers ρ, a medium can be described by its scattering cross section σ s independently of the direction of the wave: 20) with P scatt the power scattered per unit area of a uniform plane wave of intensity I 0 and µ s the scattering coefficient. Consequently, scattering in biological tissues can be characterized similarly to absorption, using the scattering coefficient µ s , which corresponds to the intensity variation in a collimated beam going through an elementary distance without being scattered :

σ s = P scatt I 0 = µ s ρ (2.
I = I 0 exp(-µ s z) (2.21)
with I 0 the initial intensity and I the remaining intensity in the collimated beam after a distance z.

The scattering mean free path l s is also defined and corresponds to the average distance a wave travels in the medium between two scattering events: Three scattering regimes are distinguished depending on the considered medium: ballistic regime, single scattering and multiple scattering. These different behaviors are described Fig. 2.6. Ballistic regime is describes wave propagation through a medium in a straight line without too many scattering events (l c << l s ). Single scattering occurs when one scattering event is happening on average (l c ≈ l s ). At last, multiple scattering describes media-wave interactions when each wave undergoes multiple scattering events (l c >> l s ). § 2.

l s = 1 µ s (2.22)

Common Medical Imaging Modalities

After the previously introduced basic notions of wave and matter interactions, the most common imaging modalities currently used in medical settings for biological tissue imaging are described in this section. The basic imaging principle, contrast and principal medical use and advantages and limitations are described for each system. X-rays, ultrasounds, MRI and nuclear imaging are discussed here.

X-rays imaging

Basic principles

The discovery of the X-rays in 1895 by W. C. Röntgen led to the development of a whole field in medicine: radiology. X-rays are considered the oldest imaging modality in medicine and allowed for the famous first image of the bones in an human hand in-vivo (Fig. 2.7). X-rays are electromagnetic waves with a photon energy range generally from 10 to 100 keV for diagnostic imaging [5].

X-rays can be produced by either atomic inner electronic transition or electron deceleration (Bremsstrahlung). This is commonly achieved by heating an cathode material to produce accelerated electrons.

The collision of these electrons with an anode material creates x-rays.

The created beam can then be directed to the region of interest.

After penetrating the tissues, the remaining xray beam can be used to retrieve structural information.

In the considered frequency range, two main interactions of x-rays with biological tissues are photoelectric effect and Compton scattering, illustrated in Fig. 2.8. Compton scattering is the inelastic scattering of photons by atoms in the matter. The photoelectric effect consists in the absorption of the energy of a photon by an electron. Both these processes, schematized in Fig. 2.9, cause the attenuation of the x-ray beam but are also responsible for the ionization of atoms due to the loss of electrons and the liberation of electrons with kinetic energy.

Medical use and accessible contrasts

The use of x-rays in healthcare can be separated in attenuation imaging, phase imaging and spectroscopy.

• Spectroscopy techniques consist in imaging the fluorescence spectrum of biological tissues after x-ray exposition. Fluorescence from x-ray exposition is caused by electronic transitions between atomic shells after ionization (loss of an electron).

The fluorescence wavelengths consequently correspond to the difference between electron binding energies between each shell. Different atoms consequently have different fluorescence wavelengths. By analyzing the fluorescence spectrum, the presence of certain substances such as lead can be detected and quantified. X-ray spectroscopes are used for compound analysis in medicine but do not fall in the domain of imaging. X-ray fluorescence imaging does exist but is mainly used for security scanners.

• Attenuation imaging is the most common in medical x-rays. The attenuation of an x-ray beam of wavelength λ can be described by Beer-Lambert's law Eq 2. 16 with the mass attenuation coefficient expressed as:

µ a = ρλ 4πα (2.23)
with α the complex part of the refractive index of the medium and ρ its mass density. This coefficient is proportional to the effective atomic number (mean number of protons in atoms of a compound).

During an x-ray scanner, each ray of the collimated x-ray beam propagates through different tissues and undergoes a certain amount of attenuation. The remaining beam is then imaged. The contrast of attenuation imaging consequently comes from the differences in atomic composition, density and thickness of biological tissues. The main contrast occurs between bones and soft tissues. X-ray scans are consequently often used for diagnostics for imaging bone fractures, tooth problems as well as breast tumors or lung diseases.

The introduction of surgical tools or dyes can also provide useful contrasts. Continuous x-ray imaging is defined as real-time x-ray imaging. It is performed to track surgical tools or bone position during radiography-guided surgery or with contrast agent (usually iodine or barium-based solutions) for angiography.

At first, radiography was made using photosensitive films. More recently, x-rays acquisitions have been digitized using the combination of a scintillator film and a visible camera or more rarely using directly an x-ray camera. This allowed for tomographic x-ray imaging starting in the 70s, which is more commonly named Computed Tomography (CT) or Computed Assisted Tomography (CAT) scan [6]. This three-dimensional reconstruction allowed for three-dimensional contrast at the cost of longer acquisition times and therefore increased radiation dose. Indeed, excessive exposition to ionizing radiations such as x-rays can cause DNA damage in the tissues. This is particularly relevant for CT-scans where the patient is exposed for up to 100 times the dose for an equivalent flat-panel image [7].

• X-rays can also be characterized by their phase. Due to the need to develop coherent sources and phase detection systems, phase contrast x-ray imaging only emerged in the 1990s. Phase imaging consists in recording the phase shift of the radiation after going through tissues. Contrary to attenuation x-rays, the contrast here comes from the real part of the refractive index of the medium. Phase-contrast images are relying on the electron density of the medium, which presents higher contrasts within soft tissues [8]. Because x-rays phase imaging requires higher energy radiations, which are less absorbed by the human body, it is a promising approach for dose reduction.

Principal advantages and limitations

X-rays were the first imaging modality able to provide physicians with in-depth and tomographic images. This modality of imaging is relatively cheap and does not require contact with the patient. It can be a great tool for cardiovascular or bone surgery. Xrays are currently mainly used clinically for diagnostics in bone, teeth, calcification or foreign matter imaging as well as for guided surgery.

The contrast x-rays offer is however still very limited in soft biological tissues. In flat panel images, regions of high contrast (like bones) can block significant diagnostic data. Recent approaches based on phase imaging have shown promising results to improved contrasts in routine breast and lung imaging. The resolution of clinical x-ray images is typically in the order of magnitude of a millimeter. This resolution is usually limited by the diffraction from the biological tissues and the scintillator used for imaging but also by possible movements from the subject. X-ray imaging is based on the use of ionizing radiations, which is a main concern for their routine use. Although x-ray imaging has been used for more than a century, the technique is still evolving to decrease the radiation dose, increase resolution and contrast but also to adapt it to endoscopic systems [10].

Ultrasound imaging

Basic principles

The use of ultrasounds for medical imaging dates from the 1950s. Unlike x-rays and light, ultrasounds are not electromagnetic waves but acoustic waves: they originate from organized movement of the matter. Acoustic waves can be divided into three frequency domains (see Fig. 2.1) : infrasound, audible sound and ultrasound. In biological tissues, sound waves travel at around 1500 m/s for soft tissues and 3500 to 5000 m/s for bone [11]. Different medium respond differently to sonic waves depending on how the particles of the medium resist mechanical vibrations. This property is the medium's acoustic impedance, defined previously in Eq. 2.9 as a function of the medium's density ρ and the acoustic wave velocity.

Medical ultrasound imaging is based on the emission and recording of an ultrasound signal through biological tissues. Piezoelectric materials are used to convert electrical voltage into mechanical pressure and act as a generator but also as a detector for ultrasound probes. A sonic wave is sent through the tissues and the echoes from partial reflections are then recorded. The time-delay and amplitude of the recorded signal inform on the tissues' impedance mismatch, which corresponds to differences in the tissue's densities and/or ultrasonic wave speeds.

Contrast in medical ultrasound imaging emerges from the back-scattered signal. Because absorption of ultrasounds in biological tissues is relatively low at lower frequencies, ultrasounds allow for in-depth imaging.

Medical use and accessible contrasts

Ultrasound systems are often employed in obstetrics, cardiology, and radiology for diagnostics, as it can provide in-depth anatomical images as well as angiography for diagnostics.

Ultrasound imaging methods can be divided into two categories: continuous waves (CW) and pulsed waves (PW). The CW mode relies on two sets of ultrasound transducers, one emitting while the other records continuously. On the other hand, the PW mode requires only one set of actuators emitting and recording successively. Ultrasound systems often present three principal acquisition modes: Amplitude-mode (A-mode), Brightness-mode (B-mode) and Motion-mode (M-mode), see Fig. 2.11. Amode consists in directly displaying the intensity of the recorded back-scattered ultrasound signal in function of time. A spike on this curve corresponds to the presence of a boundary and its intensity to the impedance mismatch. The B-mode displays the same information as a line of pixels. The spikes from the recorded A-mode signal are presented as dots, which brightness is proportional to the intensity of the echoes. The B-mode display allows the user to sweep the ultrasound beam to create images. At last, M-mode is used to generate images of moving objects: successive screenings allows to track the movement of acoustic interfaces in the direction of the ultrasound propagation.

Doppler.

Ultrasound imaging based on the Doppler effect allows for vascular imaging and blood flow quantitative assessment. The Doppler effect is based on the relative movement of the source or the receptor during wave recording. When using ultrasounds on the human body, ultrasonic waves are reflected by moving tissues or fluids. These particles, usually red blood cells, can consequently be considered as secondary ultrasonic sources. If they are moving towards or away from the detector, a frequency shift can be recorded between the emitted pulse frequency and the received echo. This frequency shift (f d ) is directly related to the speed of the receiver relative to the source. For ultrasounds Doppler, we note V the velocity of red blood cells and f t the transmitted ultrasound frequency from the probe. The Doppler shift can be expressed as [12]:

f d = 2f t V cos(θ) c (2.24)
with c the speed of sound in the tissue, θ the angle between the propagation direction of the ultrasonic beam and the propagation direction of the reflecting particles.

Figure 2.12: Schematic of the Doppler effect for blood flow imaging [12].

In the Doppler mode for ultrasounds, a tilt of the probe is usually necessary to optimize the Doppler signal. Doppler ultrasounds are usually performed using the pulsed mode so that the distance between the probe and the Doppler source can be calculated and blood flow images computed. Depending on the final image and the analyzed signal, different Doppler techniques can be found in the literature such as spectral Doppler, power Doppler or color Doppler. An example is provided Fig. 2.13 where the blood flow in the jugular vein and the carotid artery has been acquired using a Doppler ultrasound system.

Principal advantages and limitations

Ultrasound imaging is a well-established clinical tool for the morphological assessment of soft tissues which is employed frequently. Ultrasound imaging is cheap, easy to use and require little technical background. It provides in-depth images in real-time. It is widely use for anatomical images or angiography for diagnostics. Ultrasound systems are considered safe as they require not ionizing radiation. If ultrasounds can cause heating in biological tissues due to absorption and cavitation, the frequencies and amplitudes used for medical imaging are considered extremely safe.

Ultrasound medical imaging contrast emerges from partial reflections. However, great impedance mismatches (such as between air, soft tissues and bones) cause total reflection and prevent further in-depth imaging. Ultrasound systems consequently require coupling gel between the skin's surface and the probe and cannot image past bones. Ultrasound systems also suffer from a low specificity. The use of contrast agent in the form of micro-bubbles has been recently developed and has shown great potential for improving contrast [14].

Depth, resolution and field of view trade-offs. The main limitation for in-depth imaging using ultrasounds is the absorption by biological tissues. Absorption increases with the ultrasound frequency, the medium's viscosity (the measurement of friction between particles) and the medium's relaxation time (the time taken by particles submitted to an ultrasound vibration to revert to their original position). Ultrasound frequencies for medical imaging range typically from 2 to 40 MHz [11]. The axial and lateral resolution (δ z and δ x ) of ultrasounds imaging are in the order of magnitude of the chosen imaging wavelength and can be approximated by the following formula [1]:

δz ≥ λ/2 δx ∼ 3λ (2.25)
For ultrasonic imaging, the higher the frequency, the higher the resolution but the lower the penetration depth. As for examples: fetus or liver imaging is usually done around 3.5MHz to achieve 15cm penetration depth and around 1mm spatial resolution while skin imaging is usually at 20MHz for 1cm penetration depth and up to 40µm axial resolution. This trade-off can be limiting for ultrasound imaging in patients with ascites or obesity.

Magnetic Resonance Imaging Techniques

Basic principles

Nuclear Magnetic Resonance (NMR) was first exposed in 1946 by E. M. Purcell and F. Bloch. The 70s saw the rapid development of imaging systems based on the localization of NMR signals [15]. In NMR imaging, the characteristic of certain atom's nuclei (those with appropriate magnetic spins) is exploited to provide images. The contrast emerges consequently from the atomic constitution of the tissues. There are three main biomedical applications for magnetic resonance: spectroscopy, magnetic resonance imaging (MRI) and nuclear magnetic resonance microscopy. MRI is the most widely used systems in medical imaging and is detailed in the rest of this section.

MRI systems consist of a very strong electromagnet capable of producing a magnetic field of a few Tesla over the field of view. This causes the spin of certain atoms (typically hydrogen) to go from a random distribution to precessing around the direction of the magnetic field. The phenomenon is the magnetic polarization of tissues. The frequency of the oscillation of the magnetic field during polarization is named Larmor frequency (f l ). It is a function of the magnetic field's intensity (B 0 ) and the gyromagnetic ratio of the atom γ [16]:

f l = γB 0 (2.26)
Water is strongly present in the human body and its very high gyromagnetic ratio (γ = 43M Hz/T ) makes this element especially interesting for MRI. Consequently, when tissues are exposed to a strong magnetic field, their are magnetized. The magnetization vector (M) is used to measure this phenomenon: it consists in the sum of the contribution of the nuclei magnetic spins of a given voxel. The magnetization vector (M) is not precessing around the magnetic field's axis since the hydrogen protons' are considered to precess out of phase. A radio frequency (RF) pulse of a fixed frequency and intensity with a low magnetic component B 1 << B 0 can then be used to create in-phase precession of the protons, tipping the magnetization vector M orthogonal to the magnetic field's direction and generating oscillating components at a Larmor frequency f 1 = γB 1 . When the RF magnetic field is disabled, the direction the nuclei spins goes back to a random distribution around the axis of B 0 in a process named relaxation.

The speed of relaxation depends on the type of tissue as dense tissues for example limit molecule movements and lengthen the relaxation time. The measurement of the signal during the relaxation yields the contrast in MRI imaging [11]. The 3D images are obtained by using multiples pulse sequences and a gradient magnetic field. Different sequences for the acquisitions will yield different images.

Medical use and accessible contrasts

The principal MRI sequences are described here and the corresponding differences in contrast are illustrated Fig. 2.15.

T1-weighting, T2-weighting and proton density weighting. Spin-echo sequences are the most common sequences used for structural MRI. Three main contrasts can be used for spin-echo MRI acquisitions depending on the repetition time (time between each radio-frequency pulse) and the echo time (time between the radio-frequency pulse and the acquisition):

• T1-weighting consists in exploiting the time of longitudinal relaxation. T1 is the time constant characterizing the exponential decay of the longitudinal magnetization vector (M z ).

||M z (t)|| = ||M z (t = 0)||(1 -exp(-t/T 1)) (2.27)
This type of acquisition shows higher level of signal for fat or gadolinium-based contrast agents.

• T2-weighting consists in exploiting the time of transversal relaxation. T2 characterizes the decay of M xy (t):

||M xy (t)|| = ||M xy (t = 0)||exp(-t/T 2) (2.28)
This type of acquisition shows higher level of signal for water content.

• Proton density weighting consists in providing images of the number of mobile hydrogen protons per unit of volume Functional MRI. MRI is the only modality allowing for non-invasive brain imaging due to the ability of radio-frequencies to penetrate past the high density cranial cavity. Moreover, the emergence of functional MRI in the 1990s has allowed for functional brain imaging [15], Figure 2.15: Images of different modalities used for brain imaging. Functional MRI (fMRI) shows regions of brain activity due to visual stimulus. Positron Emission Tomography (PET) shows regional glucose use. MRI-T1w is an anatomical image of the brain. CT is used for detecting brain lesions. Diffuse MRI (dMRI) shows an undiffused B0 image. MRI-T2w is used for detecting white matter lesions. Images from [17]. [18]. Contrarily to structural imaging, which is used to screen tissues in the body, functional imaging aims at visualizing biological processes. Neural activity results in a local hemodynamic response. This can be tracked and image by two principal MRI techniques:

• Blood Oxygen Level Dependent (BOLD) functional MRI records changes in deoxygenated hemoglobin concentrations in the brain. These concentration changes can happen due to unregulated processes in a resting brain or can be consequent to performing a task. While fully oxygenated hemoglobin is diamagnetic, fully deoxygenated hemoglobin is high paramagnetic. Deoxygenated hemoglobin concentration can consequently be recorded using specific MRI sequences.

• Cerebral Blood Flow (CBF) can image neural activity using increases in blow flow. This technique consists in using either the injection of an exogenous contrast agent or a spin-labelled perfusion.

Diffuse Tensor Imaging (DTI). Diffusion tensors are matrices used to model local water diffusion in brain tissues. Indeed, water molecules diffuse preferentially along certain directions. DTI uses this local water diffusion to construct an image of white matter tracts, showing the connections between different regions of brain. For clinical usage, different contrasts like fractional anisotropy, mean diffusivity and apparent diffusion coefficient can be used [19]. Generally, clinical DTI acquisitions have low spatial resolution (2-5 mm). Diffusion weighted imaging assess the ease with which water molecules move around within a tissue and gives insight into cellularity, cell swelling and edema.

Principal advantages and limitations

MRI systems allow for full-body tomographic images and functional brain imaging with a unique range of contrasts. MRI are used for diagnostics in a variety of cases, from brain conditions to joint abnormalities, including the detection of tumors. Compared to a CAT scan, it has the advantage of providing a better contrast in soft tissues and doesn't require the use of ionizing radiations.

The clinical MRI images often have a moderate spatial resolution (≈mm) and signalto-noise ratio, motion artifacts and geometric distortions related to magnetic susceptibility and magnetic field inhomogenity. Even though MRI can provide a broad range of contrast allowing for 3D morphological images, full-body MRI systems are bulky and very expensive due to the need to generate controlled magnetic fields of high intensity. MRI is also time-consuming due to the long acquisition sequences necessary for tomographic reconstruction. The use of a strong magnetic field also prevents imaging if magnetic objects are present. It may be problematic for patients with artificial joints, rods, plates or screws but also with artificial heart valves or a pacemaker.

Nuclear Imaging

Basic principles

Nuclear imaging emerged in 1925 and first clinical systems were used for diagnostics in the 1950s and 1960s [20]. It is based on the targeting of a specific molecule [21]. Nuclear imaging relies on the combination of radio-pharmaceuticals and specific imaging systems. Radioisotopes are tagged to a tracer and are administrated to the patient. The nuclei of these radioisotopes emit gamma rays that allow them to be traced by using a single photon emission computed tomographic (SPECT) or positron emission tomographic (PET) scanner [22]. Common radiotracers are iodine-solution or fluorodesoxyglucose, which are used for tumor imaging where the carbohydrate metabolism is deregulated.

SPECT.

Gamma rays are created by radioactive decay where the nuclei of unstable isotopes emit highly energetic particles in order to reach a stable state. These radiations are captured by a gamma camera. Gamma cameras usually have the same design as x-ray camera. They consist of a collimator to guide the light to the scintillation crystal, that converts these highly energetic radiations into visible light.

PET.

In the case of positron imaging, the radioactive decay process is slightly different. The unstable isotope emits a positron, which travels until encountering an electron. An annihilation ensues, producing a pair of photons travelling in opposite direction (of energy 511 keV). These photons can be detected and the site of annihilation can be mapped.

Cerenkov radiation imaging.

The Cerenkov effect can be described as the emission of visible photons when a charged particle velocity is greater than the speed of light in a dielectric medium. This effect was first observed at the beginning of the 20th century as a blue glow around radioactive solutions but the theory and first applications for biomedical imaging emerged only recently. In biological tissues, the index of refraction is around 1.36 to 1.40 [23]. Consequently the energy threshold for the Cerenkov effect is around 0.2 MeV, which is below the energy of particles from most radionuclide used for PET and SPECT scans. The Cerenkov effect can consequently be observed from these existing radionuclides. This effect has be imaged in recent in-vivo studies [24], [25]. Because this phenomenon creates few photons and requires high levels of radioactivity, it is a promising tool for monitoring during radiotherapy (when administered doses are already very high) without the need for high-energy scanners such at PET or SPECT.

Principal advantages and limitations

Nuclear imaging is often coupled with CT or sometimes MRI systems to superimpose the signal from the radioisotopes on the anatomical images. The use of tracers allows to track and monitor an organ's function. Nuclear imaging is often used in diagnostics for cancers as well as cardiac, endocrinal and neurological disorders.

Nuclear imaging is very similar to CT scans in terms of its limitations: the resolution of images is of medium quality with good contrasts and since nuclear imaging is based on the use of radioactive tracers that produce ionizing radiations, it should consequently performed sparingly. The acquisition of nuclear images takes time and qualified physicians since it requires the administration of a contrast agent and successive screenings in order to assess organ function over time. § 2.4 Optical imaging Contrary to other imaging methods such as MRI or ultrasounds that saw a rapid development and integration into clinical systems during the XXth century, commercial optical-based imaging prototypes have only been developed in the last decades for medical application. Because of the scattering and heterogeneous nature of biological tissues, light's trajectory is complex and depends on a broad range of parameters. Exploiting this signal consequently requires adapted sources, detection systems and models. In this section, the optical properties of tissues and the related physical quantities are defined as they are the main source of contrast in optical imaging. A range of optical imaging systems for biological tissues are then presented. Since the set-up developed during this thesis in based on interferometry, interferometric imaging systems are not discussed here but are detailed in a following section.

Optical properties of biological tissues

Light can be defined with different levels of abstraction. Light can be described as a scalar but also as an electromagnetic wave or even as a stream of particles called photons (see Fig. 2.17). In order to understand light interactions with tissues, these different descriptions of light are used here. Biological tissues present a broad range of structures of different shapes and sizes which hinders light propagation. When light is directed at a material, it can be either absorbed by the material, be transmitted through it or be reflected back towards the source [26], [27]. These interactions can be described by different quantities defined generally in section 2.2. However, when considering light and matter interactions, other phenomena specific to the electromagnetic and quantum nature of light should be considered. These properties are discussed in the case of biological tissues in the rest of this Chapter 2. Biological Tissues Imaging section.

Refractive index

In biological tissues, the refractive index varies between 1.35 and 1.45. This variation is principally affected by the presence of water.

Absorption

In biological tissues, the absorption coefficient µ a , defined in section 2.2 is between 0.01 and 0.1 cm -1 [28]. Optical wavelength are consequently partially absorbed by tissues .

Because molecules absorb photons when their energy matches internal energy states' transitions, absorption is highly dependant on the light's wavelength. For biological tissues, the absorption spectrum (corresponding to the function µ a (λ)) reflects the tissue's molecular compositions. Some examples are display Fig. 2.18. For short wavelengths, such as UV absorption, the energy corresponds to electronic transitions (DNA, proteins...). Bigger wavelengths in the red or near-IR have been utilized in biomedical optics for therapeutic use or for oxygenation imaging as hemoglobin and oxygenated hemoglobin present different absorption coefficients in this range.

Scattering

Biological tissues are mainly considered as very scattering for optical wavelengths. The scattering coefficient µ s , defined in section 2.2 ranges around 100 cm -1 [28]. Optical wavelength are consequently rapidly scattered by tissues. This is mainly due to the presence of collagen fibers and cell constituents acting as scatterers. Indeed, their size matches the order of magnitude of the wavelengths in the optical spectrum, as illustrated in Fig. 2 Considering electromagnetism's formalism, the electric field, the magnetic field and the propagation direction of the wave are orthogonal to each other. However the orientation of the oscillation of the electromagnetic field in respect to the propagation direction can vary in space and time.

The projection of the electric field in a plane perpendicular to the propagation direction gives the polarization state of the light: polarization is generally represented as elliptical but in some cases, it can be linear or circular, as illustrated Fig.

2.21.

The polarization ellipse can be described by its geometrical parameters: its orientation and eccentricity. The Poincare parameters ψ and χ are used here to characterize the polarization ellipse, see Fig. 2.20.

Polarized light can be described using Jones vector:

J = E s E p (2.29)
with E s and E p the complex amplitudes of the electric field in the two directions orthogonal to the propagation direction.

Light can be polarized fully or partly or even unpolarized. The proportion of polarized electromagnetic waves can be quantified by the degree of polarization p. More generally, Stokes vector describes the polarization state of light with four parameters, taking into account the light intensity (I) and degree of polarization (p): Birefringence defines materials where the refractive index depends on the polarization of the incident light. The action of a medium on light's polarization properties can be described using a Mueller matrix. A Mueller matrix (M ) is simply the 4X4 transfer matrix between an initial state of polarization (with Stokes vector S in ) and the state of polarization after an optical element (with Stokes vector S out ):

S out = M S in (2.31)
Linear birefringence is found in biological tissues such as collagen, muscle fibers or retina. The application of polarimetric measures to biological tissues have consequently shown interesting contrast to investigate the anisotropic optical properties of specific structures such as skin derma, muscle, bones, elastin, myosin fibers as well as the depolarization properties of tissues [31], [32].

Fluorescence contrast

Some biological constituents have the ability to convert incident light in a different frequency. This property can be used for optical imaging. Instead of directly imaging the light forward-scattered or back-scattered by the sample to draw contrast from absorption or scattering properties of tissues, this converted light can be imaged to obtain new contrasts. Fluorescence has notably been instrumental for functional imaging.

Fluorescence describes a phenomenon where light is absorbed by a molecule through the photoelectric effect, the additional energy from the excited state then dissipates partially through the creation of a lesser energy photon. Fluorescence is an interesting contrast for biomedical imaging since it allow to track molecular processes. Additionally, endogenous or exogenous contrast agents can be used to target specific cell mechanism. Main applications of fluorescence imaging include microscopy, where cellular and molecular interactions can be imaged with high specificity, and intraoperative imaging for tumor or blood flow imaging.

Fluorescence occurs through an electronic transition from the excited state to the ground state of a molecule. The process is represented Fig. 2.22 as well as other processes responsible for energy dissipation after light absorption. Fluorescence in the visible and UV range occurs between electronic state transitions while infrared fluorescence emerges from vibrational transitions. The absorption and emission spectra of fluorochromes are usually symmetrical and the distance between their peaks is called the Stokes' shift and is illustrated in Fig . 2.23. A high Stokes shift causes less signal loss due to spectral filtering in the experimental acquisition setup and is a crucial parameter in fluorescence.

The basic optical set-up used for fluorescence imaging is schematized in Fig. 2.24. By using filters and a dichroic mirror, the excitation and the emission light can be separated to efficiently image fluorescence in a sample.

Γ f = νσ a I hν (2.32)
Consequently, good fluorochromes for optical imaging have a good absorption crosssection (σ a ≈ nm 2 ) and a high quantum yield (ν ≈ 0.8).

Quenching and photo-bleaching

At high excitation intensity, the fluorescence rate tends to reach saturation. Because of the high excitation intensity required to achieve fluorescence imaging, biological tissues can be damaged and unable to fluoresce, this process in called photo-bleaching.

Although photo-bleaching is usually not desirable in fluorescence imaging, new imaging techniques have taken advantage of this characteristic such as fluorescence recovery after photo-bleaching (FRAP), the related fluorescence loss in photo-bleaching (FLIP) or fluorescence localization after photo-bleaching (FLAP) [34].

Quenching is the loss of efficiency of the fluorescence process. It is often consequent to certain molecular interactions with the fluorochromes but many mechanisms of quenching exist. Some are presented Fig. 2.25

Biomedical optical imaging systems

A wide range of optical imaging system has been developed for biomedical tissues, with diverse contrasts based on differences in the optical properties defined previously.

When considering the highly scattering nature of biological tissues, these systems can be divided into categories depending on the scattering regime: the ballistic regime is used for thin transparent media, single scattering for surface imaging and multiple scattering for diffuse optics. Ballistic and single scattering-based optical imaging is the simplest form of imaging as it only requires geometrical approaches. It can consequently provide high resolution images but is yet limited to the penetration depth of the imaging wavelength in the tissue, usually to 100µm. On the other hand, diffuse imaging is based on the complex interactions of light with biological tissues. It is often based on a statistical description of multiple scattering to extract the scattering and absorption of tissues. The spatial resolution of these approaches is consequently limited and diffuse imaging is mainly performed at a macroscopic scale. Some common optical systems used for biomedical imaging are discussed here. Biomedical optical imaging has two principal interests: high-resolution and/or high-sensitivity images and new contrasts compared to the systems described previously. Biomedical optical imaging consequently rarely consists of a single camera: complex systems are elaborated for improving resolution and contrasts notably in microscopy. Among all the possible classification of the vast spectrum of developed approaches, it was chosen here to classify them according to the type of gating used to extract contrast or improve resolution. Indeed, the complexity of optical systems for biomedical imaging often stem from an encoding of the signal whether it is in time, space, spatial-frequency, coherence or polarization. Optical systems based on coherence-gating and interferometry are discussed in details in a later section.

Time-gated imaging

Time-gated imaging consists in extracting information from a time-resolved signal. Timegated imaging consists in sending an light pulse in a sample and recording the resulting signal resolved in time. By comparing the temporal shape of the impulse on the input and output, optical properties such as the scattering and absorption of a medium can be determined (see Fig. For a short pulse of light transmitted through tissue, the ballistic component emerges fist, followed by the diffuse component. The early par of the diffuse component is known as the "snake" component because the path of these photons approximate straight lines in the tissue. Image from from [30] 2.4.2.2 Space-gated imaging Space-gated imaging consists in confining the light used for imaging in a specific zone of the imaged sample. It is often combined with scanning to reconstruct an image. Spacegating is commonly found in confocal or light-sheet microscopy. This time of imaging solves a significant issue with imaging systems (such as in Fig. 2.24) in the imaging of thick tissues. Indeed, if infrareds take advantage of a greater penetration depth in tissues [36], fluorescence is emitted from different planes in thick tissues, including ones out of the focal plane of the objective. The quality of images is consequently deteriorated by this background signal. Space-gated systems suggest to overcome this limitation by focusing the excitation beam in the focal plane and scanning through the sample to reconstruct the image. This optical sectioning method efficiently improved imaged quality on thicker samples as highlighted in Fig. 2.28.

Two-photon-gated imaging

In 1990, two-photon microscopy was developed [37], allowing for additional increase of image quality and signal-to-noise compared to confocal systems. The principle is to use a high excitation intensity to generate the simultaneous absorption of two low energy photons in order to cause the emission of a photon with a shorter wavelength. Because this process is non-linear, it is only localized in the focal point of the laser illumination. The probability of out-of-focus emission is consequently negligible. Two-photon imaging can be considered as a form of space-gating. The principle of two-photon imaging can be generalized to the absorption of a greater amount of low-energy photons and is called multi-photon microscopy [38].

Fourier space-gated imaging

Instead of limiting the illumination in space, the spatial frequencies can also be gated. This can be performed experimentally by using spatial filters and lenses. A simple example of Fourier-space filtering is dark-field microscopy where the background of an image is extracted optically by a placing an opaque disk to mask the light that wasn't scattered by the sample. The resulting images typically show a better contrast than bright-field microscopy, especially on unstained samples. And example of dark-field image can be found in Fig. 2.30.

Spectrally-resolved imaging

All optical properties are wavelength-dependant. The spectrum considered for optical imaging is usually from the near-IR to UV-A (see Fig. 2.31). The variations in these optical properties through this spectrum are exploited in optical spectral imaging to retrieve contrasts relative to the tissue's chemical composition. The gating can occur on either the illumination or the detector: either a spectral detection system with a broadband light source or a broadband detector with successive illuminations at different wavelengths. Multi-spectral and hyper-spectral imaging have notably been developed in biomedical imaging for histopathology and surgery guidance as improvements on RGB cameras, allowing to image beyond the human visual sensitivity.

As described earlier, optical properties of tissues are highly dependent on the light frequency. Consequently, analyzing the spectrum of light after it interacted with biological tissues can provide useful information on the tissue composition and characteristics. For example, spectroscopy is used in medicine for label-free evaluation of some compounds' concentration such as drug or oxygen. In terms of imaging, the main methods presented here are Raman and Brillouin imaging.

Raman imaging.

Raman scattering occurs when an incident photon created a deformation of the electron cloud around an absorbing molecule, creating a polarity. Because the vibrational or rotational mode of the molecule is modified, when a photon is emitted by the excited molecule to return to its ground state, this photon can have a different frequency than the incident photon's. If the frequency of the emitted photon is lower than the incident photon's, the phenomenon is referred to as Stokes Raman scattering. And if the frequency of the emitted photon is higher than the incident photon's, the phenomenon is referred to as Anti-Stokes Raman scattering. Raman imaging consists in recording the frequency shift of emitted photons over a 2-D field of view in order to characterize the chemical composition and molecular structure of the tissues [41], [42]. The reconstructed images are usually colored using the known spectra of chemicals of interest.

Brillouin imaging.

Brillouin scattering is a phenomenon where light in scattered by acoustic vibrations. The direction and frequency of the acoustic wave determines the direction of the scattered photon and its frequency shift. The scattered photon can have either a lower energy (Brillouin Stokes) or a higher energy (Brillouin anti-Stokes). The frequencies for Brillouin and Raman scattering are presented in Fig. 2.32. Because Brillouin scattering is based on the interaction of acoustic waves with optical waves, Brillouin scattering allows for mechanical and viscous properties evaluation on a molecular scale. More details can be found in [40], [43]. § 2.5 Interferometry and coherent imaging techniques Interferometry was first described in the XVIIth century by I. Newton. Interferometry is the measurement of patterns created by the combination of two or more waveforms, called interferences. These pattern are created when the interfering beams have the same frequency but are out of phase. Their superposition consequently creates patterns with modified amplitudes. The analysis of these patterns can inform on the phase difference between the interfering wavefronts.

The difference of optical path between the interfering beams in regards to their wavelength is what determines the nature of the interaction. For the case of two interfering beams, the interference can be additive at a point in space if the optical path difference is a multiple of the light's wavelength: in this case, both beams arrive in phase and the resulting intensity is increased. On the other hand, if the beams arrive out of phase (the optical path difference is half the wavelength modulo the wavelength), the interference is said to be destructive, the resulting intensity is null.

Interferometric optical systems allow sub-micrometric distance differences to be visible to the naked eye. The visible spectrum being in the 400-800 nm range, exploiting light interference patterns leads to highly-sensitive systems, which are especially interesting for biomedical imaging.

The basic principles of interferometry are explained in this section, including the notions of speckle and coherence. The principle of Optical Coherence Tomography (OCT) is presented as it has been recently widely developed in biomedical imaging research. Since the focus of this work is on digital holography, the principle is detailed in the last part of this section. Some basic definitions are introduced here as they are instrumental to the rest of this Chapter.

Fourier transform.

The Fourier transform is a function commonly used in signal theory that decomposes a function in the frequency (or Fourier) domain. It is noted T F here. In one dimension, the Fourier transform, noted F of an integrable function f can be defined as:

F (ν) = 1 √ 2π f (x)e -i2πνx dx (2.33)
Coherence function. In signal theory, the coherence function between two signals x(t) and y(t) can be defined by Fourier transforms of correlation functions using Wiener-Khinchin theorem:

Γ x,y = F T [x ⊗ t y] F T [x ⊗ t x]F T [y ⊗ t y]
(2.34)

with F T the Fourier transform and t the correlation operator. F T [x ⊗ t x] and F T [y ⊗ t y] are the Fourier transform of the auto-correlations of the signals x(t) and y(t) respectively, which is also the spectral power density of these signals. The expression F T [x ⊗ t x] is the cross-spectral density.

Complex notation.

Similarly to the complex wave definition introduced in section 2.2, a complex notation is used here to describe the phenomenon of interferences. The real monochromatic electric field E is expressed here as [44]:

E(r, t) = Re[A(r)e -iωt ] = 1 2 [[A(r)e -iωt + [A * (r)e iωt ] (2.35)
with A the complex amplitude and ω the beam's angular frequency and Re designates the real part of a complex quantity.

The complex amplitude can be expressed in a Cartesian system:

A(r) =   A x (r)e ihx(r)
A y (r)e ihy(r) A z (r)e ihz(r)   (2.36) with A i and h i real functions. h i can be expressed as h i = k • r -ϕ i with k the propagation vector and ϕ i the phase constants informing on the state of polarization of light.

Alternatively, the complex amplitude can also be decomposed as:

A(r) = A 0 (r)e i(k•r+ϕ(r)) (2.37)
with A 0 (r) the real amplitude and ϕ(r) the phase.

From the previous expressions, plane waves and spherical waves can be defined from spatial symmetries.

Plane wave.

A wavefront is a surface where all points have the same phase. A plane monochromatic wave is characterized by plane wave fronts. It corresponds experimentally to a collimated beam. Using previous notation, its complex amplitude can be expressed by:

A(r, t) = Au x e i(k•r+ϕ) (2.38)
with A the amplitude u x the unit vector in the propagation direction.

Spherical wave. Spherical waves have spherical wavefronts. They correspond to the wave-field generated by an isotropic point source. Considering the center of the sphere as the origin or a spherical coordinates system, the expression of a spherical wave is:

A(r, t) = Au r r e i(k•r+ϕ) (2.39)
with u r the radial unit vector.

Intensity. Intensity is defined here as the time-averaged of the amount of energy crossing in a unit of time perpendicularly to the direction of energy flow. Considering plane electromagnetic waves, the expression of intensity can be written:

I = c 4π ϵ µ ⟨E 2 ⟩ (2.40)
with ϵ and µ the permittivity and permeability of the propagation medium.

In the rest of this chapter, the proportionality coefficient is neglected and it is considered that I = ⟨E 2 ⟩.

The optical waves have very high frequencies (375 to 750 THz in the visible spectrum), optical detectors usually recorded mean square value of the field over a large number of oscillatory periods. The recorded intensity only provides information on the amplitude of the electromagnetic field in the absence of interferences. Indeed, the intensity recorded can be written in function of the amplitude using the previous definitions [44]:

I = ⟨E 2 ⟩ = 1 2 A • A * = 1 2 (a 2 x + a 2 y + a 2 z ) (2.41)

Interferometry principle

The first case to be considered is two-waves interferometry. For two electromagnetic sources, noted E 1 and E 2 , the resulting intensity at a position r and time t can generally be written:

I(r, t) = ⟨(E 1 (r, t)+E 2 (r, t)) 2 ⟩ = ⟨E 1 (r, t) 2 ⟩+⟨E 2 (r, t) 2 ⟩+2⟨E 1 (r, t)•E 2 (r, t)⟩ (2.42)
The first two terms of this expression simply express the separate intensity contribution from each light source. 2⟨E 1 • E 2 ⟩ is the interference term, corresponding to the interactions of both electromagnetic fields.

Different cases are considered in this section based on this general expression. An analytical expression of the interference term is first provided for monochromatic light sources. The model for the source is then complexified by considerations of the source's spectrum and size. Notions of spatial and temporal coherence are introduced through these examples. A general definition of coherence is then provided.

Interferences of monochromatic plane waves at the same frequency

The simplest case in interferometry is the interaction of two perfectly coherent monochromatic electromagnetic waves. Plane waves are considered here. To simplify the notation, the light beams are considered linearly polarized waves with E in the direction of the unit vector x. The complex amplitudes A 1 and A 2 of the electric fields can be written:

A 1 (r) = A 1 xe i(k 1 •r+ϕ 1 ) A 2 (r) = A 2 xe i(k 2 •r+ϕ 2 )
(2.43)

with k 1 and k 2 the propagation vectors, A 1 and A 2 the amplitudes and ϕ 1 and ϕ 2 the phase terms.

The phase difference between the two beams can be introduced as a function of the optical path difference, noted δ:

(k 1 -k 2 ) • r + ϕ 1 -ϕ 2 = 2π λ δ + ϕ 1 -ϕ 2 (2.44)
with δ the optical path difference defined as the difference in distance traveled by the two light beams multiplied by the refractive index of the medium. The intensity from the superposition of the two fields can consequently be written:

I = 1 2 A 2 1 I 1 + 1 2 A 2 2 I 2 + A 1 A 2 2 √ I 1 I 2 ⟨cos(δ + ϕ 1 -ϕ 2 )⟩ (2.45)
In order to observe interferences, the phase terms ϕ 1 and ϕ 2 have to be constant in time relatively to the integration time. Indeed, rapid variations of the quantity ϕ 1ϕ 2 in time cause the average of the cosine to be null and eliminates the interference term. The phase difference ϕ 1ϕ 2 can often be considered constant or null. This assumption can be guaranteed experimentally by using the same source to create the interfering beams by using a division of wave-fronts (through slits) or amplitude (through semi-reflective media). The phase difference is consequently disregarded here to simplify notations.

The Eq. 2.45 can also be written using the coefficient γ to represent the interference pattern's contrast:

I = I 0 (1 + γcos(δ)) (2.46) with γ = 2 √ I 1 I 2 I 1 +I 2 and I 0 = (I 1 +I 2 ) 2 .
The interference pattern consequently depends on the phase difference between the two beams. If the propagation vectors are not collinear, the recorded intensity varies with the point of measure r. The surfaces with the same phase difference are the planes perpendicular to the direction of k 1 + k 2 as illustrated Fig. 2.33

Interferences of monochromatic plane waves at different frequencies

Let's consider now two optical plane waves with different frequencies ω 1 and ω 2 . The corresponding electric fields can be expressed similarly to the previous example:

E 1 (r) = A 1 xe i(k 1 •r-ω 1 t+ϕ 1 ) E 2 (r) = A 2 xe i(k 2 •r-ω 2 t+ϕ 2 )
(2.47)

with A 1 and A 2 the amplitude of the two fields and ϕ 1 and ϕ 2 their respective phases, considered constant.

In this case, the interference term is written:

1 (r, t) • E 2 (r, t)⟩ = 2A 1 A 2 ⟨cos(k 1 • r -ω 1 t + ϕ 1 -k 2 • r + ω 2 t -ϕ 2 )⟩ = 2A 1 A 2 ⟨cos((k 1 -k 2 ) • r -(ω 1 -ω 2 )t + ϕ 1 -ϕ 2 )⟩ (2.48)
In this interference term, the cosine varies in time with a frequency ω 1 -ω 2 . The interference term consequently depends on the integration time of the captor (T) compared to the frequency difference ∆ω = ω 1 -ω 2 . If, ∆ω >> 1/T , the second cosine is null and no interference can be observed. On the other hand, if the sources' frequency difference is small enough in comparison to the integration frequency, temporal variations in the recorded intensity can be observed. This is referred to as heterodyne interferometry. Heterodyne interferometry can be used in biomedical imaging to perform optical Doppler imaging or increase the SNR of an interferometric system.

A monochromatic description of light sources is used in these two first examples. However, the monochromatic representation is not sufficient to describe experimental sources.

Young's slit with a polychromatic point source and temporal coherence

Most natural sources cannot be considered as monochromatic. The source's model is refined here by considering a polychromatic source. This source is considered here as the superposition of monochromatic waves. The resulting intensity is consequently the sum of the monochromatic interference patterns.

In order to address the interference pattern created by two light beams from the same polychromatic light source, an example of interferometry by division of wavefront is chosen here. The principle is however the same than if division of amplitude was used to create the two interfering wavefronts.

In the case of a division of wavefront, the simplest configuration is Young's slit experiment. The case of a punctual isotropic polychromatic source in front of two pinholes is considered here with the main parameters defined in in Fig. 2.34. The axis of the pinholes is defined as the x-axis to simplify notations. In this case, the pinholes act as two secondary point sources of light. The electric fields E 1 and E 2 are defined at the pinholes O 1 and O 2 for a given frequency ω as:

E 1,ω (O 1 , t) = s 0,ω P O 1 e i(ωt+ϕω(t)) E 2,ω (O 2 , t) = s 0,ω P O 2 e i(ω(t-τ 1 )+ϕω(t-τ 1 )) (2.49)
with τ 1 = (P O 2 -P O 1 )/c the delay between the light arriving at O 1 and O 2 and s 0,ω and ϕ ω the amplitude and phase of the light at S considering a frequency ω.

The waves field from these secondary sources can be propagated to a point M of the screen:

E 1,ω (M, t) = s 1,ω (M )e i(ωt+ϕω(t)) E 2,ω (M, t) = s 2,ω (M )e i(ω(t-τ )+ϕω(t-τ ))
(2.50) The resulting intensity in M is consequently expressed for one wavelength as:

with τ = τ 1 + (O 2 M -O 1 M )/c
I ω (M ) = s 1,ω (M ) 2 + s 2,ω (M ) 2 + 2s 1,ω (M )s 2,ω (M )⟨cos(ωτ -ϕ ω (t -τ ) + ϕ ω (t))⟩ (2.51)
We consider here that time delay τ of between the two wavefronts interfering is small enough to consider the phase difference ϕ ω (t -τ ) -ϕ ω (t) negligible. The expression retrieved here can be written similarly as Eq. 2.46:

I ω (M ) = I 0 (M, ω)(1 + γ(M, ω)cos(ωτ ))
(2.52)

with I 0 = s 2 1,ω + s 2 2,ω the separate intensity contributions of the secondary sources in O 1 and O 2 and γ the contrast.

If this intensity is now integrated over the source spectrum, the total intensity can be expressed:

I(M ) = ω I 0 (M, ω)(1 + γ(M, ω)cos(ωτ ))dω (2.53)
In the case of a source with a spectrum centered on ω 0 and uniform spectral density over a width ∆ω , Eq. 2.53 can be written: In the case described here, the sinc term modulates the contrast of the interference fringes. From the equation, it can be observed that the contrast of interferences pattern in the case of a polychromatic source is a function of the time delay between τ between the interfering wavefronts. In this example of a source with uniform spectral density, interferences can only be observed with | ∆ω 2τ | < π. The corresponding time τ c can be defined as the coherence time and can be expressed in this case:

I(M ) = I 0 (M )∆ω[1 + γ(M )sinc( ∆ω 2τ )cos(ω 0 τ )] (2.
τ c = 2π ∆ω (2.55)
This coherence time corresponds to the maximum time delay between two points along the propagation axis capable of interfering.

Young's slit with an extended light source and spatial coherence

The previous examples have explored cases where the light sources were considered punctual. In order to better represent experimental sources, the case of a source with a finite size is considered here. We consider here the extended light source as a superposition of incoherent monochromatic point sources over its surface. The interference pattern consists in the sum of the interference patterns from each point on the source's surface, which decreases the resulting interference pattern's contrast. If we go back to the pinholes configuration, we can now consider a monochromatic source with a width h in the X-direction, centered on the propagation axis with uniform intensity I 0 . In order to simplify notations, only the spatial extent in the direction of the slits is considered here. The resulting pattern results once again on an integration of the interference pattern over the width of the source:

I(M ) = h/2 -h/2 I 0 (1 + γcos(2π ax λD + 2π aX λL ))dX (2.56)
with X the position of the source and x the position of M. The other quantities are defined in Fig. 2.34. Integration of Eq. 2.56 gives:

I(M ) = I 0 h(1 + γsinc( πah λL )cos( 2πax λD )) (2.57)
Similarly to the case of a polychromatic source explored previously, the interference pattern can be observed for ah λL < 1. In this case, a spatial coherence length can be defined as:

Λ c = πL/h (2.58)
This distance corresponds to the maximum distance between the pinholes to observe interferences. More generally, the spatial coherence length is defined as the maximum distance between two points of a wavefront capable of interfering.

Generalized notion of coherence

Through the previous examples, the notions of spatial and temporal coherence have been introduced and considered separately in the example of Young's slit experiment. The cases described are yet very simple (uniform spectral density or uniform intensity distribution). In this paragraph, more general definitions of coherence are provided.

As introduced previously, spatial coherence describes the ability for two points in space to interfere and is linked to the intensity distribution of the light source. On the other hand, temporal coherence is a function of the source's spectrum. Temporal and spatial coherence are however different interpretation of the same general notion.

In the case of an arbitrary field of complex amplitude A(r, t), the complex degree of coherence γ c is a complex function defined for two points O 1 and O 2 for a time delay τ between the field in O 2 compared to O 1 :

γ c (O 1 , O 2 , τ ) = ⟨A(O 1 , t)A * (O 2 , t -τ )⟩ √ I 1 I 2 (2.59)
with I 1 and I 2 the intensity at O 1 and O 2 .

This function describes both the temporal and spatial coherence, respectively

γ c (O 1 , O 1 , τ ) and γ c (O 1 , O 2 , 0).
General expressions for both the temporal coherence and the spatial coherence are presented here.

Temporal coherence.

The temporal coherence function can be expressed using Wiener-Khintchine theorem:

γ c (O 1 , O 1 , τ ) = T F [S(O 1 , ν)] (2.60) with S(O 1 , ν) = |T F [A(r, t)]| 2 the normalized spectral density of the source in O 1 .
Temporal coherence is generally determined by the source's central wavelength λ 0 and its spectral width ∆λ. The light sources' spectrum are commonly considered as a Gaussian distribution. Considering a stationary field, the coherence length l c defined as the maximal distance between two wave-fronts capable of interfering can be written:

l c = λ 2 0 n∆λ (2.61)
with n the refractive index of the medium.

Spatial coherence.

Let's consider an extended light source of surface ω contained in a plane. Each infinitesimal surface element dω at a central point S has an intensity dI = I(S)dω. Van Cittert-Zernike theorem expresses γ c for two points O 1 and O 2 of a quasi-monochromatic (δω << ω 0 ) extended light source:

γ c (O 1 , O 2 ) = 1
I(O 1 )I(O 2 ) σ I(S) e ik(R 1 -R2) R 1 R 2 dS (2.

Speckle

Interferometry allows very high sensitivity at a scale close to the light's wavelength. However, at this scale, most surfaces cannot be considered as flat [45] and in scattering media, differences of refractive index can also vary at a micrometer scale. The interaction of light with a surface or a medium that is randomly modulated at a scale close to its wavelength creates a granular interference pattern. This phenomenon and pattern is known as speckle. The appearance of speckle differs from the fringes and rings usually created by interferometers because of its random structure, illustrated in Fig. 2.37. Speckle can be observed in transmission or in reflection and is present in most coherent imaging techniques, whether it is utilized to reconstruct signal or filtered out as noise.

As described in section 2.4.1, biological tissues are very scattering in the visible spectrum of light. When illuminating a biological tissue with coherent light, the phase of the reflected wavefronts cannot be considered as deterministic as it is greatly modulated by the presence of scatterers, losing its spatial coherence properties from multiple scattering occurring within the sample [46]. In order to perform interferometric imaging on biological tissue, a basic description of the speckle phenomenon is proposed here to characterize its properties and understand its contributions in amplitude and phase imaging. A statistic approach is necessary to describe this phenomenon and detailed descriptions can be found in literature [47] [48] [45]. A basic description is provided here for monochromatic light in the context on interferometry.

Speckle can be modeled by the sum of waves with a random distribution in phase and amplitude. The complex amplitude can be expressed as the superposition of secondary sources from a monochromatic field A n . The total complex amplitude A at point M can be expressed:

A(M, t) = 1 √ N N n=1 A n (r, t)e iϕn(r,t) (2.63)
with N the number of scatterers.

First order approximation

A uniform distribution of scatterers is assumed and A n and ϕ n are considered independent variables with ϕ n uniformly distributed in ] -π, π].

According to the central limit theorem, when considering a large number of scatterers N, this expression tends to a complex Gaussian process. The real and imaginary part of the function A(r, t) are consequently statistically independent at any point and distributed with a mean value of zero. From this, it can be concluded that the modulus of the complex amplitude has a negative exponential distribution with a variance of ⟨I⟩/2 and the phase a random distribution in the interval ] -π, π] [47], [48]. Their probability densities p I and p ϕ are:

p I = 1 ⟨I⟩ e -I ⟨I⟩ p ϕ = 1/2π (2.64)
This first expression of the statistical profile of speckle corresponds to the histograms obtained experimentally (Fig. 2.38). This description is however insufficient to properly describe the granularity of the speckle pattern. The case of speckle after propagation through an optical system is considered here.

Spatial structure through an optical system

A lens is considered here between the scattering surface and the image plane, as illustrated Fig. 2.38. Considering D >> y, the phase difference between the wave arriving from P 0 at the point P (x, y, z) in regards to the wave arriving at P ′ 0 (0, 0, z) can be expressed as:

δ = ndsin(θ) ≈ ndy/D (2.65)
It is assumed here that the medium of propagation is air. If we consider P ′ 0 the center of the speckle pattern and P its border, it can be written that y = d speckle /2 with d speckle the size of the pattern and δ ≈ λ/2. Eq. 2.65 can be written:

d speckle = 2δD/d ≈ λ 2 N A (2.66)
It results that the average lateral (regarding the propagation axis) speckle grain size d speckle is proportional to the wavelength and inversely proportional to the numerical aperture N A = nsin(α).

Let's now consider the axial speckle grain size. The point P ′′ 0 (0, 0, z+dz) is considered in regards to the point P ′ 0 (0, 0, z). The path difference between both points is written:

δ = ndz α 2 2 (2.67)
Similarly as for the lateral extension, we consider here P ′ 0 the center of the speckle grain and P ′′ 0 its border. We have consequently dz = l speckle /2 and δ ≈ λ/2. The axial extension of the speckle pattern l speckle can be expressed as:

l speckle = 2 λ N A 2 (2.68)
In air, the numerical aperture is a value between 0 and 1, the speckle pattern can consequently be considered in 3D as a prolate ellipsoid shaped with the major axis in the direction of light propagation. As for the temporal variations of the speckle pattern, it depends essentially on the source's coherence properties, as described previously.

Speckle in interferometric imaging

Speckle patterns are almost inevitable when performing coherent imaging in scattering media. Because the properties of speckle vary with illumination and the sample's surface optical properties and movements, some imaging techniques take advantage of these characteristics and measure speckle statistics to extract source or surface properties or track movements within the sample (speckle tracking).

In imaging techniques such as holography or OCT, the speckle grain is considered to be encoding the information of interest. Speckle consequently has to be resolved in these imaging techniques. This is usually achieved experimentally by tuning the numerical aperture of the imaging system based on Eq. 2.66 and 2.68. Speckle can yet be a significant issue and hinder the resolution of small structures. Because the contrast of interest in these techniques emerges from the same phenomenon as speckle, the suppression of the speckle pattern is a difficult problem. Techniques found in literature often require to reduce light coherence, to acquire multiple images and/or to use numerical processing [49]. This denoising comes at the cost of either longer acquisitions times or/and more complex set-ups or loss in resolution.

Michelson and Mach-Zehnder interferometers

Optical interferometry systems emerged at the very end of the XIXth century and carry the name of their inventors. To this day, most of the interferometers found in studies are derived from these basic optical designs. The Michelson and Mach-Zehnder interferometers are described here as well as their use and characteristics as they are the most common type of two-waves interferometers in literature.

Michelson interferometer.

Michelson interferometer is described Fig. 2.39. Three configurations can be distinguished:

• In the first configuration, both mirrors M1 and M2 are orthogonal to the direction of the incident light beam but are at a different distance from the beam-splitter.

The projection of the M2 mirror on the vertical axis highlights this offset. The optical path difference δ can be written as a function of the distance difference d between M1 and M'2 :

δ ≈ 2nd(1 - θ 2 2 ) (2.69)
with n the refractive index of the medium and θ the angle between the considered position in the image plane and the mirrors' center.

In this configuration, the optical path difference is only a function of the distance from the considered point to the center of the image. The obtained interferometric figure is expected to have rotational symmetry. If the distance d is large enough, fringes are expected to be visible on the interference figure.

• In the second configuration, there is no distance offset between both mirrors but one is tilted of and angle θ. The optical path difference is now:

δ ≈ 2nθx (2.70)
with x the lateral position of the considered point in the image plane. The x-axis is perpendicular to the rotation axis of M2. The consequent interference figure is a set of evenly distributed lines.

• At last, the third configuration is for sample analysis. In this case, both mirrors have no tilt nor offset: the only source of path difference is the sample placed in one of the interferometer's arms. The interference pattern recorded on a camera can be analyzed to extract the relative phase and amplitude images of the sample. This configuration is often used in OCT imaging systems, as described in a following section. Michelson and Mach-Zehnder interferometer are both amplitude division interferometers. This optimizes the power in both arms of the interferometers compared to the use of pinholes described earlier where most of the source power did not contribute to the retrieved intensity in the image plane. The Mach-Zehnder and the Michelson interferometers are really close in principle with equivalent compactness and stability. The main difference is that the Michelson interferometer requires only a single beam-splitter where the Mach-Zehnder needs two: this reduces costs but also reduces the degrees of freedom for alignment. In the case of OCT, a Michelson configuration is generally chosen while a Mach-Zehnder configuration is preferred in this work for off-axis digital holography since one of the beam is tilted and additional degrees of freedom facilitate alignment (see Chapter 4). More details on both OCT and digital holography can be found in the following paragraphs.

OCT

Optical coherence tomography (OCT) has stormed into the scientific research in the 1990 and its rapid advances have allowed for a wide range of application, notably in biologicaltissue imaging. OCT can be defined as optical ultrasonography. Instead of recording echoes from acoustic impedance mismatches, OCT images reflected optical waves, probing the medium's back-scattering properties in-depth. Because light travels in biological tissues around 2 × 10 8 m/s, no detector is able to directly image the optical echoes. OCT differs from ultrasound techniques by the use of interferometry to achieve imaging of the optical signal. Indeed, sources with short coherence time are used to record the echoes, encoding the time-delay in the interference pattern. Contrary to holography, OCT is based on low-coherence sources and can be referred to as low-coherence interferometry or white-light interferometry.

OCT achieves imaging depth in the order of magnitude of the millimeter when using red or near-IR, but allows for micrometer resolution imaging. This technique is especially promising for imaging through transparent media like the eye. Indeed, the first OCT images were performed in 1991 [50] and five years later, the first commercial OCT system for ophthalmic diagnostics was introduced, leading to a great number of clinical studies [51]. OCT has also been adapted in a variety of medical fields: dermatology [52], cardiology [53], gynecology [54], odontology [55]... A brief description of the OCT principle and of the different configurations is provided here. This section aims as providing the basis to understand a number of publications of interest later referenced throughout this thesis.

Principle

OCT is based on a low coherence source of illumination and an interferometry set-up. The principle is presented Fig. 2.41. The OCT system described in this figure only allow for axial scanning: a galvanometric mirror has to be added to the sample illumination beam in order to reconstruct 3D-images. Both beams recombine on the camera to form an interference pattern. This interference pattern varies depending on the optical path difference between both beams. For low coherence light, the interference signal is maximal when this optical path difference is null and is indistinguishable when the path difference is greater than the coherence length of the source.

In this simple configuration of OCT, we consider here the following signals for the object and reference waves [45]:

E r (t) = α r s(t -2l r /c) E o (x, z, t) = α o (x, z)s(t -2l o /c -2nz/c) (2.71)
with s the signal from the low-coherence light source, α r and α o the attenuation of each beam through the beamsplitter, n is the average refractive index and l o and l r are respectively the distance travelled by the beam between the object and the beamsplitter and between the reference mirror and the beamsplitter.

In order to simplify notations, we consider here α o (x, z) = α r = 1/2. The recorded intensity can consequently be expressed as:

I(x, z) = I r + I o (x, z) + 1 4 [⟨s(t)s * (t + τ )⟩ + ⟨s * (t)s(t + τ )⟩] (2.72)
with τ = 2(l r -l o )/c -2nz/c and I r and I o the intensities of the reference and object beam respectively.

By defining the normalized source function N s (τ )e -iωτ = ⟨s(t)s * (t+τ )⟩ ⟨|s(t)| 2 ⟩ , Eq. 2.72 can be written [45]:

I(x, z) = I r + I o (x, z) + Is 4 |N s (τ )|cos( 4π λ 0 (l r -l o -nz) -ϕ s (τ )) (2.73)
with I s = ⟨|s(t)| 2 ⟩, λ 0 and ω the source's mean wavelength and center angular frequency and N s (τ ) = |N s (τ )|e iϕs(τ ) .

The normalized source function N s (τ ) can be expressed as:

N s (τ ) = T F [S(ν)] (2.74)
with S(ν) the normalized power spectral density of the source. Let's consider a mirror scanning in-depth at a speed V r and a scatterer at a position (0, z k ). Considering x=0, the corresponding intensity is expressed:

I(0, z k , t) = I s 2 [1 + |N s (τ ′ )|cos( 4π λ 0 (V r t -nz k ) -ϕ s (τ ′ ))] (2.75) with τ ′ = 2 Vrt c -2n z k c .
The interference term consists of two components: one depending on the source's spectrum (N s ) and one depending on the scatterer (cos). This results in a modulated signal in time as illustrated Fig. 2.42. The shape of the envelope of the signal depends on the source while the rapid oscillations in signal are due to mirror motion. The maximum of the signal is reached when the optical path difference between the mirror and the scatterer is null. Considering a source with a Gaussian power spectrum of bandwidth ∆ω, N s (τ ) can be expressed as [45]:

N s (τ ) = e -( ∆ω τ 4 √ ln(2) ) 2 (2.76)
OCT axial resolution δz can be determined from combining Eq. 2.75 and Eq. 2.76

δz = 4c ln(2) n ∆ω (2.77)
By expressing ∆ω ≈= 2π c∆λ nλ 2 and adding the expression for the lateral resolution δx and depth of field (DOF) [56], the following general characteristics of OCT are obtained:

     δz = 2 ln 2λ nπ∆λ δx = 4λf πd DOF = nπ(δx) 2 λ (2.78)
with ∆λ the bandwidth of the light source, f the focal length of the objective and d the spot size of the incident beam on it.

As seen in the previous equations, one of the main advantage of OCT is that the axial resolution is only dependent on the source spectrum while the lateral resolution depends on optical system. This decoupling between the axial and lateral resolution is one of the main advantage in OCT. This has allowed cross-sectional images with a 10 to 100 times better resolution compared to ultrasound imaging [51].

Types of OCT

OCT techniques can be divided into time-domain OCT and frequency-domain OCT [57]. A brief description of those term is introduced here.

Time-domain OCT.

Time-domain OCT utilized broadband light sources. As described in the previous example, the image can be reconstructed by scanning the reference mirror along the optical axis to acquire the signal using a single detector and recording the envelope of the interference signal. The amplitude of the successive interference signals corresponds to each scattering layer detected in the sample arm.

Fourier-domain OCT.

Fourier-domain OCT consists in encoding the depth information in the frequency, allowing for faster imaging. Fourier-domain OCT can be divided into spectral-domain OCT and swept-source OCT. Spectral domain OCT consists in using a continuous light source and a spectrometer to disperse the interference signal into separate frequency components. As for swept-source OCT, the source's emission wavelength is swept to separate frequency components.

Full-field OCT. Contrary to conventional OCT that requires a lateral scanning to achieve tomographic imaging, full-field OCT provides two-dimensional images orthogonal to the optical axis without any scanning of the illumination light beam. The entire field of the image is illuminated with low spatial coherence light and a camera records interferometric images [58]- [60].

Digital holography

Holography was first introduced as wavefront reconstruction by its inventor D. Gabor in 1948 [61]. The principle is that from the recording of an unfocused diffraction pattern from an object, the image of the object can be reconstructed. The recording of this diffraction pattern is called an hologram. Holography can be applied to the full electromagnetic range but only optical holography is considered here.

Optical holography is based on the recording of the interference pattern created by the superposition of the diffuse light from a sample and a reference beam. The resulting interference pattern is then recorded, containing amplitude as well as phase information.

In order to record the hologram, a photosensitive film was first used. It's only more than 20 years later that digital holography theory and first set-ups were built, trying to record holograms and reconstruct images [62]- [64]. Digitization of holography began in the 1970s with the works of Konrod [65].

Since the beginning of this century and the development or high resolution imaging sensors, optical digital holography has seen a great range of application, notably in polarimetric and quantitative phase-contrast microscopy and in biomedical applications [66]- [68]. These holographic set-ups can be divided into two categories : in-line and off-axis.

The basic principle of holography is described in this section, including the differences between in-line and off-axis configurations as well as the methods for digital holography image reconstruction. The first configuration explored is in-line holography. Fig. 2.43 presents a simple in-line configuration. To ensure coherence, the light beam emerging from the source is separated into two beams: the reference beam, that directly illuminates the photosensitive film or camera and the object beam, which illuminates the sample that then diffuses light. Assuming perfect coherence between the two beams, an interference pattern is formed when they overlap on the image captor, whether is it a camera or a photosensitive medium. The case of a photosensitive film is first described here.

In-line holography

Once the hologram is recorded, if the photographic film is then re-illuminated with the same reference beam, the wavefront is reshaped as if it was originally coming from the object. Consequently, by observing the wavefront, the presence of the object is perceived. On the other hand, if the photosensitive film is illuminated from the back with the conjugate of the reference beam (the wave fronts and the propagation direction is opposite to the reference beam during the recording), the reshaped waterfronts are now converging to the real object [69].

To understand the physics behind holography, the light intensity can be described using the tools presented earlier for interferometry. The intensity (I) at the hologram plane can be described in function of the light field arriving from the reference beam (E 1 = A 1 e iϕ 1 ) and the object beam (E 2 = A 2 e iϕ 2 ):

I ∝ (E 1 + E 2 )(E 1 + E 2 ) * = ||E 1 || 2 A 2 1 + ||E 2 || 2 A 2 2 + E 1 E * 2 + E * 1 E 2 2A 1 A 2 cos(ϕ 1 -ϕ 2 )
(2.79)

with * designating the complex conjugate.

We consider here the photosensitive plate to a have a linear intensity response.

If we consider now the first case, where the hologram is being illuminated with the same reference beam, the transmitted intensity I t can be expressed:

I t ∝ E 1 (||E 1 || 2 + ||E 2 || 2 ) + E 1 E 1 E * 2 + ||E 1 || 2 E 2 (2.80)
The last term in this equation is indeed proportional to the diverging wavefront initially coming from the object.

On the second case, the conjugate of the reference beam is used to illuminated the hologram. The resulting transmitted intensity

I t ∝ E * 1 (||E 1 || 2 + ||E 2 || 2 ) + ||E 1 || 2 E * 2 + E * 1 E * 1 E 2 (2.81)
The second term of the equation here is proportional to the conjugate of the diverging wavefront initially coming from the object. This correspond to a converging wave creating a real image of the recorded object.

First optical holography experiments were carried out in Gabor's original in-line configuration [69]. This configuration is also designated as in-line holography or on-axis holography. In Gabor holography, the observer is in a configuration where it sees the superposition of the illumination source, a virtual image of the object and a real image of the object. This overlap is the principal drawback of in-line holography since it hinders contrast in image reconstruction.

Off-axis holography.

In order to solve the issue with image superposition of in-line holography, Leith and Utpatnieks presented the first off-axis holography set-up in 1962 [63]. The principle lies in the introduction of an angle between the reference beam and the object beam to separate the real image from the twin-image and zero-order diffraction in the spatial frequency domain. This method was first limited by the short coherence length of the available light sources. Holography was indeed mainly performed using Mercury lamps (l c ≈ mm). The progressive introduction of lasers greatly contributed to the development of interferometric methods, including holography. The angle introduced in off-axis holography between the beams allow the image, the twin image and the zero order of diffraction to propagate in different directions. As illustrated Fig. 2.46, this is equivalent to a separation of the components of the hologram in the spatial domain. In the case of off-axis holography, the reference beam can be described as

E 1 = A 1 e -i 2π
λ xsinθ . The intensity recorded in the hologram plane can consequently be expressed as [70]:

I(x, y) = A 2 1 + A 2 2 (x, y) 0 order + A 1 e i 2π λ xsinθ E 2 real image(+1 order) + A 1 e -i 2π λ xsinθ E * 2 virtual image(-1 order) (2.82)
The choice of the angle between the beams has to be made carefully in order to avoid any overlapping between the spatial components. Indeed, the angle (θ) between the beams modulates the carrier frequency of the interferogram, it must not exceed the Nyquist frequency of the detector [70]:

θ ⩽ arcsin λ 2∆x (2.83)
with ∆x the pixel's size and λ the light source's wavelength.

In order to optimize the use of the space bandwidth of the captor, θ must also verify [71]:

θ ≥ arcsin( 3 N A 2 √ 2πM ) (2.84)
with NA the numerical aperture of the system and M its magnification.

Digital image reconstruction

Instead of using photographic plates, digital holography is based on the use of a camera and numerical methods for image reconstruction. The principle is the same as described in the previous examples for recording the hologram, a camera just replaces the photosensitive film. As for the image reconstruction in digital holography, the process is entirely performed numerically. Instead of experimentally re-propagating the light to image the object, Fresnel's diffraction or near-field formalism can be used to simulate light propagation and compute the object image from the real image part of the recorded hologram. Fresnel diffraction equation can be interpreted as the sum of elementary wavelets coming from all points of the z=0 plane. The Fresnel integral describes the electric field diffraction in free space. Using notations presented Fig. 2.47, the field in the reconstructed plane (x, y, 0) can be written from the field in the hologram plane (ξ, η, z) as [72]:

E(x, y, z) = -i z λ ∞ -∞ E(ξ, η, 0) e i2πr λ r dξ dη (2.85) Under Fresnel's approximation (z 3 >> 1 8λ ((ξ -x) 2 + (η -y) 2 ) 2
), the vector r can be developed to obtain the following expression of the integral and Fresnel transform F is defined:

E(x, y, z) = -i 1 λz exp(iπ x 2 + y 2 λz ) ∞ -∞ E(ξ, η, 0)exp(iπ ξ 2 + η 2 λz )exp(-2iπ ξx + ηy λz ) dξ dη (2.86)
E(x, y, z) = -i 1 λz exp(iπ x 2 + y 2 λz )F x λz , y λz [E(ξ, η, 0)e iπ ξ 2 +η 2 λz ] (2.87)
In the case of off-axis digital holography, the image reconstruction consists first in isolating the real image from the rest of the signal on the interferogram and then numerically re-propagating it to extract the object's image. Different reconstruction methods can be used depending on the recording conditions. Fourier-based approaches (using only one fast Fourier transform) are well-adapted to far extended objects while convolution approaches are better for small objects close to the image sensor [73]. Digital holography is discussed in Chapter 4 when describing the experimental set-up. § 2.6 Multi-wave imaging Optical imaging in medicine present some limitations due to the nature of light propagation biological tissues such as low penetration for depth imaging and degraded resolution due to scattering. To overcome these limitations, new approaches have been developed at the interface with ultrasound imaging. Different imaging methods have indeed been combined to achieve enhanced performances either by using light and sound interactions or by translating acoustic medical imaging methods to optics. The resolution in optical deep tissue imaging was improved with photo-acoustic or acousto-optic approaches and new contrasts were imported from ultrasounds to optics such as Doppler angiography or elastography.

Acousto-Optic Imaging

Contrary to visible light, acoustic waves undergo few scattering in biological tissues and allow for depth imaging with millimeter spatial resolution. Combining acoustic and optical waves could allow imaging systems to benefit from the advantages of ultrasound imaging but with optical contrasts. Acousto-optic methods emerged in the 1980s [74]. They combine a highly coherent illumination on tissues with an ultrasonic beam. The ultrasonic beam introduced local refractive index changes in the tissues, which is visible on the speckle figure. Because the acoustic pulse can be focused in a small volume ( ≈ 1mm 3 ), only the photons going through this volume are scattered by vibrating light-scattering structures. These photons are said to be tagged by the acoustic wave. Tagged photons are characterized by side-bands on their spectrum corresponding to the optical frequency plus or minus the ultrasonic frequency (see Fig. 2.48). Synchronized detection then allows to recover only tagged photons, corresponding to a very localized region in the tissues [75]. The recover light intensity is quite small in comparison to the incoming laser intensity and a scanning is required to obtain depth imaging. But compared to purely optical approaches like diffuse optical tomography, acousto-optic methods offer higher resolution and using interferometric approaches provide high-sensitivity [76]. The photo-acoustic effect was first observed by A. G. Bell in 1881: it was noticed that after sending a focalized and modulated beam in gas, an audible signal could be heard at the modulation frequency. The principle of the photo-acoustic effect is that the illuminating electromagnetic wave is partly absorbed by the medium, producing heat and an expansion of the medium, creating an acoustic wave.

Photo-acoustic Imaging

In biomedical imaging, a light source is used on an absorbing tissue.

The created acoustic signal is then recorded using a matrix of transducers. An absorption map of the tissue can then be created. Non-invasive imaging can be performed using the photoelectric effect as well as functional imaging [77].

Doppler Optical Imaging

As described before for ultrasound imaging, Doppler blood flow imaging is mainly based on the reflection of waves on moving red blood cells in the blood stream. The reflected wave undergoes a slight frequency shift which is proportional to the speed of the moving cell. Consequently, the frequency shift provides quantitative information on the velocity of red blood-cells and the intensity of the recorded Doppler-shifted signal provide the number of red blood cells. If the imaging system is able to image the blood vessel's cross section, the blood flow can be determined quantitatively. Optical Doppler images have mainly been performed using interferometry methods such as holography [78] or OCT [79]- [81]. In these case, inferences between light which frequency has been shifted due to the Doppler effect and a reference beam is observed. By exploiting the linear relation between the intensity of the interference pattern and the cosine of the frequency shift presented in Eq. 2.48, the velocity of red blood cells can be mapped. OCT-based Doppler and OCT angiography both provide non-invasive three-dimensional blood flow visualization with high resolution. OCT angiography is based on motioncontrast imaging achieved from phase-sensitive set-ups with high temporal resolution. The decorrelation signal (differences in the back-scattered OCT signal intensity or amplitude) between sequential OCT B-scans provide this motion-contrast. The application of this technique to retinal or cerebral blood flow imaging seems especially promising since standard techniques are based on the use of dyes such as fluorescein and indocyaninegreen [79], [80], [82], [83].

Elastography Imaging

Elastography is the measurements of mechanical properties of a medium. Elastography contrast has been used for centuries in medicine through palpation. Increased stiffness in biological tissues is linked to structural anomalies such as fibrosis, calcifications or lesions. Because palpation only screens for lesions with an already considerable size and cannot provide accurate monitoring or quantitative assessment, providing physicians with elastography-contrast imaging systems could be of great help in diagnostics. Evaluating the stiffness of biological tissues was first achieved quantitatively through mechanical testing. However, this invasive method was not well-suited to medical applications. The term elastography was introduced at the beginning of the 90s [84] with the first non-invasive methods of stiffness assessments for biological tissues. Elastography imaging methods were developed first in ultrasound imaging but MRI quickly followed. It's only in 1998, with the introduction of elastography in OCT By Schmitt [START_REF] Sandrin | Transient elastography: A new noninvasive method for assessment of hepatic fibrosis[END_REF] that the field of optical elastography emerged. Elastography methods are defined and detailed in the following Chapter. § 2.

Conclusion

In this Chapter, the common modalities for medical imaging have been discussed with their principal clinical applications and limitations. As this thesis concentrates on in-terferometric optical imaging, optical imaging of biological tissues is especially detailed, including coherent imaging and the principle of digital holography, which was implemented experimentally. The experimental set-up is composed of two parts: an imaging system and an elastography approach. The basic principle and context surrounding the chosen imaging method for the experimental set-up have been presented in this Chapter. In order to have a more comprehensive understanding of the set-up, mechanical properties of tissues and the existing elastography approaches are presented in the following Chapter.

[85] J. M. Schmitt, "OCT elastography: Imaging microscopic deformation and strain of tissue," Optics Express, vol. -Chapter 3 -

Mechanical properties of biological tissues and their measurements § 3.1 Summary

Mechanical properties of materials describe their behavior under all types of forces and consequently their stability in specific environments. Mechanical properties are closely related to a material's composition and structure at different scales. Assessing the mechanical properties of a sample can consequently provide precious information. This information can be used in a variety of applications such as the assessment of new materials for specific applications, the mapping the Earth's structure or biomedical imaging.

In order to set the motivation for this work, physical quantities are defined in order to qualify and quantify a material's mechanical properties. From these quantities' definitions, models have been established for biological tissue's. These models describe the behavior of tissues when subjected to external forces. The second section of this Chapter defines the different mechanical properties used to describe biological tissues and general notions necessary to understand elastography, including wave propagation in tissues and the models used for biological tissues.

Based on the previously introduced notions, different approaches have been developed for determining the elastic properties of materials. These elastography approaches all rely on the same principle. At first, a mechanical constraint is applied on the sample of interest. Secondly the response of the sample to this stress is recorded. Then, using the models defined previously, mechanical properties such as stiffness, anisotropy or viscosity can be extracted. These different types of mechanical testing and elastography imaging are introduced in a third section. Their respective principle, advantages and limitations are discussed.

For biological tissues, an abnormal stiffness has been notably linked to a number of pathologies [1]- [3]. At a cellular scale, different cell properties and mechanism can be highlighted [4]. Retrieving stiffness information is a promising tool for diagnostics and the set-up developed in this thesis aims at developing such technologies. Potential applications of elastography in medicine and motivations for this work are presented in the last section. § 3.2 Theory, definitions and models of mechanical behavior

The notion of stiffness discussed previously characterizes the material's ability to resist deformation. The larger the stiffness, the greater the force required to cause a given deformation. The notion of stiffness is yet qualitative and does not correspond to a defined physical quantity. In this section, the physical quantities corresponding to this notion of stiffness are presented as well as some models that are the cornerstone of most elastography approaches. The application of these definitions and models to biological tissue is also specifically examined.

General definitions of mechanical properties

Biological tissues present a high diversity in terms of mechanical behavior. Their complex structures and various molecular compositions lead to a diverse range of materials: bones, muscles, fat... These different biological tissues all show very different mechanical behaviors characteristic of their function. Definitions are provided in this section to characterize these different mechanical properties through some physical quantities.

Mechanical properties are defined by a material's response to a constraint. It is worth noting here that a wide range of mechanical constraints can be applied on a material, as illustrated in Fig. 3.1 and the defined mechanical quantities depends on the type of loading the material undergoes. In this Chapter, the focus is mainly on solids under tensile, compressive and shear stresses. We consider the medium of interest to be a solid homogeneous continuum in order to attribute it average mechanical properties.

Stress and strain

The first notions that need to be introduced here are the notions of stress and strain. We consider a force, noted P , applied on a surface. Cartesian coordinates are chosen so that the normal to the material's surface corresponds to the axis z. This force is represented in this space by a vector which norm represents the amount of force applied on the surface and its direction the direction of the force. This vector can be described as the sum of a normal component P z and two tangential (or shear) components P x and P y in regards to the loaded surface. The tensile or compressive stress, noted σ z here, is defined as the force applied perpendicularly to the surface of the material per unit area [5]:

σ z = dP z dA (3.1)
with P z the normal force applied on a surface A of a cylindrical free material.

If the stress is assumed to be uniformly distributed, Eq. 3.1 can be written:

σ z = P z A (3.2)
The shear stress can be described similarly by using the tangential components of the applied force.

As for the strain, it defines the deformation of a material subjected to a loading. The strain, noted ϵ, is defined as a displacement per unit distance in a given direction. When an sample is deformed due to a force, it can contract and elongate along different directions in regards to the loading it undergoes. Axial strain considers deformations in the direction of loading and lateral strain in the transverse direction. The true (or logarithmic) strain can be defined as the elementary change in dimension (δl) of a sample relative to its length (l) in the considered direction:

ϵ = δl l (3.3)
Integration of Eq. 3.3 over an arbitrary length variation ∆l provides the following expression:

ϵ = ln l 0 + ∆l l 0 (3.4)
with l 0 the initial length and ∆l the variation in length.

For small strains (ϵ < 0.1), the following approximation can be found:

ϵ = ∆l l 0 (3.5)
The relation between stress and strain in a material is highly characteristic of its mechanical behavior. The stress-strain curve is a graphic representation of the behavior of a material, often under tensile or compressive loads, as illustrated in Fig. 3.2. This curve highlights the behavior of a material depending on the amount of force applied. In most biological tissues and under low forces, a linear behavior is observed. However, beyond a certain amount of stress, the stress-strain curve cease to be linear. This point in the stress-strain curve is the proportional limit. If the load increases even more and the elastic limit is exceeded, the material won't return to its initial shape if unloaded. The resulting irreversible change in the material is called plastic deformation. Beyond this point a stronger curvature of the curve is observed and the material ends up rupturing [6].

Elastic and Hookean materials

In order to later introduce some necessary definitions, notions of elastic and Hookean materials should be first defined. An elastic material is defined as a material that goes back to its original shape when unloaded. If the stress in a material is directly proportional to the strain for strains up to the elastic limit, the material can be called Hookean or linear elastic [5], [7]. This is illustrated in Fig. 3.3. When the loads applied on a linear elastic material are relatively small, as it is the case in elastography, knowledge of the full stress-strain curve is not necessary. It can be considered in this case that the material's behavior is described by the linear part of the stress-strain curve. In this region, the curve is fully characterized by its slope. In the case of tensile or compressive stress, this slope corresponds to the elastic modulus or Young's modulus, noted E. Its definition is indeed:

σ z = Eϵ z (3.6)
This definition is known as Hooke's law and is later generalized using tensor notation.

The previous definitions were provided based on tensile or compressive stress, meaning that the component of the stress normal to the surface were considered. If the tangential component of the stress is considered, the shear modulus, noted G or µ, can be defined as the shear stress to lateral strain ratio below the proportional limit in an elastic material:

σ x = µϵ x (3.7)
When a tensile load acts on an elastic isotropic material, it elongates axially and contracts in the two transverse directions. On the other hand, under a compressive load, the material contracts axially and elongates transversely as illustrated in Fig. 3.1. It can be observed that the axial and transverse strains are proportional. The ratio between the two strains is known as Poisson's ration (ν):

ν = - ϵ x ϵ z = - ϵ y ϵ z (3.8)
with ϵ z the axial strain and ϵ x the lateral strain.

Time and history-dependent mechanical properties

The previous notions have laid the foundations for characterizing mechanical properties of solids based on the stress-strain relationship. However the complexity found in biological tissues makes it necessary to introduce more notions to describe more accurately their behavior. The time and history dependence of mechanical properties is notably examined here.

A significant notion for characterizing soft biological tissues is the notion of viscosity and viscoelasticity. A viscous material's behavior is time-and history-dependent. Viscosity is defined in liquids as a resistance to flow. Viscoelasticity defines a material which is both viscous and elastic. A viscoelastic material can be characterized by different physical quantities such as stress-relaxation, creep, strain-rate sensitivity and hysteresis. Stress or force-relaxation describes the decrease of stress over time in a material that is continually stretched at a fixed length. Creep describes the elongation of a material over time when subjected to a constant amount of stress. Retardation time is the rate at which the elastic deformation takes place in a material creeping under a constant load. Strain-rate sensitivity measures the differences in stress-strain behavior under different strain rates. At last, hysteresis can be defined as the history-dependence of the material's loading on its mechanical properties.

To simplify these considerations, a linear viscoelastic behavior is often considered in models. It considered when the stress to strain ratio in a material is a function of time alone [8].

In order to investigate time-dependent properties of a material and characterize viscoelastic behaviors, dynamic loading can be realized experimentally. The phase and amplitude difference between the loading signal and the recorded response indicates the behavior of the material. A phase lag of zero corresponds to a purely elastic material and a phase lag of 90°to a purely viscous behavior. The notion of complex shear or elastic modulus, noted G * or E * is introduced to characterize this behavior. The complex modulus characterizes the complex stress-strain relation in viscoelastic materials under an periodic strain [9], [10]. Elastic or storage modulus (the real part of G * or E * ) corresponds to the portion of the total rigidity attributable to elastic behavior. Viscous modulus or loss modulus (the imaginary part of G * or E * ) corresponds to the portion of the total rigidity attributable to viscous behavior.

Direction and spatially-dependant mechanical properties

The last important considerations to be made in defining mechanical properties in a material is its homogeneity and anisotropy. It was indeed assumed previously that the material was homogeneous and isotropic to define its mechanical properties. As it is described later, it is rarely the case in biological tissue elastography.

Anisotropy defines a direction-dependent behavior. The orientation of the considered medium in regards to the applied constraint needs to be taken into account in describing an anisotropic mechanical behavior. Anisotropy often emerges from the material's layered or fibrous structure. Anisotropy of mechanical behavior can be accounted for by using direction-dependent mechanical properties in the modeling of their behavior.

Homogeneity is assumed in most mechanical testing methods and as well as in models.

Homogeneity is necessary, at least at reasonable scale, to define mechanical properties. Heterogeneity and discontinuity in mechanical properties is a major challenge in stiffness imaging. Heterogeneity is indeed the source of contrast in imaging but can also challenge the assumptions of the models used.

Biological Tissues Models

The type of biological tissue and the scale considered greatly affect the model considered for mechanical behavior. Majors differences in mechanical properties allow the separation of materials into groups such as solids or liquids. As for solids, which we are focusing on, viscoelastic behavior is mainly observed in biological tissues. Living biological tissues are often considered incompressible, which means that their volumes remains unchanged when deformed. When considering soft biological tissues, the mass density is generally considered as identical to water (ρ = 1.0g/cm -3 ) [11].

Biological tissues are generally viscoelastic, anisotropic and non-homogeneous. They are by design hierarchically structured: changes in dimensional scale entails different mechanisms of deformation. Biological tissues can be described by their structural and material properties. Structural properties characterize the tissue in its intact form whereas material properties characterize the behavior of the material in the tissue [7].

The intricate structure and varied composition of biological tissues lead to complex mechanical behaviors. The models described in this section depict a relatively simple approach to the sometimes complex mechanical behaviors of biological tissues. It is however important to note that more refined and specific descriptions have been made in literature. A number of sources in literature have described in details these structures and mechanical behaviors for different tissues categories [10] [12]. According to the classification of materials established in [13], biological tissues can be described as elastomer (muscles, skin, soft tissues), polymer or polymer composite (ligaments and tendons) or ceramic or ceramic composite (bones and teeth) [14].

Considering biological tissues' properties, two main models are be considered here: linear elastic behavior and linear viscoelastic behavior. Two types of materials are described: isotropic medium and transverse isotropic medium. These models are the most common for elastography experimental studies.

Purely linear elastic models

When observing a truly elastic behavior in an homogeneous unbound solid, Hooke's law characterizes the linear relation between the mechanical constraint applied on a medium and the resulting deformation. Hooke's generalized law can be written as:

σ i,j = c i,j,k,l ϵ k,l (3.9)
with σ i,j the stress tensor, ϵ k,l the strain tensor and c i,j,k,l the proportionality coefficient. the indices i, j, k, l represents the three spatial dimensions considered.

Due to the symmetry of the stress and strain tensors, behind the 81 components of the proportionality tensor, 36 elements are to be determined for a given material. A 6 by 6 matrix c α,β is consequently often used to represent this tensor. The indexes α and β represent the possible values of (i,j) and (k,l) respectively. Hooke's law can consequently be written using matrices:

σ α = c α,β ϵ β with σ α =         σ i=1,j=1 σ i=2,j=2 σ i=3,j=3 σ i=2,j=3 σ i=1,j=3 σ i=1,j=2         and ϵ β =         ϵ k=1,l=1 ϵ k=2,l=2 ϵ k=3,l=3 ϵ k=2,l=3 ϵ k=1,l=3 ϵ k=1,l=2         (3.10)
The isotropic case. When considering a tissue is isotropic, it means its behavior and properties are the same regardless of the considered axis. The elastic tensor defined previously can consequently be expressed using only two variables λ and µ, called Lamé coefficients. Its corresponding matrix is written:

c α,β =         λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ 0 0 0 0 0 0 µ         (3.11)
Under the assumption of mechanical isotropy, Lamé coefficients λ and µ can be expressed using Young's modulus and Poisson's coefficient:

λ = Eν (1+ν)(1-2ν) µ = E 2(1+ν) (3.12)
Using Eq. 3.9 and 3.11, we can express the stress-strain relation using Lamé coefficients:

σ = λT r(ϵ)I + 2µϵ (3.13)
with Tr the trace and I the identity tensor.

In the case of a linear-elastic and isotropic medium, σ can be expressed as:

σ ij = 2µϵ ij + λϵ kk δ ij (3.14)
with δ ij the Kronecker delta and λ and µ Lamé coefficients.

The transverse isotropic case. In the case of anisotropic soft biological tissues, we can assume a symmetry axis along the direction of the fibers. Under this assumption, the medium can be considered as transverse isotropic (equivalent to an hexagonal system) [15]. Using Voigt's form and the transverse isotropic assumption, the stiffness tensor can be expressed using 5 different elastic constant:

C =         c 11 c 11 -2c 66 c 13 0 0 0 c 11 -2c 66 c 11 c 13 0 0 0 c 13 c 13 c 33 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 66         (3.15)
The transverse shear wave speed (c ⊥ ) can be defined when the direction of propagation is perpendicular to the fiber's axis and the parallel shear wave speed (c ∥ ) where the direction of propagation is along the fibers. The coefficient of the stiffness matrix can be expressed as [16]:

c 44 = ρc 2 ∥ c 66 = ρc 2 ⊥ (3.16)
Consequently, when considering the surface of a transverse isotropic sample, the local stiffness can be described as a function of the angle between the considered direction of stiffness evaluation and the local main fiber orientation (θ). If we suppose no coupling between the transverse and longitudinal propagation directions, it can be considered the following expression for Young's modulus :

E(θ) = 3ρ(c 2 ⊥ sin 2 θ + c ∥ cos 2 (θ)) (3.17)

Linear viscoelastic models

In viscoelastic models, the relation between stress and strain is more complex than for purely elastic behavior. Depending on the model, they can involve the creep and relaxation function as well as elastic moduli at different times and stress/strain values. Only the case of linear viscoelasticity is considered here.

Voigt and Maxwell models. Viscoelastic behavior can be also modeled by the combination of springs and dampers. Springs represent the elastic behavior, where the energy from the load applied on the material is entirely stored in the compressed spring as potential energy. On the other hand, dampers model the viscous behaviour as they represent the dissipation of energy with viscous friction. Voigt (or Kelvin-Voigt) model and Maxwell model are both used to represent the behavior of viscoelastic materials using a spring and a damper. In Voigt model,both element are connected in parallel whereas in Maxwell model, they are in series, as illustrated Fig. 3.5. This gives the following expressions for the stress:

σ voigt = Eϵ + ηdϵ (3.18) σ maxwell + η E dσ maxwell = ηdϵ (3.19)
These models can be generalized in the three dimensions of space for anisotropic materials by using tensors and combined with Newton's law of motion.

Theory of mechanical wave propagation

This section aims at linking the previously defined mechanical properties with measurable physical quantities, notably the propagation speed of mechanical waves. This is the foundation for the elastography approaches presented in a following section. Based on previously defined biological tissues models, this paragraph introduces the theory of mechanical wave propagation in such materials.

Basic relations govern displacement and motion of tissues. These basic relations are summarized in this section. An elastic medium defines a medium's capacity to go back to its initial state after the external constraints that deformed it are suppressed. This behavior is due to internal constraints. At equilibrium, those internal constraints are negated by the external forces. Using continuum mechanics methods for a deformable medium, conservation of momentum is expressed as [17]:

d dt V ρ du dt dV = S T (n) dS + V ρb dV (3.20)
with ρ the medium's density, u the displacement vector, b the body force per unit mass vector and T (n) = σ • n the traction vector on the surface S on the volume V in the normal direction n .

Using the divergence theorem, the differential form of Eq. 3.20 can be expressed [18]:

ρu = ∇ • σ + ρb (3.21)
Considering small deformations and neglecting the effect of body forces on the measured material's response, the deformation (∂u) can be expressed in function of the strain sensor ϵ ij :

ϵ ij = ∂u i ∂x j + ∂u j ∂x i 2 (3.22)
Combining expressions 3.21 and 3.22, an equation governing wave propagation in media can be obtained. Most often, assumptions such as the ones are made to simplify its terms. The case of elastic and viscoelastic behaviors for isotropic media are described here.

Homogeneous linear-elastic isotropic medium

Combining expressions 3.20, 3.21 and 3.14 for a homogeneous medium, the general dynamic response of an isotropic and linear-elastic medium can be written:

∇ • σ = λ∇ • (Tr(ϵ)I) + 2µ∇ • ϵ = ρ ∂ 2 u ∂t 2 (3.23)
By expressing the strain sensor ϵ as the symmetric of the gradient of the displacement grad(u), we obtain:

(λ + µ)∇(∇ • u) = ∇ • (T r(ϵ)I) + 2∇ • ϵ = ρ ∂ 2 u ∂t 2 (3.24)
By decomposing the vector u in transverse and axial displacements in respect to the propagation axis (respectively u S and u P ), it can be written:

u = u P + u S with ∇ × u P = 0 and ∇ • u S = 0 (3.25)
Injecting Eq. 3.25 in Eq. 3.23, two separate equations are obtained:

(λ + µ)∆u P = ρ ∂ 2 u P ∂t 2 µ∆u S = ρ ∂ 2 u S ∂t 2 (3.26)
From these relations, the speed of compression waves (corresponding to axial displacements) and shear-waves (corresponding to transverse displacements) can be expressed:

c 2 P = (λ+2µ) ρ c 2 S = µ ρ (3.27)
These expressions are the link between wave propagation (that can be recorded experimentally) to quantitative stiffness estimators characteristic of the material.

Homogeneous viscoelastic isotropic medium

A more complex description is based on a viscoelastic model. Considering viscous behavior as a loss of energy in viscous friction, the partial differential equation (PDE) governing mechanical waves in a viscoelastic isotropic material is:

ρ ∂ 2 ∂t 2 U -γ ∂ ∂t U = E 2(1 + ν) ∆U + E 2(1 + ν)(1 -2ν) ∇∇ • U (3.28)
with γ the attenuation coefficient and U the complex displacement U = u(x)exp(iwt).

Considering a sinusoidal excitation, the expression Eq. 3.28 can be written:

(-ρw 2 -iwγ)u = E 2(1 + ν) ∆u + E 2(1 + ν)(1 -2ν) ∇∇ • u (3.29)
In this equation, Young's modulus and Poisson's ratio are directly linked to displacements. From this equation, the attenuation coefficient can also be calculated to estimated the viscosity. § 3.

Elastography Modalities for biological tissues

The term elastography was introduced in 1991 [19] and the first review article on the subject followed rapidly [20]. However, first studies of biological tissue's response to mechanical vibration dates back from the 1950s using strobe light and photography [21] and first non-mechanical quantitative stiffness measurements were reported in 1987 [22] using an ultrasonic probe. In order to determine the mechanical properties of a sample, different methods can be used but they are all based on the same principle: a stress is applied on the sample and its response is then recorded. From this recording, mechanical properties can be extracted, usually with prior knowledge of the medium's density and geometry.

Elastography modalities include imaging systems able to provide information on the mechanical properties of a sample. Mechanical testing such as tensile tests or compressive tests are using force transducers and an actuator to provide structural mechanical properties of tissues. Mechanical testing is a useful tool for establishing models and providing measurements comparison with elastography but is not included in elastography modalities since it does not provide images. Mechanical testing of biological tissues is rarely done in-vivo since it often requires homogeneity of a specific sample geometry. The use of tonometer probes is the main mechanical testing method in medicine. It relies on different principles: applanation, indentation, dynamic contour or rebound [23]. Tonometer usage has been mainly reported to evaluate intraocular pressure (IOP) [24] but some studies also mention arterial pressure monitoring [24] or brain stiffness measurements [25]. Some other techniques that have been utilized in medicine include transducers, pressure sensors and strain gauges to monitor tendons, ligaments and joints.

Ultrasound, MRI and optic-based elastography have also been developed with different modalities of testing, yielding different values for the same mechanical parameter. The experimental set-up used for testing mechanical properties are indeed based on different amounts of stress and different frequencies. It was chosen to classify theses modalities according to the stress used for testing: compression, harmonic or transient stress.

Compression, static and quasi-static elastography

Principle

Static and quasi-static elastography methods consist in using a static or a low-frequency sinusoidal mechanical loading. An image of the sample is taken before and after the loading. The under the assumption of homogeneous loading, the strain image can then be computed between the two images. In quasi-static methods, additional frequencydependent viscous information can be provided using Eq. 3.29. This technique is especially well-adapted to in-depth imaging since compression on the tissues can be achieved easily in the imaging plane. This technique was first adapted to ultrasound imaging [19] but MRE and optical setups quickly followed [1], [26]. An illustration of compression elastography is presented Fig. 3.6

Advantages and limitations

Static and quasi-static methods are easy to implement and provide high-resolution images. Indeed, in the absence of artifacts and noise, the resolution of compression elastography is theoretically the same than the imaging system's. The first tests were achieved using a handheld ultrasound scanning transducer pressed against the region of interest. Clinical results in breast cancer on margin assessment or screening systems show promising results [28]. Different forces can be applied using compression elastography: with different amplitude and different frequencies in quasi-static elastography. This allows investigation of the stress-strain curve and viscous behavior of the material, which is especially relevant for biological tissues.

Static and quasi-static elastography don't achieve quantitative stiffness images because of unknown non homogeneous loading and complicated boundary conditions, as demonstrated for OCE [29]. Yet, quantitative compression elastography has been demonstrated [30] by quantifying of the stress distribution using a reference sample for calibration, measuring the surface distribution of compression force using a silicone sensor between the sample and the compression plate [31], using a pressure chamber [32] to apply uniform loading or using inverse problem methods [33]. Quantification is however challenging and assuming uniform or known depth distribution of the stress within heterogeneous tissues under unknown boundary conditions in and out of the imaging plane makes compression elastography contrast still considered qualitative in the main elastography reviews in ultrasounds [34], optics [35] and MRI [36].

Another possible limitation of compression elastography is the need of contact due to the generation of compressive stress, which limits the potential applications. The use of Acoustic Radiation Force (ARF) to compress tissues uniformly has been demonstrated recently [37] but no in-vivo experiment was performed yet. The elongation of the material under compression can also introduce biases in image reconstruction and compression elastography using an external compression source is limited in depth.

Harmonic elastography

Principle

Vibration sonoelastography.

Vibration sonoelastography describes techniques utilizing the sample's response to lowfrequency vibrations (20-1000 Hz) using ultrasonic imaging. The response in amplitude, frequency or phase can be used to determine viscoelastic properties. Vibration sonoelastography is the oldest elastography modality and the first study of biological tissues response to vibration was published as early as the 1950s [21]. The first sonoelastography system is reported in 1987 [38] based on pulsed Doppler detection. The contrast and values provided by vibration amplitude sonoelastography are qualitative and require at least a calibration step. The development of phase gradient sonoelastography, where the phase and amplitude of the sample's mechanical response is recorded, has allowed for quantitative measurements [17] based on Eq. 3.29.

MRI-and optic-based harmonic elastography.

Harmonic elastography also consists in using an oscillating sinusoidal excitation and recording the frequency-dependence of the recorded material signal, except using optical imaging or MRI. Harmonic elastography relies on the same principle as sonoelastography, using local frequency estimations or direct or iterative inversion algorithm to retrieve viscoelastic properties. Harmonic elastography was introduced for MRI in 2000 [39], [40] and is now the main approach in MRE [41]. Harmonic elastography has also been reported on Rayleigh waves in optical elastography [42], [43], using different frequency to probe the sample in depth.

Advantages and limitations

Harmonic elastography shares a lot of the advantages and drawbacks of compression elastography. The viscoelastic behavior of tissues can also be accessed but at higher mechanical frequencies. Harmonic elastography usually require contact with a mechanical transducer to generate sinusoidal oscillations within the tissues, which might limit some biomedical applications. The boundary conditions of heterogeneous tissues may cause modal patterns and losses due to reflections which hinders image quality. Contact is also required by the vibration source and the imaging transducer for ultrasound imaging. The principal advantage of harmonic elastography in regards to static or quasi-static elastography is that it is experimentally easier to generate a sinusoidal shear-wave field in depth that it is to generate compression. MRE is principally developed based on harmonic excitation signals and successfully achieves elastography tomography in full slices of the body [44].

Transient elastography

Principle

Transient elastography was first introduced in MRE [45]. It is based on direct imaging of the mechanical wave speed in the medium. Considering the shear and compressionwave speed ranges in biological tissues (1-10 m/s and 1500 m/s, respectively), transient elastography is often performed using shear-waves. The straightforward link between shear-wave speed and elastic modulus in isotropic linear elastic media presented in Eq. 3.27 is promising for quantitative stiffness measurements without the need for calibration, inversion algorithms or sensors. Classic transient elastography is performed by sending an impulse through the tissue and using an imaging system to track its propagation. The speed-map is then computed using inverse problems [46], [47] or time-of-flight [42], [43], [48]- [57]. Transient elastography is the base of most commercial ultrasound elastography systems : Logiq E9 by GE Healthcare, Fibroscan TM by Echosens, Virtual Touch TM by Siemens Healthineers, Aixplorer TM by Supersonic Imagine or ElastPQ TM by Philips.

Advantages and limitations

Since transient elastography relies on the local and instantaneous measure of the propagation speed of shear-waves, it is less dependent on the sample's geometry than harmonic, static and quasi-static methods. The direct link between shear wave speed and shearmodulus also suggest that no prior knowledge on the shear-wave amplitude is necessary to achieve quantitative stiffness estimation. As transient elastography relies on shearwave speed estimations, the resolution in transient elastography is not as good as with static or quasi-static approaches. The problem of resolution in transient elastography is a complex problem. A recent study has yet found that the resolution in shear transient elastography is not limited by the mechanical wavelength [58] and often depends on the algorithm used for speed processing and the imaging system's resolution and sensitivity.

In classic transient elastography, prior knowledge of propagation direction or filtering methods [59] is often necessary when performing transient elastography. This makes it especially sensitive to local changes in the propagation direction caused by refraction or reflections.

Transient elastography is particularly developed in ultrasound imaging but a recent increase in research published in optical transient elastography provides very promising methods for high-resolution imaging. Methods utilizing contact-less mechanical stimulation based on ARF, air-pulse [49], [51], [57] or heartbeats [60] are especially interesting for ophthalmology application and intraoperative imaging. Lately, passive methods based on noise-fields have been introduced to ultrasounds [61] and then to optical elastography [4], [62], [63] to overcome the limitations or a controlled shear-wave source. Passive elastography is described in more details in Chapters 5 and 6. § 3.4 Motivations for biological tissues' mechanical properties evaluation

The structure of biological tissues can be altered by certain pathologies. This change of structure leads to the modification of the mechanical properties of tissues whether it is local or affects the whole organ. A vast range of scientific studies have found correlations between the tissue's stiffness and an abnormal function in various cases: malignant tumors [64], metastasis [65] or inflammatory diseases [66]. When the tissues are accessible, experienced physicians can detect those structural anomalies through palpation, which provides only imprecise and qualitative information. Having quantitative stiffness contrast on medical imaging systems could be a strong aid in diagnosis, particularly in cases of small lesions or inaccessible tissues. For MRE, the first commercial system was marketed by GE Healthcare in 2010. The market for elastography systems in medicine is growing due to its novelty and the promising applications of elastography. This includes diagnostics in ophthalmology, liver diseases and breast tumor screening, which are described here. The recent development of optical elastography and possible applications for intraoperative imaging or endoscopic imaging in minimally invasive surgery are also explored in this section. The eye's main function is to refract light to the back of the eye and convert it into electrical signals that are sent to the brain through the optic nerve to be interpreted and understood.

Ophthalmology applications

The cornea is responsible for refracting light by providing a rigid protection and maintaining the shape of the eye.

Corneal structural abnormalities such as inflammation or lesions can significantly reduce visual acuity.

The study of corneal biomedical properties has significant importance in ophthalmology as the properties of pathological corneas differ from that of healthy ones. Corneal stiffness, corneal thickness and intraocular pressure (IOP) measurements have been linked to a number of ocular diseases such as glaucoma, age-related macular degeneration, presbyopia or keratoconus [2], [68]. Developing imaging techniques for the mechanical properties of ocular structure could be a great help for diagnosis but also for monitoring refractive eye surgeries. Corneal biomechanical properties are indeed not often taken into account when performing refractive surgery, which can cause residual refractive errors due to inter-individual variability [69], [70]. Refractive defects, glaucoma and age-related macular degeneration are among the most common eye disorders. For these reasons, imaging the biomechanical properties of cornea could be a great tool in clinical ophthalmology applications, whether it is for diagnostic purposes or for monitoring the effects of drug treatments or surgeries.

IOP is defined as the pressure of aqueous humor inside the eye. IOP is clinically measured using tonometers. However, the measured IOP with such devices has been shown as highly correlated to other factors such as corneal geometry: corneal stiffness [71], thickness [72], but also astigmatism or corneal curvature can lead to errors in IOP measurements [73]. The introduction of tonometers using estimations of corneal hysteresis and corneal resistance factor to correct the influence of corneal biomechanical on IOP have recently allowed more accurate measurements of IOP [74]. Elastography using is a very promising technology for corneal stiffness imaging as it allows non-invasive accurate quantitative images of the elastic structural properties in-situ where tonometers only allowed for one measurement. Cornea elastography was mainly investigated using high frequency ultrasound systems (50-100 MHz) or OCT.

The emergence of high-frequency ultrasound imaging in the 1990s led to the development of ultrasound microscopy [75]. In the following decade, elastography images of the cornea were performed ex-vivo [76] and in-vivo [77], [78]. Ultrasonic systems using spectrometry-based [79] or time of flight-based [80] methods were also reported to perform mechanical properties measurements but these methods did not allow full-field images or depth resolution. ) of IOP using an a phase-sensitive OCT system and non-contact stimulation based on ARF, published in [55] OCT was developed for ophthalmology in the 1990s as it is particularly well-adapted to mostly transparent ocular tissues. OCT-based elastography offers a promising alternative to ultrasonic elastography, allowing for better resolution and potentially more comfortable non-contact imaging. Optical imaging is indeed more sensitive than ultrasound imaging and can achieve stiffness estimation with minimal mechanical constraints. Elastography using OCT on the cornea was performed successfully based on ultrasound methods [81]. Recent development of quantitative non-contact OCT elastography systems using air-puffs or acoustic radiation force (ARF) have shown particularly great potential for clinical application [57].

Fibrosis and cirrhosis diagnostic and monitoring

The earliest applications of ultrasound elastography was for liver fibrosis and cirrhosis diagnostics. Fibrosis can be defined as the formation of permanent scars on an organ, consequently modifying its structure and interfering with its function. Cirrhosis is the most advanced stage of liver fibrosis. Both diagnostics are characterized by a stiffening of the liver tissues [3] and can be classified using a scoring system. Fibrosis and cirrhosis are common and growing diagnostics. The liver is however inaccessible to palpation and liver biopsy, although considered as the gold standard, is invasive and flawed. Indeed, the scoring-based evaluation of the degree of fibrosis of a sample has shown biases due to the sample's localization and size and to intra-and inter-observer variability [82] Ultrasound imaging, CT scans and MRI have been investigated as an alternative to liver biopsy for cirrhosis diagnostic. Because of the higher-speed, lower cost and lower expertise required, ultrasound imaging is considered in this case more adapted for non-invasive liver imaging. Morphological ultrasound, Doppler or contrast-agentbased imaging have been tested for cirrhosis diagnostic. Ultrasound-based elastography have shown good sensitivity and specificity [83]. Ultrasonic elastography based on transient, ARF and supersonic shear-wave have been investigated and comparison between the techniques showed similar accuracy [84]. Elastography was performed clinically on commercial systems such as the Fibroscan (Echosens, Paris, France) since 2003 based on transient elastography [START_REF] Sandrin | Transient elastography: A new noninvasive method for assessment of hepatic fibrosis[END_REF] and the Aixplorer (SuperSonic Imagine, Aix-en-Provence, France) since 2008 based on ARF and supersonic shear-waves imaging [START_REF] Dhyani | Validation of shear wave elastography cutoff values on the supersonic aixplorer for practical clinical use in liver fibrosis staging[END_REF]. Even though these systems have been demonstrated to provide reliable non-invasive liver assessment, the penetration depth limitation of ultrasound imaging makes this image technique less reliable in some cases, such as ascites and obesity. In these cases, Magnetic Resonance Elastography (MRE) has shown significantly reliable results [START_REF] Chen | Mr elastography of liver disease: State of the art[END_REF], [START_REF] Yin | Ultrasound or mr elastography of liver: Which one shall i use?[END_REF] as illustrated Fig. 3.10.

Characterization of structural anomalies for breast cancer diagnostic

Manual palpation has been used routinely by physicians and patients to screen the breast for variations in the tissue's elastic properties. Palpation is a cost-effective and highly practical tool but is limited to accessible tissues and sizeable lesions [START_REF] Hill | Self examination of the breast: Is it beneficial? meta-analysis of studies investigating breast self examination and extent of disease in patients with breast cancer[END_REF]. The practice of routine mammographic screening have been a great help in early detection of breast cancer and has decrease its mortality [START_REF] Kerlikowske | Efficacy of screening mammography: A meta-analysis[END_REF]. Mammography is however limited in the detection of early-stage or node-negative cancers, especially in higher density breasts. When an abnormality is detected through a mammography, the tumors should be classified using a grading systems in order to provide a prognosis and treatment to the patient. Tumor grade is a classification based on the degree of differentiation of the tumor tissue and evaluated invasively by performing a biopsy.

Quantification of breast-tissue stiffness assessment has shown great contrast for tumor localization and characterisation with up to a one to ten ratio between Young's modulus of healthy breast tissue and malignant tumors using mechanical testing [START_REF] Samani | Measuring the elastic modulus of ex vivo small tissue samples[END_REF]. This high variability of Young's modulus in soft biological tissues makes elastography a promising tool for non-invasive breast screening and could be a great non-invasive indicator of the tumor's grade. Ultrasound elastography imaging was performed on breast tissues with promising contrasts [1], [START_REF] Nightingale | Shear-wave generation using acoustic radiation force: In vivo and ex vivo results[END_REF], [START_REF] Wojcinski | Multicenter study of ultrasound realtime tissue elastography in 779 cases for the assessment of breast lesions: Improved diagnostic performance by combining the bi-rads®-us classification system with sonoelastography[END_REF]. Breast cancer is the most common type of cancer worldwide, which explains the vast amount of clinical research but similarly to the case of breast tissues, elastography images of human specimens have outlined successfully a number of pathological tissues such as thyroid [START_REF] Dighe | Differential diagnosis of thyroid nodules with us elastography using carotid artery pulsation[END_REF], prostate [START_REF] Castaneda | Performance of three-dimensional sonoelastography in prostate cancer detection: A comparison between ex vivo and in vivo experiments[END_REF] or liver tumors [64]. These studies have mainly been performed using MRE or ultrasound elastography and demonstrate promising results for medical imaging.

Detection of structural anomalies for intraoperative margin assessment

Imaging techniques such as MRI or metabolic imaging using nuclear medicine are able to localize, identify and monitor tumors in a large number of cases. These systems, although they demonstrated impressive results, are unable to perform real-time acquisition and not well-adapted to surgical imaging. X-rays and ultrasounds are used in surgery as imaging techniques but the low contrasts of radiology for soft tissues and the required contact of ultrasound probes limit their use.

During tumor removal surgery, it is crucial to remove the cancerous lesion with a small amount of normal tissues surrounding it, called the margin. The margin ensures that no cancerous cell was left behind. However the delineation between healthy tissues and lesions is not obvious in all cases. Surgeons are usually mostly guided by their sight and touch to judge for tumor margins. This subjective assessment is backed by ulterior analyses. The gold standard for intraoperative margin assessments is a cytology or histology analysis of the sample, illustrated Fig. 3.12. This analysis is yet time and resource consuming and negative results often implies the need for a second surgery weeks later [START_REF] Pradipta | Emerging technologies for real-time intraoperative margin assessment in future breast-conserving surgery[END_REF]. Probes have been developed as well as artificial intelligence to perform fast histology analysis. Optical imaging techniques are also emerging in surgical fields to provide realtime non-contact feedback during tumor resection surgeries, based on spectral, OCT or fluorescence imaging [START_REF] Brown | Optical assesssment of tumor resection margins in the breast[END_REF], [START_REF] Boer | Optical innovations in surgery[END_REF]. Elastography has been recently introduced for margin assessments using OCT-based elastography. Studies have notably demonstrated that abnormalities and heterogeneity in mechanical contrast reflected the presence of a tumor in breast [START_REF] Allen | Wide-field optical coherence microelastography for intraoperative assessment of human breast cancer margins[END_REF] and prostate biopsies [START_REF] Li | Detection and characterisation of biopsy tissue using quantitative optical coherence elastography (oce) in men with suspected prostate cancer[END_REF]. This is illustrated Fig. 3.13 where the malignant tissue can be distinguished from the stroma. Elastography images showed a increased contrast between healthy and diseased tissues. It has been observed that use of elastography greatly improved the sensitivity and specificity of margin assessments compared to existing methods such as OCT [START_REF] Kennedy | Diagnostic accuracy of quantitative micro-elastography for margin assessment in breast-conserving surgery[END_REF]. Elastography appears in these studies as a promising technique to help surgeons in differentiating malignant from non-malignant tissue biopsies during surgery and consequently reduce the incidence of residual tumor remaining after a tumor-removal surgery.

Detection of structural anomalies for endoscopic systems

An endoscope is an imaging system mounted on a thin metal tube, allowing it to be slid inside the body for diagnostic purposes. Laparoscopy is type of minimally invasive surgery using only small incisions on the patient. By sliding surgical instruments and an endoscope through these incisions, some surgeries can be performed with minimal scaring. Since the 1990s, laparoscopy has become more and more common for a growing number of surgeries on the abdomen or pelvis [START_REF] Lau | History of endoscopic and laparoscopic surgery[END_REF]. Even if laparoscopy requires specifically trained physicians, it displays major advantages compared to open-body surgery, notably a decrease risk of nosocomial infections, shorter hospital stays and better patient remission.

The recent combination of endoscopic OCT imaging systems has shown promising results by providing high resolution imaging of inaccessible biological tissues [START_REF] Tearney | In vivo endoscopic optical biopsy with optical coherence tomography[END_REF]. An example of a steerable OCT endoscopic system is presented in Fig. 3.14. The adaptation of optical coherent elastography to endoscopic systems [START_REF] Chau | Mechanical analysis of atherosclerotic plaques based on optical coherence tomography[END_REF], [START_REF] Frolov | Compression elastography and endoscopic optical coherence tomography for biomechanical properties evaluation of cerebral arteries walls with aneurysm and their phantoms[END_REF] are promising for diagnostics and laparoscopic surgery guidance.

Conclusion

Elastography contrast can be instrumental in helping physicians for a range of challenges in medicine, addressing issues in diagnostics as well as in surgery guidance. Motivations for developing elastography in biomedical imaging have been introduced in this Chapter, highlighting the role of different imaging modalities in elastography. The basic notions and the principle of elastography have been presented as well as brief state of the art in terms of biological tissue elastography, including ultrasounds, MRI and optic-based elastography approaches. These first two chapters aimed at laying the foundation and context for presenting the experimental work of this thesis in the next chapters. The digital holography elastography system and the quantitative stiffness images it provides are presented in the following Chapters.

-Chapter 4 -Experimental set-up and control measurements § 4.1 Summary With recent developments in ultrasound elastography, physicians are now provided with commercial systems providing a quantitative and non-invasive stiffness images. If ultrasound systems are very useful in diagnostics, their measurements require contact which is limiting in a number of cases, such as surgery and delicate tissues. The introduction of optical methods for elastography in the late 1990s has allowed contactless and precise measurements, which are especially promising for applications to ophthalmic imaging or surgery imaging.

Full field elastography methods using optics have been researched, using controlled sources of mechanical perturbations. Oldenburg's team [1], [2] achieved full field depth elastography with digital holography using harmonic perturbations. A. Nahas et al. [3] demonstrated the feasibility of full-field transient elastography using full-field Optical Coherence Elastography to track shear wave propagation in real time.

The objective of this thesis is to design an experimental imaging system to perform full-field stiffness images. By choosing optical imaging, a fast, non-contact and high-resolution imaging system is built. Based on interferometry, the set-up achieves high sensitivity displacement imaging. The resulting images can be exploited to perform shear-wave elastography imaging with low-power mechanical excitation. The optical elastography system is presented and characterized as well as the different elastography approaches performed in the first sections.

In order to ensure the quantitativity of the developed approaches, results should be compared. However, the excitation frequencies used in mechanical testing (≈ 1 Hz) and commercial ultrasound systems (≈ 50Hz) are very different from the ones used here (≈ 1 kHz). A previous study [4] had indeed compared elasticity measures from mechanical testing (0.1 to 4 Hz) and ultrasound elastography on the liver. This study concluded that no direct and quantitative comparison could be made between both measurements techniques because of the different mechanical excitation frequency ranges. Control measurements were consequently performed with the experimental set-up. Harmonic and classic transient elastography approaches were used here for validation of the algorithms and control measurements of the manufactured samples. The focus of this Chapter is on control measurements performed using harmonic and transient elastography based on time-of-flight. These control measurements were used to set the benchmark for the new 99 approaches introduced in Chapters 5 and 6 based on noise-correlation.

Testing of the experimental set-up was performed on simulation images, polymer samples and ex-vivo biological tissues. The simulation and the homemade tissue-mimicking samples protocol are also described. The results of control measurements are presented in the last section of this Chapter § 4.2 Optical experimental Set-up

The elastography set-up was developed in order to perform single-shot full-field quantitative phase images. This set-up allowed for different types of mechanical shear-wave sources to be used: impulses, harmonic signals, chirps or noise. The objective of the developed design is to provide slow and fast imaging modes. The slow mode functions on the 500Hz camera frame rate where displacements are mostly uncorrelated between successive frames. On the other hand, the fast imaging mode reaches a minimum frame rate of 25 kHz that is necessary for tracking propagation in order to perform transient elastography in soft biological tissues. The details of this set-up are provided here.

Recording quantitative phase images with digital holography

The set-up presented here is an off-axis interferometer based on a Mach-Zehnder, illustrated Fig. 4.1 [5]. The light source is a 671 nm continuous single mode LASER (CNI MSL-FN-671-S) with a total power output of 250 mW and up to 11 mW/cm² power density on the sample. A 1440 by 1440 pixels Complementary Metal Oxide Semiconductor (CMOS) camera (Adimec, Q-2HFW) records image-plane holograms at 500 frames/s. The field of view in the object plane is 11×11 mm².

The propagation of surface waves in soft tissue being in the order of one meter a second, a stroboscope-approach was used to reach an equivalent 25 kHz frame rate to image propagation over the field of view. This fast mode is necessary for classical transient elastography. Stroboscopic acquisition is achieved by synchronizing our mechanical source with our camera using a Field Programmable Gate Array (FPGA). Each acquisition consists in 500 frames, which are acquired 10 at a time. Between each 10-images acquisition, the time-delay between the generation of the mechanical impulse and the image acquisition is slightly extended. This time-difference corresponds to the period of the equivalent frame rate, here 25 kHz. The frames acquired chronologically then need to be reordered to reconstitute the propagation movie. The sensor of the camera acquires image-plane holograms (see section 2.5.2) containing phase and amplitude information. The holograms acquired are then processed to extract the relative phase and amplitude from each frame. A two 2D Fourier transform method is chosen here [6]. From the 1440 × 1440 raw hologram, the retrieved phase and amplitude images are 600 × 600 pixels. Compact and adaptable version for future preclinical trials. Most of the presented images are acquired on the set-up presented in Fig. 4.1. This set-up provides numerous degrees of freedom for the alignment with the two free-space beamsplitters. It is used for the control measurements, the proof of concept of optical noise-correlation elastography on a full-field system and for performing a majority of the experimental images presented here. The set-up is however limited to a fixed field of view by its design. Another set-up design was consequently developed around a single free-space beamsplitter and Lefevre polarization loops to achieve a more compact system with the possibility to easily adjust the field of view and numerical aperture by changing the objective. The second set-up design is presented in Fig. 4.1. Contrary to the previous design, this one piece set-up relies on two principal posts and can easily be fitted for a box to protect it from dust and to protect users from the laser radiations.

From quantitative phase images to displacement movies

The previously presented set-up provides full-field images of the phase and amplitude. The phase is especially interesting here as it can be directly linked to displacements within the sample. If the vertical displacement at the surface of the sample is noted ψ z , the phase ϕ retrieved by the holography set-up can be expressed as:

ϕ(x, y, z) = 4πnψ z λ LASER + ϕ noise (x, y, z) + ϕ speckle (x, y, z)[2π] . (4.1)
with n the refractive index of the medium at the interface with the sample, λ LASER the illumination wavelength, ϕ speckle the speckle noise and ϕ noise other 'standard' noises, notably from the light source and the camera. Speckle is distinguished from other noises due to its specific spatiotemporal characteristics. The noise is characterized experimentally in a following paragraph.

The phase difference dϕ between successive images is often used in interferometry to directly link the phase difference to the axial displacement ψ z , thus cancelling the noise considered constant between successive frames. This provides displacement images that can be used to compute shear wave speeds and perform elastography.

Using the phase-difference is quite successful in the slow mode where images are chronologically ordered 2ms apart. However, in the fast imaging mode, the "successive" images of the reconstructed 25kHz movie are actually taken more then 20ms apart. Therefore in Eq.4.1, the speckle pattern is not stable between successive frames. In consequence, the resulting phase difference in fast mode imaging are very noisy. The noise contributions are discussed in more details in the following paragraph.

Characterization of the experimental set-up

Noise estimations

In order to characterize the noise in the experimental images, two types of noise have to be considered here: the 'standard' noise of optical systems and the speckle noise from the coherent imaging.

Standard sources of noise.

The standard noise at the camera level consists of photon noise (from the photon detection process), electronic noise (from the sensor) and quantization noise (from analog to digital conversion). Values for these noise are characteristics of the camera. Another source of noise to consider is the laser fluctuations, especially the stability of its phase and power. This standard noise can be quantified using manufacturer data-sheets or experimentally by performing measurements.

In order to characterize the standard noise of the set-up, two experimental measurements have been made: the spatiotemporal characteristics of the dark noise and the shot noise fluctuations.

A typical acquisition was performed without the laser source with the usual integration time (200µs), camera frequency (500Hz), duration (1s) and ambient lighting. The spatial characteristics of this dark noise presented in Fig. 4.4 reflects the four different panels of the CMOS sensor. Each quarter of the captor has a Gaussian distribution of noise. As for the temporal characteristics of the dark noise, they are presented in Fig. 4.5. The average temporal offset recorded on the captor was a Gaussian at maximum 25% of the 11 bits captor. The dark noise is yet constant in time over a typical acquisition.

The second set of acquisitions consisted in using the reference beam to illuminate the captor with gradually a increasing intensity to characterize the shot noise. Shot-noise is often the dominant source of noise in sensors and can be modeled by a Poisson process. In order to ensure that the system was limited by shot noise, thirteen acquisitions were made with different laser intensities. The intensity and intensity (I) fluctuations were then calculated for each acquisition. Fig. 4.6 presents the intensity fluctuation in function of I/ √ I. The observed curve shows a linear tendency with a coefficient of determination of R 2 = 0, 97. This is characteristic of a Poisson process and indicates that shot-noise is the dominant source of noise of our imaging system.

Speckle.

The speckle noise has a very different structure than standard noise. As discussed previously in Chapter 2, speckle noise is a statistical process which depends on the light source coherence, on the numerical aperture and on the roughness of the considered sample. As described in Eq. 4.1, it is considered here that the recorded phase is the sum of the speckle pattern and the axial displacement of interest. The typical speckle grain for the demonstrated set-up is 5 pixels, corresponding to 40 µm. This was determined on a typical acquisition using the auto-correlation of the image [7]. The result of the auto-correlation is shown in Fig. 4.7. Fig. 4.8 a. and b. display respectively typical speckle on a phase image and its resulting artifacts on phase difference. Filtering was consequently used on the phase difference images to reduce these artifacts.

Median or median-based filters are relatively simple and used in ultrasound imaging [8] and radar images. More advanced methods have been shown to outperform the filtering in some cases [9]. These investigations on speckle denoising have yet mainly been performed on speckle images in amplitude. However, it has been shown in section 2.5.2 that the profile of speckle varies from the amplitude to the phase. In optics, various methods have been implemented in literature to filter out the speckle grain in coherent imaging [10] [11]. With this experimental set-up, we aim at filtering directly the phase difference images. Algorithm by [7].

Two simple numerical approaches were used with this set-up: a median filter and a threshold, linear interpolation and Gaussian filtering. A two-dimensional median filter was chosen with a size of 11×11 pixels, corresponding to twice the speckle's size. This approach is yet non-linear and time-consuming. Another approach was consequently implemented in order to reduce computation times. Considering the histogram of the phase difference images, it seemed that the signal could be easily extracted from the noise using a threshold. The value of the threshold was set depending on the histogram of the first image of the data-set. The filtered out pixels were filled-in with linear interpolation of the neighboring pixels and a Gaussian filter of size 9×9 is then applied to smooth the image.

These two methods are illustrated Fig. 4.8d. and e. and display similar results. The threshold approach is 7 times quicker than the median filtering. The threshold value is however more challenging to set as it can depend on the recorded displacement amplitudes. In practice, the method using median filtering is preferred in the cases where a lot of different data had to be treated without constraint of time. For quicker results on images with consistent displacement values, the threshold approach is used. ) with a better suited color-map to compare with the filtered images (e.) and (f.) using respectively a median filtering and a threshold and Gaussian filtering. The threshold used is illustrated in (c.) by blue dotted lines.

Sensitivity and resolution.

In order to determine experimentally the sensitivity and resolution of the experimental set-up, a resolution target was imaged as well as a mirror.

A mirror was put in the object plane of the set-up and a quater wave plate was placed in front of it in order to image its reflection. This was made in order to evaluate the sensitivity of the system. Indeed, this configuration allow the observation and measurements of phase shifts caused by perturbations in the imaging system and not due to axial displacements in the sample of interest. These perturbations in phase and amplitude measurements can notably be caused by fluctuations of the laser source or mechanical vibrations on the set-up, which we consider here the most limiting factor in terms of sensitivity. The acquisition with a mirror and a quarter wave plate in the object arm of the interferometer was performed on a typical acquisition of 500 frames at 500 Hz with a 200 µs integration time. Similarly to the camera noise estimations, the variations are observed over the reconstructed phase and amplitude. The average variation between two successive frames in laser intensity is 2.7 (in gray levels). This is equivalent to 0.1 % of the captor dynamic range. As for the phase variations, the average and maximum variations were respectively 0.34 and 0.51, corresponding to an average 20 nm phase error and 30 nm at maximum. As for the resolution, it was assessed experimentally using an United State Air Force (USAF) resolution target. An image of this target is displayed Fig. 4.9 at the resolution limit of the system. The smallest lines clearly visible are indeed the 20 pairs of line per mm, which corresponds to a 50 µm lateral resolution. ) is a zoom centered on the resolution limit of the presented optical system.

It should be yet noted that these values characterize the optical performances of the system. The resolution and sensitivity in terms of elastography images are a lot more difficult to evaluate as it is also dependent on the method and algorithm used for elastography. For shear-wave elastography, the resolution is in theory limited by the laser wavelength and the sensitivity by the SNR of the imaging system [12]. These considerations are discussed in following experimental results, notably in Chapter 5. § 4.

Implemented elastography approaches

In order to perform elastography using the presented set-up, a mechanical source is used to generate mechanical waves in order to probe the sample. The recorded response of the medium to this stimulation is then interpreted using the models described in the previous Chapter to extract mechanical properties.

For control measurements, a synchronized piezoelectric actuator (Thorlabs, AE0505D08F) was used at a 5V peak voltage to transmit mechanical impulses of hundreds of nanometers in amplitude to the sample. The impulses were transmitted to the sample through a metallic plate with a 10×1mm² section in contact with its surface. Harmonic and impulse signals were utilized for performing control measurements.

For noise-correlation elastography, at first three to four piezoelectric actuators were used around the sample or on its base to create a diffuse mechanical field with a controlled frequency range. Other types of uncontrolled sources were then tested such as fingers tapping, phone vibrations and a more powerful actuator (Cedrat Technologies, APA100M).

Harmonic approach

Harmonic elastography is performed here by sending a continuous harmonic signal in homogeneous sample. The average wavelength in the sample can then be retrieved in the slow imaging mode for different excitation frequencies. This approach has the advantages of only requiring one frame. Experimentally wavelength measurements are performed by binning the image in the propagation direction and extracting the dominant wavelength of the resulting signal. This method provides measurements which have been used to investigate the frequency-dependant properties of the homogeneous test samples.

Harmonic excitation was also often used experimental as a quick visual assessment of the proper functioning of mechanical source's and an overview of the sample's mechanical properties. Harmonic elastography was yet not fully investigated with the set-up as this work focuses on transient elastography. Digital holography and harmonic elastography have yet been successfully combined in previous studies [1].

Classical transient elastography

Classical transient elastography is implemented using the stroboscopic mode of the experimental set-up and a time-of-flight approach. An algorithm is implemented in Python to process propagation movies into speed-map images. This algorithm uses temporal correlations along the direction of propagation to determine the distance travelled by the wave front between the different frames.

Another approach for processing classical elastography propagation movies was implemented with the digital holography set-up based on spatio-temporal correlations. It has been shown to provide improved SNR on experimental images and less artifacts than time-of-flight [13]. It is introduced in Chapter 6.

Noise-correlation elastography

Noise-correlation elastography suggests to analyze a diffuse mechanical field in the sample in order to extract the medium's mechanical properties. This approach is performed with the experimental set-up using the slow acquisition mode. Different sources are tested experimentally for a range of mechanical sources, from controlled high-frequency piezoelectric actuators to fingers tapping randomly on the surface of the samples. A noise-correlation approach is then implemented on the captured images of the diffuse mechanical field in the sample to create quantitative stiffness maps. This approach is the focus of Chapter 5. More details on the method can be found there. § 4.

Simulation images

In order to test the elastography approaches and image processing algorithms, finite difference simulations were coded on Python. The propagation of shear-waves at the surface of a purely elastic non-absorbing medium was modeled based on the wave propagation equation 2.2 presented in Chapter 2 with an added source term:

∆u(r, t) - 1 c 2 ∂ 2 u(r, t) ∂t 2 = S(r, t) (4.2)
with u(r, t) the displacement and S(r, t) the source term.

This propagation equation was first written with Cartesian coordinates in two-dimensional space and time :

∂ 2 u(x, y, t) ∂x 2 + ∂ 2 u(x, y, t) ∂y 2 - 1 c(x, y) 2 ∂ 2 u(x, y, t) ∂t 2 = S n (x, y) (4.
3)

The derivatives are approximated using central differences for space and backward differences for time, transitioning to discrete coordinates. The following notation is introduced to simplify the expression of the space and time-dependency of the displacement : u n (x, y) = u(x, y, t = n) and S n (x, y) = S(x, y, t = n) with n an integer representing the temporal step. The final expression is: un(x+1,y)+un(x-1,y)-2un(x,y,t) dx 2

+ un(x,y+1)+un(x,y-1)-2un(x,y)

dy 2 -1 c(x,y) 2
un-2(x,y)+un(x,y)-2un-1(x,y)

dt 2 = S n (x, y) (4.
4) The calculus of the displacement field at a time t can consequently be calculated from the displacement field at the time t -1 and t -2. The boundaries are considered reflective. The speed-map of the considered medium is determined depending of the desired stiffness and the sources' position and temporal profile is fixed. Two main configurations are tested : a one-dimensional propagation and a diffuse field.

For the simulations of one-dimensional sources, the mesh is usually 4096×2048 pixels with a varying number of iterations depending on the source's signal's duration. The source is a vertical line positioned at a about a quarter of the simulated surface. Different functions are tested for the source: sinusoidal, hann, sine-cube and random functions are generated with different frequencies and durations. The central 500 × 500 pixels zone of interest is then selected to crop out the reflections at the interfaces.

The mesh for the simulation is typically set to 1024×1024 pixels over 8192 iterations for the diffuse field. Random uncorrelated noise sources with the same bandwidth are placed at a regular interval around the central 500 × 500 pixels final field of view. The resulting movies are saved starting at the 1500 th iteration to ensure that the created waves fully covered the region of interest. § 4.

Polymer Phantoms Manufacturing

As for the first tests of the experimental set-up, polymer phantoms were manufactured in order to obtain samples with controlled mechanical properties in order to test the experimental set-up. Details on their fabrication are provided here.

Agarose samples

Tissue-mimicking samples were made using agarose (A9539, Sigma-Aldrich, St. Louis, MO, USA) for mechanical properties and titanium dioxide (277370010, Acros Organics, Morris Plains, NJ, USA) for optical scattering properties. The mass concentration of agarose in the samples ranges from 0.5 to 2% depending on the desired stiffness. The reduced scattering coefficient was measured at 15 cm -1 at 650nm using a time-domain scattering measurement [14]. As for mechanical properties, the shear-wave speeds in such samples typically range from 1 to 10 m/s (corresponding to shear-moduli between 1 and 100 kPa), similarly to biological soft tissues [15]. The samples were 50 mL and solidified as approximately 65 × 65 × 15 mm 3 gels. All the homogeneous samples were prepared independently following the same protocol.

The agarose inclusion and heterogeneous samples were prepared by either cutting the gels or placing plastic rod vertically in the mixture while still liquid. Plastic rods with different cross-sections were used to create different shapes. After the sample cooled down to room temperature, the plastic rod could be removed and warm agarose gel with a different agarose mass concentration could be poured into the cavity, thus creating a softer or stiffer inclusion.

Structured samples were also created to mimic the presence of fibers in tissues and create anisotropic mechanical behavior. White cotton thread and rice noodles were used and incorporated at 1 to 2 mm below the samples surface.

PVA samples

Agarose samples are rapid to manufacture but not durable. The mechanical properties of agarose samples were also found to be very temperature-dependent. PVA samples were chosen as an alternative as they are easily found in elastography literature for comparison. PVA has another advantage over agarose as it can be used to create mechanically anisotropic samples without the need for a structuring [16].

PVA samples were made by using a mixture of 6% mass concentration of PVA powder, titanium dioxide and water. 150mL was poured into molds to form 75 × 150 × 15 mm 3 gels. The mixture underwent between 3 and 6 cycles of freezing at -25°C during at least one hour. The number of freezing cycles determines the obtained stiffness.

Because PVA solidifies when undergoing freezing cycle, it has been found that anisotropic PVA samples could be manufactured by applying tensile stress during the freezing cycle. Anisotropic samples were realized by stretching the gels during both the 4th and 5th freezing cycles.

Structured PVA samples were also realized by using custom molds to create cylindrical PVA fibers between 1 and 3mm in diameters. Three cycles of freezing were performed. When the fibers were formed in the mold, they were extracted to be aligned in a new mold where they were immersed in non-solid PVA. This fibrous sample then underwent four additional freezing cycles. § 4.6 Validation of the experimental set-up and the control measurements

Since common elastography methods couldn't be used for accurate control measurements, classical transient elastography and harmonic elastography were performed with the set-up as control measurements during this thesis. In this section, the classical transient elastography algorithm is validated and the manufactured polymer samples are characterized using both harmonic and classical transient elastography. Results on heterogeneous samples and ex-vivo biological tissues are also presented to show the limits of the classical transient elastography.

Validation of the classical transient elastography algorithm on simulation images

Validation for homogeneous medium

Homogeneous media with shear-wave speeds of 2, 4, 6, 8 and 10 m/s were simulated (250×250 pixels, 500 frames). Different functions were tested for the mechanical stimulation: sinusoidal, hann and sine-cube functions were generated with durations of 1, 2 and 4 milliseconds each. The shear-wave speed was calculated from the simulated movies using the time-of-flight algorithm. Each estimated shear-wave speed was calculated over the simulated area. The error was negligible as it was below 0.001 % and due to the spatiotemporal sampling used for the simulation. This first testing validated the timeof-flight algorithm.

Results on a medium with an inclusion

The time of flight approach was then tested on simulated propagation movies in heterogeneous media. Star-shaped inclusions (500×500 pixels, 1000 frames) were used to assess the effects of discontinuity and heterogeneity on speed maps. A shear-wave speed of 6 m/s was simulated for the background, and 3, 4.8, 5.4 and 5.7 m/s for the soft star-shaped inclusion. Conversely, stiffer star-shaped inclusions were also generated for a 3 m/s background. A propagative filter was used to improve the performances of the time-of-flight algorithm and reduce artifacts [17]. The results as displayed Fig. 4.10 for the case of a 3 m/s inclusion in a 6 m/s background, which is the case where artifacts were the strongest.

The inclusions can can distinguished in both cases. The directional filtering has significantly improved the reliability of the reconstructions. Fig. 4.10 still shows high variations in stiffness estimations especially in the cases where the inclusion is in the order of magnitude of the mechanical wavelength. Higher speeds are notably estimated around the inclusion and this effect seem to be increasing with the length of the mechanical excitation. These first experiments validated the time-of-flight algorithm coded for classical transient elastography acquisition but also highlights some characteristic artifacts in the speed map reconstruction.

Characterization of agarose homogeneous samples

Polymer samples are used throughout this thesis to test the performances of the developed elastography approach. In this section, homogeneous agarose samples are characterized using harmonic excitation in order to provide some preliminary speed measurements and assess viscoelastic behavior in the frequency range used for transient elastography approaches. Classical transient elastography is then performed experimentally to validate the approach and the experimental fast-imaging mode.

Harmonic testing

Agarose samples are prepared with agarose mass concentrations from 0.5% to 2.0%. To characterize the samples, sinusoidal excitation signals are used with the digital holography set-up. From images of the sample's surface, the spatial wavelength can be directly linked to the shear wave speed. As displayed in Fig. 4.11, the same 1 kHz mechanical source generates very different wavelengths in samples with different stiffness. This acquisition mode does not require synchronization and allows for quick mechanical testing of the samples.

The shear-wave speeds provided by these measurements are provided Fig. 4.12. The agarose gels used here can be considered as purely elastic in the frequency range of test-ing. Indeed, for each sample, the measured shear-wave speed shows low dependency on the excitation frequency. The resulting curve is monotonous, as expected, indicating that the measured shearwave speed increases the with agarose mass concentration of the sample. The shear-wave speeds provided by the transient elastography results are aligned with the harmonic measurements. This result validates the stroboscopic approach and the experimental set-up for recording propagation movies and perform transient elastography measurements. Increased variability of measurements with greater speeds also show the limitations from the temporal sampling. Indeed, wider wavelengths and fewer frames on the propagation movies decrease the time-of-flight measurements' accuracy.

The experimental set-up provides amplitude as well as phase propagation movies. For classical transient elastography, amplitude images are sometimes presented to illustrate some acquisitions as they are often more pleasing to the eye. This is justified as the recorded speed on the amplitude and phase movies are similar, as displayed in Fig.

Classical transient elastography on agarose heterogeneous samples

After first testings on homogeneous samples, samples were prepared with two different stiffness. An acquisition of a 1.0 % and 0.5 % sample is displayed Fig. Inclusions samples were then tested to evaluate the robustness of the experimental set-up to reflections and refraction. A star-shaped inclusion with similar contrast that in the simulation was created. The propagation recorded in this sample is shown Fig. 4.16.

These tests on controlled samples with heterogeneous mechanical properties showed the ability of the set-up to provide reliable shear-wave speed measurements using classical transient elastography. The acquisition in fast mode are yet too noisy to achieve stiffness images using the time-of-flight algorithm, even with a directional filter. The coded time-of-flight algorithm is not robust enough to retrieve images from experimental data in most cases due to the bad SNR in the fast acquisition mode. Indeed, because of the speckle decorrelation between successive frames of the acquired digital holography movies, the quality of speed maps is deteriorated. This algorithm's use was consequently limited to simulation images and control speed measurements were performed manually on the experimental propagation movies.

Classical transient elastography on biological tissues

Since the focus of this work is biomedical imaging, it was crucial to ensure that the experimental set-up and control methods were able to provide measurements in biological tissues. In this objective, porcine liver as well as murine nerves were tested ex-vivo. Because of the previously discussed limitations of the classical transient elastography algorithm, only stiffness estimations are presented here from spatiotemporal images of the acquisition.

Ex-vivo porcine liver

An ex-vivo porcine liver sample was also tested. The sample was cut into two parts: half of it was cooked and the other half was left raw. Both parts were then put back together for recording. The shear-wave speeds measured on each side in presented Fig. 4.17. The retrieved shear-waves speeds indicate that the cooked liver is stiffer than the raw, which is coherent with our expectations.

Ex-vivo murine nerves

Ex-vivo murine nerves were then tested. These nerves where tested under no stress and under tensile stress. Two clamps were used to maintain the sample in free space and apply tensile stress in the second case. Because the nerves presented were a lot stiffer than previous agarose samples, the stroboscopic approach was adapted to reach an equivalent 1.25MHz frame rate for these acquisitions. Fig. 4.19 present some of these results.

In these tests, the stroboscopic acquisitions presented a worst SNR compared to the previous porcine liver images. Even though the stroboscopic frequency was high and the images noisy, measurements of wave speeds could yet be performed. A higher wave speed was indeed recorded from the nerve under a tensile constraint. These quantitative result should however be interpreted carefully as the wavelength and the nerve diameter are close in size and phenomena of guided waves cannot be excluded. Fig. 4.19 shows indeed guided modes in the nerves observed under harmonic excitation with the same frequency as the pulse used for transient elastography in the previous figure. This shows that particular attention should be paid to the geometry of the sample at the scale of the mechanical wavelength as it can strongly affect speed-measurements.

Conclusion

The experimental set-up was presented in this section, including its characteristics and performances, as well as the image processing algorithms. This experimental set-up aimed at creating a benchmark for testing different elastography approaches on simulation as well as experimental images of manufactured polymer samples with controlled mechanical and optical properties. This chapter describes all these components. The results of harmonic and classical transient elastography control measurements are presented. The results presented for classical transient elastography highlighted the limitation of this approach when dealing with heterogeneous or discontinuous media or with the speckle decorrelation of present stroboscopic imaging.

The control experiments rely on the generation of a controlled and sometimes synchronized wavefront, which is a significant constraint for clinical application. In order to overcome this limitation, a noise-correlation approach was implemented. The introduction of the noise-correlation approach and the proof of concept of its adaption to full-field digital holography is presented in the following Chapter.

-Chapter 5 -

Noise-Correlation Elastography § 5.1 Summary

Transient elastography approaches for biological tissues were first introduced based on the tracking of a single mechanical impulse. For optics, the major limitation of this classical approach is the generation of a controlled mechanical wave and its tracking over smaller fields of view. Indeed, classical transient elastography requires tracking shearwaves propagating at 1-10 m/s in soft tissues. This entails the use of a fast camera or a phase-lock synchronization between the camera and the mechanical source as shown in the previous Chapter of this thesis. Another limitation of the classical transient elastography approach is that it is based on ballistic propagation of the mechanical waves. As displayed in the previous Chapter on simulation images and ex-vivo samples, the presence of scattering and interferences from tissues inhomogeneities or boundaries can hinder shear-wave tracking. Also, the presence of other mechanical waves from natural vibrations (heartbeats and other muscular activities) can limit the clinical translation of this method. New passive elastography methods have consequently been developed to work with complex elastic fields, whether it is generated externally [1], [2] or naturally present in-vivo [3], [4].

Approaches using noise-correlation from seismology [5] and time-reversal (TR) [6] have been recently adapted to elastography of soft biological tissues [7] [8] to overcome the main constraint of classical methods by exploiting diffuse shear wave fields. These approaches are based on a numerical refocusing of the shear-wave field by using noise-correlation algorithms. At each pixel, the refocused field gives access to the local shear-wave wavelength which is directly related to the local stiffness. This method was first developed for ultrasound systems but in 2016, Nguyen et al. [9] adapted noise-correlation to optical elastography. The proof of concept of the noise correlation approach on optical images was made on OCT images. Because of the nature of OCT images, correlation was calculated line by line. The full-field set-up presented in this thesis offers an adaptation of this noise-correlation method to full-field images where correlation can be done along two dimensions. The elastography method presented here is based on a mechanical noise field to image elasticity, which eliminates the need of an controlled impulse generator.

The noise-correlation approach is presented in this Chapter as well as the results combining it with the full-field digital holography set-up introduced in the previous Chapter.

123

The proof of concept is presented as well as results on inclusion samples and ex-vivo tissues [10]. § 5.

Principle of Noise-Correlation Elastography

The definitions and principle of correlation functions, time-reversal and noise-correlation are provided in this section.

Green's function

At first, an important notion to define in order to describe time-reversal and noisecorrelation is Green's function. In the case of linear time-invariant system, a Green's function of the system can be viewed as the impulse response of the system. A general definition of the Green function is presented here. If we consider a differential linear operator L, the Green's function G is defined as a solution of the equation 5.1, where δ(x) is the Dirac distribution in the x position

LG(x, x 0 ) = δ(x -x 0 ) .

(5.1)

If we apply this definition to the wave propagation equation, G is defined as [11]:

(∆ - 1 c 2 ∂ 2 ∂t 2 )G(r, r 0 , t) = -δ(r -r 0 )δ(t) . (5.2)
In free space, the classical form of the Green's function (corresponding to a diverging spherical wave from an impulse) is expressed as:

G(r, r 0 , t) = 1 4π|r -r 0 | δ(t - |r -r 0 | c ) . (5.3)
This function is a a solution of the wave equation for a point source placed in r 0 . With the knowledge of this function, the solution to the wave equation with an arbitrary source field S can be resolved. The solution of the general problem LΦ = S is the convolution of the Green's function with the source function. Considering a linear time-invariant system, the notations for G and the impulse response of the system h are used equivalently in this section.

Correlation in the case of a single source in one dimension

Let's first consider only one direction in space and a single source. This example aims at demonstrating the link between the correlation of the signal between two spatial positions and the impulse response of the medium. The temporal convolution operator ⊗ t is a commutative operator defined here as:

f ⊗ t g(t) = f (t -τ )g(τ )dτ .
(5.4)

When considering a mechanical source of shear-waves e at the spatial position x s , the axial displacement field ψ z within a lossless medium at a time t and position x can be expressed using the temporal convolution ⊗ t as ψ z (x, t) = e(x s , t) ⊗ t h(x, x s , t) .

(5.5) with h (x, x s ,t) the solution of the wave equation in the impulse regime for a receiver and a source localized respectively in x and x s . In this case, it corresponds to the Green's function of a medium between two points and represents the record we would obtain at one point if an impulsive source is applied at the other.

Let's consider the correlation, noted C, of the temporal signals from two separate points noted A(x A ) and B(x B ). It can be expressed as

C(x A , x B , t) = ψ z (x A , -t) ⊗ t ψ z (x B , t) .
(5.6)

By using the general expression of the axial displacement field in Eq. 5.5 and commutative properties, Eq.5.6 can be expressed:

C(x A , x B , t) = e(x s , -t) ⊗ t e(x s , t) ⊗ t h(x A , x s , -t) ⊗ t h(x B , x s , t) .
(5.7)

In a lossless one-dimensional medium, the space reciprocity principle states that the cross-correlation of the signal at two points r A and r B is equivalent to the response at one of these points as if the second point was a source [12]:

h(x B , x A , t) = h(x B , x s , t) ⊗ t h(x A , x s , -t) .
(5.8)

Eq.5.7 can consequently be rewritten as:

C(x A , x B , t) = S(t) ⊗ t h(x B , x A , t) .
(5.9)

with S(t) = e(x s , -t) ⊗ t e(x s , t) the auto-correlation of the source.

This equation is valid for all types of sources, including noise. If we consider an infinitely broadband source, the auto-correlation of the source function converges toward a Dirac distribution in Eq. 5.9 and the correlation function C gives direct access to the medium impulse response h.

The case presented here is the most basic case but the argumentation can also be developed in more complex cases (multiple sources, multiple spatial dimensions, anisotropic or viscoelastic behaviors) in order to link the recorded field to the medium's mechanical properties.

Time-reversal principle

Time-reversal is based on the recording and re-emission of an impulse through a medium. Receivers are placed to form a closed cavity around the source at a point r S in order to record the response (g i (t)) of a medium to an impulse source signal, as illustrated in Fig. 5.1. The recorded signal can then be reversed in time (g i (-t)) and re-emitted by the receivers. As a consequence, the wave equation states that the produced mechanical wave propagates through the same path but in opposite direction and refocuses at the source position. This concept is described here in order to link time-reversal with correlations and Green function retrieval. Notations of time-reversal are also introduced in this section.

In previous research, the time-reversed field has been generally written in function of the field ψ z for any spatial and temporal profiles of sources [4]:

ψ T R z (r, -t) = ψ z (r 0 , -t) ⊗ t ψ z (r, t) .
(5.10) with r 0 the emitter's position.

The definition is equivalent to the definition of the correlation function C of Eq. 5.6. The time-reversed field at a point r can consequently be expressed as:

C(r 0 , r, t) = ψ T R z (r, t) .
(5.11)

Noise-correlation principle

The closed cavity model described earlier doesn't exactly fit the elastography set-up and instead of the single impulse source described earlier, multiple sources, reverberation and diffusion are to be considered in noise-correlation. Also, the re-emission step of TR can be performed numerically based on cross-correlation computations, hence the name of the method.

Noise-correlation was first developed in seismology. Indeed, traditional methods in seismology are based on earthquake records and present short-comings similar to those encountered in transient elastography: earthquakes are localized and rare which limits the number of studies and the resolution in the areas where earthquakes are absent [13]. The use of ambient noise has allowed to overcome these limitations in seismology. The translation of such approaches to elastography is consequently based heavily on previous works in geophysics. However, elastography of soft tissues proceeds with very different scales and sometimes entails different hypotheses. For an imaging system, the receivers are the pixels of the camera, the correlation is computed here between the signal recorded at each pixel. Axial detection of the displacement ψ z (r, t) induced by shear-wave sources is considered. The principle is described here in different configurations from the most simple case to the one closest to the experimental considerations.

Elastic medium in one-dimension

We consider here only one spatial direction (x) and one temporal dimension (t). without affecting the generality of the foregoing, the source signals can be expressed as two sources: one propagating in the direction x and one in -x. The source's functions are respectively defined at the spatial positions x s and x ′ s . The corresponding wave-field ϕ z in a point x A can be noted:

ϕ z (x A , t) = h(x A , x s , t) ⊗ t e(x s , t) + h(x A , x ′ s , t) ⊗ t e(x ′ s , t) . (5.12)
Consequently, the correlation function between the points x A and x B can be expressed:

C(x A , x B , t) = [h(x A , x s , -t) ⊗ t e(x s , -t) + h(x A , x ′ s , -t) ⊗ t e(x ′ s , -t)] ⊗ t [h(x B , x s , t) ⊗ t e(x s , t) + h(x B , x ′ s , t) ⊗ t e(x ′ s , t)] (5.13)
If the sources are considered fully independent, as it is the case for a diffuse field, the correlation between the two sources can be considered null. Eq. 5.13 can be written:

C(x A , x B , t) = h(x A , x s , -t) ⊗ t e(x s , -t) ⊗ t h(x B , x s , t) ⊗ t e(x s , t) +h(x A , x ′ s , -t) ⊗ t e(x ′ s , -t) ⊗ t h(x B , x ′ s , t) ⊗ t e(x ′ s , t)] (5.14)
Since the medium is considered here as lossless and one-dimensional, we can use Eq. 5.8 in this case:

C(x A , x B , t) = h(x B , x A , t) ⊗ t S(t) + h(x B , x A , -t) ⊗ t S ′ (t) .
(5.15)

with S and S' the auto-correlation of the source signals e(x s , t)s and e(x ′ s , t).

Considering infinitely broadband sources, the auto-correlation of the source functions converges toward a Dirac distribution and the correlation function C directly yields the impulse response. This description is yet not adapted to the full-field displacement images provided by the optical set-up and more spatial dimensions should be considered.

Case of noise in a lossless scattering medium in two-dimensions

Let's consider here a two-dimension open lossless medium with scatterer randomly distributed. The correlation between the signal at two points A and B can be written as the sum of the contributions from each scatterer.

C(r A , r B , t) = C C ′ h(r A , r C , -t) ⊗ t e C (-t) ⊗ t h(r B , r ′ C , t) ⊗ t e C ′ (t) (5.16)
with e C the source contribution from C.

In order to obtain the link with the impulse response of interest (h(r B , r A , t)), a few assumptions have to be made here. First, we consider here that the medium is immobile so that propagation can be considered reciprocal. A continuous distribution of sources on a surface surrounding the medium is also assumed in order to form a perfect time-reversal medium. In this case, it can be written that [14]: If we consider the signal from the scatterer to be independent from each other, the terms in the correlation where C ̸ = C ′ will be null. This yields a single sum in the correlation expression:

C(r A , r B , t) = C h(r A , r C , -t) ⊗ t h(r B , r C , t) = h(r B , r A , t) + h(r B , r A , -t) (5.
C(r A , r B , t) = C h(r A , r C , t) ⊗ t h(r B , r C , -t) ⊗ t S C (t) (5.18)
with S C the auto-correlation of the source contribution from C.

If we consider the integration time of the correlation function superior to the correlation time of the source's noise, we have S C (t) = δ(t). Using Eq. 5.17, Eq. 5.18 can be written:

C(r A , r B , t) = h(r B , r A , t) + h(r B , r A , -t) .
(

These results from [14] show that the Green's function of the medium can be recovered from the cross-correlations in an open multiple scattering medium. This reasoning is valid under the assumptions of an immobile lossless isotropic medium and a continuous distribution of independent scatterers on a surface surrounding the medium. In the next section, the case of a viscoelastic medium is examined.

Viscoelastic medium with a homogeneous distribution of noise sources

We consider here the three dimensions of space and a homogeneous distribution of sources in the sample, all uncorrelated and acting at all times. The medium is considered here as viscoelastic. We considered the following governing equation for elastic wave propagation in a viscoelastic medium:

Lψ z - ∂ 2 ψ z ∂t 2 -2a ∂ψ z ∂t = -S (5.20)
with L an operator defined as L = c(r) 2 ∆, a an attenuation constant characteristic of the medium and S the source field.

The Green function of the medium is by definition the solution of this propagation equation with S = δ(t)δ(r -r 0 ). Its analytical expression for Eq. 5.20 is [15]:

G(r, r 0 , t) = Π(t)e -at sin(t) √ -L -a 2 √ -L -a 2 (r, r 0 ) (5.21)
with Π the Heaviside function and [|.|] the integral kernel of an operator, corresponding to the integral function equivalent to the operator. It can be defined for an operator M and function u as:

[M u](x) = [|M |] (x, y)u(y)dy (5.22)
If we consider now multiples isotropic sources at positions r i around in the region of interest ∆X, the resulting field in r A can be expressed as the following sum:

ψ z (r A , , t) = ∞ τ =0 r i ∈∆X h(r A , r i , τ ) ⊗ t s(r i , t -τ )dτ dr i (5.23)
with s the source function, representing in this case the white noise.

With the same consideration of independent sources and the previous equations, the correlation between the displacement field at r A and r B can be written [15]:

C(r A , r B , t) = σ 2 e a|t| 4a [|L ′ |](r A , r B ) . (5.24)
with L ′ an operator defined based on the operator L: its analytical expression is provided in [15].

We consider here white noise (containing all frequencies with random phases) acting everywhere in the medium at any random time in the form of homogeneously distributed uncorrelated sources. In this case, the correlation function can be directly linked to the Green function of the medium by derivation of Eq. 5.24 and using Eq. 5.21 to obtain the Green function [15]:

∂ ∂t C(r A , r B , t) = - σ 2 4a [G(r A , r B , t) -G(r A , r B , -t)] . (5.25)
with σ the variance of the diffuse field.

Consequently, if the attenuation coefficient is small enough, the time-derivative of the cross-correlation computed between the wave-field recorded at two points is proportional to the Green's function of the medium.

General considerations

A few cases have been discussed where, under different conditions on the medium, the elastodynamic Green's function can be retrieved from the cross correlation of two recordings of field at different locations. These approaches can be generalized in order to limit assumptions on the medium parameters [16].

The assumption required to obtain the link between the correlation function and the Green's function presented in Eq. 5.25 are difficult to obtain experimentally. However, the simple assumption of equipartition of the sources at boundaries if sufficient to have a direct link between cross-correlation and Green function.

More generally, the following equation is used as the base for noise-correlation elastography:

∂ ∂t C(r A , r B , t) ∝ G(r A , r B , t) -G(r A , r B , -t) . (5.26)
The proportionality coefficient depends here on the initial assumptions on the medium. By considering that the bandwidth of the considered displacement field is finite, we have [17]. Consequently, most noise-correlation approaches are based on the following general equation:

∂ ∂t C ∝ C
C(r A , r B , t) ∝ G(r A , r B , t) -G(r A , r B , -t) .
(5.27) § 5.

Algorithms for noise-correlation elastography

Different methods have been used in elastography in order to extract a medium's properties from the diffuse field based on Eq. 5.27. Four approaches have been identified [18] based on either phase velocity [19], focal spot size [20], amplitude or derivation ratio. The focal spot size and derivative ratio approaches are introduced here as they are both utilized experimentally and well-suited to low frame-rate imaging systems, allowing shear-wave speed computations in the absence of sufficient temporal sampling.

Focal spot size method

This method consists in refocusing the diffuse field at t=0 using cross-correlations. This yield a function with a specific spatial width, which can be linked to the shear wave speed. For isotropic lossless media with an equipartition of sources, Eq. 5.28 describes this link [18]:

r L (θ) = c s 2πf 10(1 -L) 2 -cos 2 (θ) .
(5.28) with θ the polar angle, r L the distance from the focusing point to the contour at which the correlation field takes a fraction L of its maximum value and c s the shear-wave speed and f the effective frequency of the field.

For the Full Width at Half Maximum (FWHM) at the maximum angle, we have:

F W HM = √ 5 c s 2πf . (5.29)
The method is illustrated on simulation images in Fig. 5.2.

Derivation ratio method

In order to exploit the link between the correlation function and the Green function of the medium, a method based on derivatives of the time-reversed field was principally implemented during this thesis. This method was first introduced in ultrasound elastography [4] and then applied to OCT images [9]. In this thesis, it was adapted to full-field digital holography images.

Two cases can be considered here: the case of a low frame-rate scanner where the different acquisitions of the diffuse field are uncorrelated and the case of a high speed scanner where consecutive wave-field images are correlated.

High frame rate scanning

In this first approach, the spatial and temporal derivatives of the diffuse field ψ z are used to extract the shear-wave speed. They can be approximated using plane wave decomposition as [4]:

ξ z (r, t) = ∂ψz(r,t) ∂r ≈ ikψ z (r, t) V z (r, t) = ∂ψz(r,t) ∂t ≈ -iωψ z (r, t) . (5.30)
At r 0 and t=0, their time-reversed equivalent are expressed :

ξ TR z (r 0 , r 0 , t = 0) ≈ -k 2 ψ TR z (r 0 , r 0 , t = 0) V TR z (r 0 , r 0 , t = 0) ≈ -ω 2 ψ TR z (r 0 , r 0 , t = 0) . (5.31)
The local shear-wave speed can be expressed using the ratio between both the spatial and temporal derivatives [4]:

c s (r 0 ) = ω k ≈ V TR z (r 0 , r 0 , t = 0) ξ TR z (r 0 , r 0 , t = 0)
.

(5.32)

Low frame rate scanning

In the context of a low-frame rate set-up, the temporal incoherence of the captured noise field prevents the calculus of temporal derivatives. Only the time-reversed strain field ξ TR z is used here. Considering its expression in Eq. 5.31 and the wave number k = 2π/λ, at a given time t 0 and location r 0 , the wavelength depends on the time-reversed displacement field ψ TR z and its derivative, the time-reversed strain field ξz TR :

λ(r 0 ) ≈ 2π - ψ TR z (r 0 , t 0 ) ξ TR z (r 0 , t 0 ) . (5.33)
This equation is equivalent to calculating the derivative of the correlation function at t=0. In Fig. 5.2 c., the corresponding calculus is the curvature at the peak of the curve (t=0). Compared to the previous approach, the frequency cannot be calculated here. The link between the wavelength in Eq. 5.33 and the shear-wave speed consequently required the knowledge of the central noise frequency.

Implemented algorithms for full-field images

Three different algorithms were implemented experimentally for noise correlation using MATLAB. Each algorithm corresponds to one of the processing methods presented previously: focal spot size, derivation approach for low and high frame-rates. The input of these algorithms is the full-field phase movie acquired with digital holography. This input can be expressed as a matrix Ψ z (x i , y i , t i ) with x i and 'y i two integers representing the spatial dimensions and t i the frame number of the phase movie. For the presented setup, we have 1 ≤ x i ≤ 600, 1 ≤ y i ≤ 600 and 1 ≤ x i ≤ 500.

• For the focal spot width approach, the algorithm simply consisted in tracing the cross-correlation function C in two-dimensions for t=0. The discretized expression implemented for computing the focal spot (PSF) at a point (x 0 , y 0 ) was:

P SF x 0 ,y 0 (x, y) = |Ψ z (x 0 , y 0 , t) × Ψ z (x, y, t) T | 2 .
(5.34)

with T the transpose of the matrix or vector.

The focal spot reconstructed from this equation is centered on the position (x 0 , y 0 ). Consequently, in order to avoid cropped images and fully reconstruct the focal spot, the coordinates of x 0 and y 0 were restricted to the center of the image (150 < x 0 < 450 and 150 < y 0 < 450).

The computation of the two-dimensional correlation takes 226 seconds for one position. For this reason as well as other reasons described previously, the focal spot is only performed here for verification purpose. This approach is mostly performed experimentally in homogeneous media. In this case, the resulting focal spot image is averaged over a hundred points.

• For the high frame rate deviation approach, Eq. 5.32 was implemented. This algorithm is mostly used combined with the stroboscopic approach (see in the following Chapter) as the frame rate of the camera doesn't allow sufficient temporal sampling to compute the temporal derivative of the correlation function.

The two main quantities in this Eq. 5.32 (the spatial and temporal derivatives) are computed from the phase movie utilizing the MATLAB function imgrad to process the central gradient with GPUarray. This yields two images: a wavelength map and an angular frequency map. The square root of the ratio of these images can then be computed to obtain the shear-wave speed map.

• For the low frame rate deviation approach, Eq. 5.33 was implemented. The algorithm is similar to the high-frame rate approach except that only the spatial derivative is computed to obtain the wavelength map. The GPU processing allowed for typical processing times of 2.3 seconds. Since the the high-frame-rate approach cannot be used with the slow-mode of the experimental set-up, this algorithm was used in this Chapter for the proof of concept of noise correlation. § 5.4 Results

The noise-correlation algorithm was tested on simulation images and with the experimental digital holography setup presented in the previous Chapter. The proof of concept on agarose samples and results on heterogeneous samples as well as ex-vivo biological tissues are also presented. Wavelength measurements and maps are presented as their are the direct results of the noise correlation algorithm.

Validation on homogeneous agarose samples

In order to assess the reliability and the quantitativity of the noise correlation elastography, homogeneous samples have been prepared with gradual agarose concentration (1.00, 1.25, 1.50, 1.75 and 2.00 % mass concentration). These agarose samples are tested consecutively with noise-correlation and with classical transient elastography.

For the proof of concept of the method, high frequency excitation was chosen in order to be able to fully reconstruct the PSFs for using the focal spot width approach. Piezo-electric actuators are placed randomly around a region of interest and a function generator delivered chirps between 800 and 1600 Hz with 0-3V voltage and with random durations (133-501 ms). The generated diffuse field at the surface at the sample was then imaged during one second (500 frames). Five different zones were imaged for each sample. From these acquisitions, the mean wavelength was then calculated from the phase images over the field of view and averaged. The measurements are then compared to the control method that provided the shear wave speed (see Chapter 4 for more details on the control method).

The results are presented in Figure 5.3 show the mean wavelength obtained with noise correlation depending on the wave velocity given by the control experiment. The corresponding error-bars are calculated from the different zone of the sample acquired. The red line represents the theoretical wavelength in the medium for a mean actuator frequency of 1.1 kHz. The wavelength assessments are all within an 8 percent dispersion of the theoretical curve. This deviation could be explained by an imperfect coupling between the piezoelectric actuators, the plastic box holding the sample and the agarose gel, possibly modifying the frequency content of the excitation. Independently, the full impulse response of each sample was traced by correlating the temporal signal of one pixel to the temporal signals of the neighboring pixels [4], [9], [19]. This impulse response calculation was made randomly on different points of the homogeneous samples on each acquisition to average it over its surface. The impulse response of each sample further confirm the results given by noise correlation as the PSF size visually increases with stiffness.

Inclusion agarose samples

The contrast and resolution of the full-field noise correlation method was then validated on samples with heterogeneous mechanical properties. Results on a 1% agarose sample with a 2% agarose inclusion are presented here. The parameters for the mechanical sources are first chosen identical to the previous experiment: random sources were positioned around the inclusion and a function generator delivered uncorrelated chirps (800 to 1600 Hz, 0-3V, 133 to 501 ms). The generated diffuse field at the surface at the sample was then imaged during one second (500 frames).

A description of the sample and mechanical excitation can be found in Fig. 5.4a. A white light image of the sample's surface shows its homogeneous optical scattering properties Fig. 5.4b. and four of the successive 499 phase difference images obtained during one acquisition are shown in Fig. 5.4c. It can be noted that the low temporal sampling (500 Hz camera frame rate) doesn't allow for any correlation of the field between successive frames. The reconstructed wavelength map in Fig. 5.4d. shows high contrast and super-resolution. The reconstructed PSF for the two zones of the sample show that the FWHM are comparable to the ones obtained in Figure 5.3.

These first experiments validated the experimental system with controlled sources of noise. The obtained shear-wave speed were highly correlated to the shear wave speeds obtained with the control measurements. Images of heterogeneous samples showed contrast and sub-wavelength resolution. The reconstructed PSF also showed comparable results with the images obtained with the derivation approach. The wavelength map image is retrieved using the noise correlation method. Averaged images of the correlation functions for a zone in the softer background (e.) and for a zone in the inclusion (f.).

Ex-vivo biological tissue

As the perspective for this set-up is biomedical imaging, preliminary testing on biological tissue was made. Figure 5.5 presents the wavelength map on a pork liver sample and on chicken skin.

For the liver, the tissue was cut into two and half of it was cooked. This stiffened the tissue, creating mechanical contrast. The two pieces of liver were then placed next to each other and imaged. For the chicken, a nylon rod of a few millimeters in diameter was inserted in the thigh of a chicken leg to create the contrast. The image were obtained with mechanical noise generated by actuators around 1.5 kHz.

The liver has a low scattering compared to other organs like the brain or the skin [21], which can explain lower contrast compared with the previous images of agarose samples. The image in Fig. 5.5a. shows yet good contrast between the two sides : the mean wavelength value on the raw side is 4.3 mm and on the cooked side 9.3 mm, indicating as expected that the cooked liver is stiffer than the raw liver. Considering no dispersion, these values would correspond to a shear modulus of 41 and 193 kPa, respectively.

In Fig. 5.5b. the nylon inclusion was highlighted with dotted lines. The contrast between the nylon rod and the chicken skin is clearly visible on the noise-correlation image.

Figure 5.5: Wavelength maps obtained on a pig liver and on chicken skin using noisecorrelation and digital holography. (a.) presents results on a porcine liver where the upper half was left raw as the other part was cooked. The average full cross-correlation is traced for each red square on the image for a 300 x 300 pixel window : (b.) and (c.) correspond respectively to the raw and cooked liver. (d.) images a circular nylon rod inserted in a chicken thigh. The nylon rod is highlighted with a red dotted line.

Unknown noise on agarose inclusion samples

The high-frequency controlled noise sources used for validation are very different from natural noise sources in-vivo in terms of their frequency range. The experimental set-up was consequently tested with unknown low-frequency noise-sources.

In order to access the quantitative stiffness from the wavelength maps, the mean frequency of the generated noise is necessary to link the wavelength map to the shearmodulus. When noise is generated by piezoelectric actuators, the mean frequency of excitation is a priori known, although deviations can occur due to the coupling between the source of noise and the sample. For this technique to be translated to physiological noise or any unknown broad-band noise, the mean noise frequency must be acquired precisely to quantify the stiffness. The camera can only capture the noise-spectrum in a limited range determined by its frame rate and acquisition's length, here (1-250Hz). In the following experiment, the generated noise frequencies can be considered in this range.

Fingers were used to tap on the sample's surface to generate a noise field with an unknown frequency content, expected with low-frequencies. An inclusion sample was chosen and a one-second acquisition was acquired. The sample and mechanical excitation is described in Fig. 5.6a. Because of the low frequency content, the PSFs were too wide to be imaged with the experimental system but the derivation approach for processing the noise-images was well-adapted and the reconstructed shear-wave wavelength map shows good contrast and the inclusion can be delineated. The low-frequency content allowed the computing of the noise spectrum from the acquisition. Its average over the field of view is presented in Fig. 5.6c. The central frequency is ten times smaller than for the actuators used previously. As a consequence, the retrieved wavelengths are ten times bigger.

Conclusion

This Chapter introduced the noise-correlation approach used to perform 'passive elastography' with digital holography. The experimental results show that the designed imaging system is able to provide full-field one-shot quantitative phase images, which is well-adapted to perform two-dimensional noise-correlation and obtain quantitative stiffness images. High-frequency sources were first used in order to present the focalized full correlation function to visually assess the focal width as well as the derivation approach for the proof of concept. The adaptability of the method to unknown low-frequency noise and ex-vivo biological tissues was then demonstrated.

Compared to classical transient elastography, this approach is compatible with a lowframe rate and lower frequency physiological noise (< 200 Hz) as the resolution of this elastography method has been shown to not be limited to the mechanical wavelength [22]. Even though the signal to noise ratio decreases at lower frequencies, noise down to 110 Hz has been used with satisfying contrast on polymer samples. This demonstrates the possibility for this technique to be totally passive, quantitative and contactless using broadband noise source, from physiological noise to unsynchronized external vibrators.

The setup is using correlation on 2D data which intrinsically gives access to the directional displacement of the sample. As confirm by the circular PSF retrieved in agarose, its mechanical properties are isotropic but in-plane anisotropy could be measured with this method on adapted samples by correlating the noise and obtaining the directiondependent Green's functions. This is discussed in a following Chapter.

One of the main assumption in order to perform noise-correlation is the equipartition of sources around the field of view. This is however not always the case experimentally. For example, when only one actuator is used in the sample to generate noise, a strong direction-dependence is observed in the stiffness estimation and influences the PSF's central symmetry. Considering isotropic media, the two-dimensional PSF provided by the experimental setup can inform on the energy repartition of the sources. The following Chapter explores how this information could be exploited.

-Chapter 6 -Source and Medium Anisotropy in Noise-Correlation Elastography § 6.1 Summary

The elastography approaches presented previously have described two main methods to perform transient elastography: one using one-dimensional sources of mechanical excitation for control measurements (Chapter 4) and one based on noise fields (Chapter 5). Both these methods were considered only for isotropic media. This section aims at bridging the gap between these two approaches by looking at the influence of source and medium anisotropy in speed-map reconstructions. Indeed, unlike ultrasound or OCT systems previously used for elastography, the fullfield approach of digital holography enables two-dimensional correlations. This provides directional information on the response of the sample. From this data, information on the source or the medium's anisotropy can be extracted.

The first part of this Chapter focuses on the use of a noise-correlation-based algorithm for non-noise sources such as spatially coherent shear wave sources used in classical transient elastography. By exploiting the two-dimensional correlation, the medium's response to the impulse can be processed to extract the local propagation direction and shear-wave speed. The results on full-field images are shown to be more robust to reflection and refraction phenomena than a time-of-flight approach [1].

The following section presents how to utilize the full-field system to extract the in-plane anisotropic properties of our samples using the noise-correlation elastography approach presented in the previous Chapter 5. Simulation as well as preliminary experimental results on anisotropic samples are presented here. § 6.2 Noise-Correlation-inspired Approach for Spatially Coherent Shear-Wave Sources

Elastography based on spatially coherent shear-wave sources such as the impulses used in classical transient elastography is often based on the tracking of wavefronts. Common approaches used to extract stiffness from propagation movies are detailed here as well as a new approach inspired by diffuse field elastography.

Existing algorithms and models for image processing for spatially coherent shear-wave sources

Algorithms found in literature can be grouped into two categories for classical transient elastography: inverse problem and time-of-flight approaches. These two approaches of classic transient elastography are described in this section. A third method, introduced during this thesis, is presented in the next section.

Inverse problem algorithms

Inversion algorithms were first used to retrieve the shear-wave speed from acquired spatiotemporal data. The first method performed a straightforward inversion of the Helmholtz equation 3.27 by calculating the stiffness as the ratio between the second temporal derivative and the Laplacian operator [2][3], [4] under the assumption of a linear, isotropic, incompressible, and piece-wise homogeneous material. This approach is however very sensitive to noise due to the temporal and spatial second-order derivatives.

It is also sensitive to changes in the propagation direction.

Global inversion approach based an the adjoint method inspired from seismology was shown to improve the image quality and sensitivity on simulation images [5]. However, time-of-flight algorithms are the main image processing approaches found in recent studies.

Time-of-flight algorithms

To retrieve the local speed of the waves from a propagation movie, time-of-flight is the common approach in transient elastography. The principle is quite simple: by looking at the spatial displacement of the propagating wave front between different frames, the shear-wave speed can be directly measured from the spatiotemporal data or computed using correlation-based methods. Indeed, by cross-correlating two adjacent lines in time to find the local propagation delay, the shear wave speed can be easily computed.

The time-of-flight approach has been successfully implemented in a commercial scanner as well as in most recent transient OCE papers to process propagation movies [6]- [17]. The major limitation of this method is the need for a known direction of shear-wave propagation and a clear wave-front, which reduces the possible mechanical sources, usually one-dimensional narrow-band pulse sources, and creates artifacts in heterogeneous media caused by the change of direction of shear-waves on the interfaces. Under the assumption of a low contrast, a directional filter can be applied to limit the artefacts from reflections [18] as shown in Chapter 4.

The resolution and accuracy of time-of-flight approaches depends on the kernel size of the algorithm (the distance between two considered points for the time-delay computation). The larger the kernel, the smaller the relative error on the speed estimation but the lower the image resolution. While time-of-flight is generally used with narrow temporal pulses, it has lately been extended to broadband pulses by using inverse filtering. This approach has allowed to significantly increase the signal-to-noise (SNR) ratio [12]. There is however a trade-off to be considered between increasing the shear wave bandwidth and the SNR as higher frequency shear-waves undergo higher attenuation in biological tissues.

Noise-correlation inspired (NCi) approach

Based on previous studies in elastography using shear-wave noise-fields for ultrasound imaging [19] and optics [20]- [22], a new approach was introduced to retrieve the shearwave speed from propagation movies. This approach is consequently referred to as Noise-Correlation-inspired (NCi).

The principle of noise-correlation elastography is described in the previous Chapter. Noise-correlation uses spatiotemporal correlations to extract the medium's stiffness properties from acquisitions of a diffuse mechanical field. We present in this section how the algorithm developed for diffuse field elastography [22] can be used with spatially coherent sources, for example in the shape of pulses or chirps used in classic transient elastography.

It has been demonstrated for seismic applications that correlation-based approaches can be used for both one-sided transient sources [23] and two-dimensional noise-fields [24]. The method used is explained in detail in the rest of this section.

In a lossless medium, the general expression of the temporal correlation between two points in space (r A and r B ) can be written using Eq. 5.9. Its expression is reiterated here:

C(r A , r B , t) = S(t) ⊗ t h(r A , r B , t) . (6.1)
with S(t) the auto-correlation of the source and h the impulse response of the medium.

Contrary to noise-correlation elastography, the source cannot be considered infinitely broadband here. When considering spatially coherent sources, the contribution of the auto-correlation of the source has to be taken into account in the expression of crosscorrelation. This means that in order to use temporal cross-correlation for transient elastography, knowledge of the source temporal shape is necessary to retrieve the impulse response from the correlation C.

If we assume an isotropic lossless medium with a constant elastic-wave velocity c, the impulse function h can be expressed as a Dirac distribution δ [25]:

h(r 0 , r A , t) = δ(t - |r 0 -r A | c ) . (6.2)
By using Eq. 6.2, Eq.6.1 can be rewritten as:

C(r 0 , r A , t) = S(t - |r 0 -r A | c ) . (6.3)
Eq. 6.3 shows how temporal and spatial variations of the correlation function can be directly linked by the mechanical wave speed c. Consequently, the functions C and S are only linked by a change of scale. The measure of this change of scale gives quantitative access to the speed of the perturbation, which is here the shear-wave speed c shear-wave .

The NCi approach suggests to take advantage of the simple relation between the cross-correlation C and the auto-correlation S to show that the algorithm developed for passive elastography can be applied to transient elastography. Two approaches are presented in the rest of this section, inspired by the focal spot size and the derivation methods used in noise-correlation elastography and introduced in Chapter 5.

The source auto-correlation S(t) can also be expressed using the correlation function at the position r = r 0 . The notation C(r 0 , r 0 , t) is preferred for the auto-correlation S(t) from this point.

The FWHM method

A simple way to quantify the scale between correlation function variations in time and space is to use the FWHM of each function. This is equivalent to the focal spot size approach presented in section 5.3.1. By tracing both the cross-correlation function at t = 0 and the auto-correlation at r = r 0 , the ratio of their FWHM can be estimated. With the developed experimental setup, we have access to the 2-dimensional cross-correlation. The FWHM must consequently be taken along the propagation direction. This first approach is robust but requires a full spatial correlation reconstruction for each spatial position, which is time-consuming and also deteriorates the spatial resolution. In our case, the FWHM method is particularly limiting as the size of the field of view with optical imaging is often too small to access full spatial correlation functions. A second method is preferred, making use of an algorithm for passive elastography based on the derivation approach [19] developed in our previous work [22].

The NCi method

The method used to perform shear-wave elastography using noise-correlation with highframe rates is detailed in Chapter 5. This approach provides the following formula for shear-wave speed estimation:

c s (r 0 ) = ω k ≈ V TR z (r 0 , r 0 , t = 0) ξ TR z (r 0 , r 0 , t = 0) . (6.5)
with ξ TR z and V TR z respectively the spatial and temporal second derivatives of the displacement field. These quantities are also called time-reversed strain field and timereversed particle velocity field [19].

The time-reversed strain field and the time-reversed particle velocity field are respectively the spatial and temporal second derivatives of the time-reversed displacement field. These fields can also be expressed in function of the correlation function:

       V TR z = ∂ 2 C ∂t 2 , ξ TR z = ∂ 2 C ∂r 2
.

(6.6)

The relation 6.5 can consequently be interpreted for transient elastography as the ratio between the temporal and the spatial curvature of the correlation function at t = 0 and r = r 0 . From this ratio, the change of scale between the temporal and the spatial correlation can be calculated, which corresponds to the shear-wave speed. Using Eq. 6.3 can be applied, V TR z and ξ TR z are linked directly by the shear-wave speed:

V TR z = ∂ 2 C ∂t 2 = ∂ 2 C ∂r 2
.

∂r ∂t

2 = ξ TR z .c 2 s (6.7)
The equation 5.32 used in diffuse-field elastography for high-frame rates consequently remains true for one-dimensional signals and the algorithm developed for noise-correlation elastography [22] can be used to process transient elastography images. The NCi approach is based on correlations and can consequently be performed on one-dimensional mechanical perturbations with any temporal profile.

In contrast to the FWHM method that requires the full cross-correlation and autocorrelation curves, the derivation approach of the NCi method only requires a few points to compute the spatial derivative of the auto-correlation function at t = 0. The second approach is consequently be preferred for processing stiffness images.

Adaptation to full-field images

When adapting this NCi approach to a full-field imaging system with one dimensional mechanical sources, the propagation direction is intrinsically calculated by the NCi algorithm as it corresponds to the direction for which the gradient of the two-dimensional spatial correlation (or the width of the focal spot) is minimal. This was implemented on the derivation approach by simply using the trace of the Hessian matrix corresponding to at the considered pixel.

The two main approaches are tested experimentally : the time-of-flight processing and the NCi approach using both the FWHM and derivation approaches for validation measurements. However only the derivation approach is used on images. These imageprocessing methods are summarized in Fig. 6.1.

Results

Simulation images

At first, the NCi approach was tested on finite difference simulations described in Chapter 4. Propagation in homogeneous media with shear-wave speeds of 2, 4, 6, 8 and 10 m/s was utilized with different mechanical excitation functions: sinusoidal, hann and sine-cube functions were generated with different durations (1, 2 and 4 milliseconds). The shear-wave speed was calculated from the simulated movies using the NCi approach with the derivation algorithm and FWHM. The estimated shear-wave speed was estimated over the simulated area : the retrieved values were correct and homogeneous with negligible ( < 0.6 %) difference between the measured and the simulated shear-wave speeds with both methods.

Non-differentiable source functions such as rectangular or triangular pulses were also tested The results from the NCi approach were not conclusive in this case as the method is based on the calculus of derivatives. In this case, the FWHM approach was successfully used to retrieve the shear-wave speed accurately on the central pixels.

Simulations of shear-wave propagation through inclusions were tested to assess the effects of heterogeneous media on speed maps. Two configurations were tested : soft and stiff inclusions. A shear-wave speed of 6 m/s was simulated for the background, and of 3, 4.8, 5. Since the artifacts were more visible with greater stiffness contrast, shear-wave speedmaps are shown Fig. 6.3 and Fig. 6.4 for the maximum simulated contrast. The NCi algorithm is tested with 4 and 1.33 millisecond-long half sine excitations functions. In order to demonstrate the adaptability of the NCi algorithm to any one-dimensional excitation signal, the NCi approach was also tested with a 4 millisecond-long chirp excitation function from 500 to 700 Hz and a 4 millisecond noise signal from 400 to 800 Hz. Time-of-flight is also used implemented on the images with the half-sine excitation as the reference (approach described in Chapter 4).

The images obtained show that the inclusion can be resolved with both methods using half-sine mechanical source signals. The obtained images display a resolution much smaller than the mechanical wavelength, which was expected as super-resolution has already been demonstrated for transient [3] and passive elastography [25].

Even with a directional filter, the images from the time-of-flight algorithm display specific artifacts. These artifacts probably stem from reflections, refraction and scattering occurring at the interface between the background and the inclusion. On the NCi image, artifacts are greatly reduced. A possible explanation is that the predominant direction of propagation is inherently chosen locally by the NCi algorithm. Changes in the propagation direction consequently do not impact the speed estimations as much. The images from the NCi algorithm all show the same resolution regardless of the source signal for both stiff and soft inclusions.

These findings validate the NCi method and showed its potential for applications in experimental images. The NCi algorithm was indeed able to perform quantitative full-field stiffness measurements from heterogeneous simulation images using different excitation signals with minor artifacts compared to time-of-flight. The simulation images demonstrate that the NCi approach is a promising alternative to classic time-of-flight algorithms for heterogeneous media. ) and (d.) for the 4ms-long pulse, (c.) and (e.) for the 1.33 ms-long pulse), chirp function (f.) and noise (g.).

Experimental images on agarose samples

Two types of excitation signals were used experimentally in this section. The first excitation signal consists in a repeated 1 ms-long half-sine function. In this case, the propagation of two to five pulses was recorded on each acquisition. The second type of excitation signal is a 2 ms-long chirp function, with frequencies ranging from 500 Hz to 2 kHz. The propagation of a single pulse is imaged on each acquisition. These acquisitions are respectively referred to as pulse and chirp acquisitions to simplify notations.

To validate the NCi approach experimentally, samples with gradually varying agarose mass concentrations (0.5, 0.75, 1.0, 1.25, 1.5, 1.75 and 2.0 %) were imaged using the digital holography setup in the fast acquisition mode (25 kHz). Pulse acquisitions and chirp acquisitions were obtained on these samples. Rayleigh wave speeds were measured consecutively with a time-of-flight method as well as with the FWHM method and the derivation based algorithm. Six acquisitions are acquired for each sample. The Rayleigh wave speed values obtained for the samples of 0.5, 0.75, 1.0, 1.25, 1.5, 1.75 and 2.0% agarose mass concentration are respectively 2.3 ± 0.1, 3.7 ± 0.3, 4.7 ± 0.2, 5.5 ± 0.2, 6.9 ± 0.2, 7.1 ± 0.2 and 8.9 ± 0.5 m/s with the reference time-of-flight measurements.

Results for the FWHM measures are provided in Fig. 6.5 for pulse acquisitions. As for the derivation method, the corresponding curves on pulses and chirps are in Fig. 6.6a. and b. respectively. Even if errors can be observed on the experimental data depending on the method, all the curves display a linear tendency, which means that the measurements are similar to those of time-of-flight with comparable tendencies. The correlation between time-of- flight and the FWHM measurements is indicated by a R2 coefficient of 0.93 while results on both excitation signals for the NCi approach have a R2 coefficient of 0.98. These results validate the NCi approach for quantitative experimental measurements of the phase velocity. The adaptability of the NCi algorithm to any one-dimensional excitation signal is also demonstrated by the experimental results from broad and narrow-band sources.

To illustrate the feasibility of the method on inhomogeneous samples, a 2.0 % agarose sample with a 1.0 % agarose inclusion was tested with a pulse acquisition. Fig. 6.7a. illustrates this experiment. Images of the pulse propagation are shown in c. as well as the reconstructed stiffness image in d.

Although the sample has homogeneous optical properties, the soft inclusion can be distinguished from the propagation movie by using the NCi algorithm developed for diffuse fields. A quantitative speed map is retrieved with contrast corresponding to the values obtained previously for homogeneous samples. Because only one wave front in one direction is used here, an overestimation of the shear-wave speed value can be noticed in the upper corner of the image.

Experimental images on ex-vivo biological samples

A pig liver sample was tested with the NCi method. Similarly to what was done in the previous Chapter for diffuse-field elastography, half of the sample was cooked to alter its mechanical properties and the other half was left raw. A 1 ms-long half-sine function was used for the mechanical impulse. Fig. 6.8a. and b. display the spatiotemporal images of the wave fronts obtained for the propagation of the impulse.

In the raw and cooked liver, the Rayleigh wave speeds were respectively 3.6 and 6.3 m/s using time-of-flight. The equivalent wavelength map was computed using the NCi approach, as shown in Fig. 6.8 c. The obtained speed map shows contrast between both parts of the samples with mean equivalent Rayleigh wave speeds of 3.5 m/s and 4.5 m/s. If most elastography techniques assume soft biological tissues as isotropic, the structure of a number of biological soft tissues such as muscle [26], [27], kidney [28], skin [29] or brain tissues [30], [31] makes their mechanical properties dependent on the direction of probing. In these cases, the estimation of the mechanical properties of soft biological tissues can be biased if using an isotropic model [32].

In order to account for soft biological tissues' anisotropic behaviour, a linear transverselyisotropic, quasi-incompressible model in the axis of the tissues' fibers has been proposed as described in Chapter 3 [33]- [35]. This model has been validated by a number of mechanical measurements on soft biological tissues ex-vivo [34], [36], [37]. Experimental imaging of soft biological tissues' anisotropic properties have mainly been performed using ultrasound imaging techniques [28], [29], [38]- [44] or Magnetic Resonance Elastography (MRE) [26], [30], [31], [45]- [48], allowing non invasive measures. These methods have been very successful at providing direction-dependant stiffness assessments for tracking aging [31], characterization of breast tumors [41], [42], [44], [45] or kidney fibrosis [28] but also for diagnostic and clinical treatment of muscle injury or muscle atrophy as a consequence of a pathology [49]- [51].

In optical elastography, experimental studies on mechanical anisotropy have only been performed by rotating the sample to obtain stiffness measurements at different angles [52], [53]. Because classical transient elastography consists in tracking the propagation of a wave front through the sample, only a single direction of propagation can be probed at a time. This limit was recently overcame by recent methods based on diffuse fields [12], [22]. Contrary to classical transient elastography, which is based on a single direction of propagation, with noise-correlation transient elastography, the shearwave mechanical wave fronts are probing the sample in all directions. Consequently, if this approach is combined with a full-field imaging system, this opens the possibility of retrieving in-plane stiffness anisotropy in-real-time without needing to rotate the sample.

Retrieving in-plane stiffness anisotropy images using noise-correlation

Let's consider here an equipartition of noise sources. In the previous Chapter, spatiotemporal correlation were computed in two-dimensions. The average width of the focal spot or the average derivative was used in order to provide stiffness estimations. In the case of an isotropic medium with equipartition of the noise sources, this averaging allowed for an increase in signal as the derivatives and focal width could be extracted from two measurements. However, in the case of transverse in-plane anisotropic medium, both directions should be computed separately to obtain a direction-dependent stiffness estimation.

This was implemented experimentally in the derivation algorithm used previously for noise-correlation in Chapter 5. The MATLAB function imgrad used in isotropic samples was changed to the function imgradxy that computes separate gradient images in both axis of the image. Two shear-wave wavelengths can then be obtained from each pixel, one for each direction. In order to retrieve the main axis of anisotropy, the wavelengths are extracted from the eigen values of the corresponding Hessian matrix.

Considering transverse anisotropy in the plane of imaging, we assume here that the main axis of the ellipse corresponds to the parallel direction in regards to the fibers. Consequently, the two shear-wave wavelengths retrieved using the previous approach are noted λ ⊥ and λ ∥ with λ ∥ corresponding to the direction of the fibers and λ ⊥ to the orthogonal direction. The two wavelength retrieved correspond to the shear-wave speeds c ⊥ and c ∥ introduced in Chapter 3:

c ⊥ = f 0 λ ⊥ c ∥ = f 0 λ ∥ (6.8)
with f 0 the central noise frequency.

The ratio between the shear-wave speeds can be used to define shear-wave anisotropy in two-dimensions. This ratio is referred to here as the anisotropy factor:

f a = c ∥ c ⊥ = λ ∥ λ ⊥ (6.9)
The anisotropy factor is by definition always superior or equal to one: f a = 1 in isotropic media and f a > 1 in the case of anisotropy. This quantity shouldn't be confused with fractional anisotropy, noted FA and usually defined in three spatial dimensions (i = 1, 2, 3) as:

F A = i̸ =j i<j (λ i -λ j ) 2 2 i λ 2 i (6.10)
This quantity is mostly used in volumetric diffuse tensor imaging (see Chapter 2, section 2). In two-dimensional imaging and using the previous notations, its expression is:

F A = 1 2 (λ ∥ -λ ⊥ ) 2 λ 2 ∥ + λ 2 ⊥ = 1 2 (1 - 2f A 1 + f 2 A ) (6.11)
This quantity is between 0 and 1 with 0 corresponding to a perfectly isotropic medium.

Results

Simulation images

In order to validate our approach, the modified algorithm was first tested on simulated movies. A finite difference simulation was ran to create images of a diffuse wave-field in a medium whose properties were anisotropic. Media with gradually increasing in-plane stiffness anisotropy were simulated. The shear wave speed in the vertical direction of the image was fixed at 6 m/s while four different shear-wave speeds were simulated in the horizontal direction (6, 4, 3 and 2 m/s). 500 images of 250 × 250 × pixels were created for each medium. The focal spots are presented for each medium in Fig. 6.9. Figure 6.9: Refocused shear-wave fields at t=0 from finite difference simulation images. An isotropic noise field was generated in homogeneous media with different mechanical properties: the shear wave speed in the y-axis was fixed (6 m/s) and different shear-wave speeds were simulated in the x-direction (6, 4, 3 and 2 m/s respectively for a., b., c. and d.). The corresponding shear-wave speed ratio is indicated on the top of each image.

In order to validate the noise-correlation algorithm's ability to image stiffness anisotropy, a contrasted medium with an isotropic upper part and anisotropic bottom part was simulated. The simulated mechanical properties were the one simulated previously with the correlation functions presented in Fig. 6.9a. and b.. The lambda maps along the vertical and horizontal direction retrieved using the derivation approach were then used to create a map of the fractional anisotropy presented in Fig 6 .10.

As expected, the retrieved focal spots in Fig 6 .9 do present an increase asymmetry with the increasing medium anisotropy. This demonstrates that the effect of medium anisotropy is visible on the two-dimensional cross-correlations on full-field images of a noise field. The retrieved stiffness values and image of Fig. 6.10 show that the modified noise-correlation algorithm for anisotropic media can indeed provide anisotropy contrast on the simulated data.

Experimental results on agarose samples

As seen in the symmetrical point spread functions presented in the previous Chapter, agarose gels have a very isotropic behavior. Creating samples with mechanical anisotropic behavior is challenging as it requires either a structure or a mechanical constraint. In order to modify the behavior of agarose and create anisotropy, glass rods were introduced in the sample with the same orientation to create a structure. These rods were a millimeter in diameter and uniformly distributed within the sample at a few millimeters from its surface. This first experiment shows that the structured agarose samples shows indeed a different correlation function from unstructured agarose samples. The anisotropy from the focal spot corresponds to an anisotropic behavior. This demonstrates that the developed set-up can detect anisotropic in-plane stiffness experimentally on controlled samples.

In order to create anisotropy contrast, cooked rice noodles were also used to structure the agarose samples in different direction. At first a sample was created with the noodles all aligned in one direction, see Fig. 6.12. As the sample could be considered homoge-neous, the focal spots were averaged over the field of view for two types of mechanical noise: fingers tapping (unknown and low frequency) and piezoelectric actuators (chirps of 300 to 1100 Hz). It can be observed that the focal spots present an asymmetry for both sources. Additionally, the sample was rotated in order to check that the anisotropy was indeed from the medium's structure and not from the mechanical sources' positions. The rotation of the focal spots accordingly to the sample's orientation indicates that the anisotropy was indeed characteristic of the medium and not of the source. andd. represent the focal spots calculated from a one-second acquisition using fingers taping as the source of mechanical noise. e., f., and g. represent the focal spots calculated from a one-second acquisition using chirps generated by piezoelectric actuators (300-1100 Hz). The acquisisitons were performed with the noodles' direction horizontal (b. and e.), vertical (c. and f.) and diagonal (d. and g.) on the images.

Anisotropy contrast was then created using the same type of sample. On half of the sample, the noodles were aligned in one direction and on the other half, a perpendicular orientation of the noodles was chosen. The sample is described Fig. 6.13a. The corresponding wavelength maps according to both and each orientation of the noodles are presented in the rest of the Figure . It can be observed on Fig. 6.13b. that the sample seems homogeneous when anisotropy is not considered. The average stiffness is uniform. On the other hand, the wavelengths, and consequently the shear-wave speeds, are very different when considering the axis separately. The wavelengths are bigger in the direction of alignment of the fibers. These results show that the experimental setup can image anisotropic behavior in contrasted samples.

First results in-vivo

To conclude this Chapter, first in-vivo images of stiffness and anisotropy measures are presented in this last experiment. A vibrating phone was used to create a diffuse field in the forearm, as displayed in Fig. 6.14a. In this configuration, two different zones were imaged: one close to the wrist (images b. and d.) and one closer to the elbow (images c. and e.). The stiffness images in the form of wavelength maps averaged over the the different directions of propagation are displayed in b. and c. As the sample is homogeneous, the anisotropy can be calculated over the field of view. The averaged correlation functions are presented in d. and e. for each zone. Shear-wave speed ratio of 1.2 and 2 were calculated for Fig 6 .14 d. and e. respectively. These values are comparable to recently published results [53] where a ratio of 2 was found between the shear wave speeds in perpendicular directions on a similar region of the forearm. These images show the influence of anisotropy on the stiffness estimations performed by averaging the retrieved mechanical wavelengths over the different directions. It also displays the strong anisotropy of biological tissues such as skin and muscle. These first images in-vivo show that the set-up developed during this thesis can successfully achieve full-field stiffness images from a diffuse field as well as anisotropy measurements.

The experimental images are only showing anisotropy measurements from focal spots on homogeneous samples as samples with anisotropy contrast couldn't be obtained yet to create anisotropy images as the one presented in Fig. 6.10 for simulations. This is due to the complexity of manufacturing samples with mechanical anisotropy. Consequently, anisotropic stiffness contrast at this scale is even harder to create from manufactured samples. § 6.

Conclusion

This Chapter presented how the noise-correlation approach could be used in the cases of either a one-sided source or medium anisotropy. The noise-correlation approach presented previously was adapted to one-dimensional sources to perform classical transient elastography image processing. Compared to timeof-flight approaches commonly used, the NCi approach is especially well-adapted to full-field images and heterogeneous media. Experimental images of transient elastography on heterogeneous media could be presented using this approach where time-of-flight could only provide measurements due to the bad SNR. These results are encouraging for further testing of the NCi algorithm on different imaging systems.

The case of in-plane anisotropy was also tackled using a similar approach to take advantage of the full-field setup. The anisotropy of structured samples could be extracted experimentally. First results in-vivo were also presented, including anisotropy measurements as well as stiffness images. These results are however all obtained under the assumption of an equipartition of the noise sources. It can yet be challenging to know experimentally if this hypothesis is valid, especially in the case of natural noise. Without this assumption, the source of asymmetry in the focal spot can stem from the source or the medium, or both.

Simulations ran and displayed in Fig. 6.15 show the separate influence of medium or source anisotropy in noise-correlation. From these images, it seems that computing the correlation function at t ̸ = 0 could separate the source from the medium anisotropy by showing the direction-dependent energy repartition. In the case of an anisotropic sample (e.), the wavefront seem to have the same amplitude regardless of the direction with an ellipsoid shape. On the other hand, in the case where noise sources are onesided (f.), the energy repartition reflects that and seems to be circular in shape. This shows that stiffness images can be highly biased if the condition of source equipartition is not respected in noise-correlation elastography, especially when providing anisotropy measurements. Based on these simulation, it seems that a criteria could be considered to solve that issue and ensure the reliability of the method and/or correct biases due to an anisotropic energy distribution of noise sources. -Chapter 7 -

Conclusion § 7.1 Summary

Elastography in biomedical imaging was first developed based on conventional imaging systems like ultrasounds and MRI. Performing elastography with optical imaging is yet promising for accessing stiffness images quantitatively, with high resolution and without contact. Optical elastography is consequently well-suited for clinical applications in diagnostics as well as in surgery imaging. The presented work is part of the current research investigating the feasibility and perspectives of optical elastography systems for biomedical imaging. In this thesis, a full-field optical imaging system based on digital holography was combined with different elastography approaches to perform quantitative elastography on controlled samples and biological tissues.

Comparing measurements and images with other elastography modalities is tough due to the different frequencies used for the mechanical excitation and the different imaging contrasts between ultrasounds, optics and MRI. The optical set was consequently designed to perform control measurements based on different elastography modalities. The developed optical set-up is highly sensitive and provides single-shot full-field quantitative phase images, compatible with a slow or fast imaging modes. Using different mechanical excitations, various elastography approaches were tested on finite difference simulations and experimentally: harmonic and classical transient elastography for control measurements as well as noise-correlation approaches.

Elastography approaches based on diffuse fields appear as good candidates for clinical translation and are especially well-adapted to the developed full-field imaging system. Compared to the controlled impulses used for classical transient elastography, a diffuse field is indeed easier to obtain with external vibrators or directly from the noise naturally present in-vivo. The proof of concept of noise-correlation elastography was performed on the experimental full-field digital holography system. Stiffness images could be obtained successfully with varying sources of noise and in contrasted biological tissues.

By taking advantage of the full-field set-up, the algorithm developed for diffuse fields was adapted to one-dimensional sources. The method was successfully demonstrated on propagation movies from classical transient elastography with improved image quality in experimental images compared to conventional time-of-flight approach, especially in the case of heterogeneous media.

169

A diffuse field is a superposition of mechanical waves with random and uncorrelated frequency contents, phases, amplitudes and directions in a medium. Consequently, the mechanical excitation probes the sample of interest in all these directions. Since the developed optical system allows displacement imaging in two spatial dimensions, the medium's in-plane stiffness anisotropic properties can also be retrieved. This properties are especially interesting as most biological tissues present anisotropic elastic behavior, which can bias stiffness measurements. In-plane anisotropy measurements were performed experimentally in biological tissues with promising results in-vivo. § 7.

Perspectives

The results presented in this manuscript present a baseline for the future development of an optical elastography system suited for pre-clinical trials. The method and set-up have yet to be refined and developed. Some perspectives are provided here in terms of the elastography method, the optical set-up and the experiments.

As for perspectives concerning the method, a rigorous description of noise-correlation specific to full-field optical imaging would be a great help in considerations of medium and source anisotropy. Indeed noise-correlation methods were first developed for seismology, where a three-dimensional space is considered with large distances between the receivers recording the signal. Noise-correlation theory for soft tissues elastography was mainly developed for ultrasound systems, often based on in-depth images acquired line by-line. In the case of optics, the image acquisition is yet different as it consists in our case of a CMOS sensor acquiring full-field images. Problems specific to full-field imaging emerged during this thesis when adapting existing approaches developed for OCT or ultrasounds to our system. If an analytical expression of the correlation can be obtained in the case of a full-field detection in anisotropic media, it might be possible to assess or even quantify the contribution of the source energy partition in stiffness estimations.

Concerning the experimental optical setup, modifications could be considered to improve the optical signal. The optical sensor of the camera is yet not saturated with a 200µs integration time. This makes the imaging of biological tissues with absorption more challenging. To solve this issue, the power of the laser could be directly increased or alternatives could be found for the optical elements responsible for power losses such as the polarization plates. Another experimental limitation of the set-up is the speckle noise. This noise creates a grain on the resulting stiffness images and a loss of resolution due to numerical image processing. Speckle is inevitable in coherent imaging but with a sufficient laser power, heterodyne interferometry could be further investigated to suppress optically this pattern by inducing a slight frequency shift in the reference arm. The developed set-up can currently provide stiffness images within 4 to 5 seconds using the noise-correlation derivation-based approach. One second accounts for the acquisition time and the rest is image processing and display. This cannot be considered as realtime but the image processing time could be greatly reduced by improving the parallel computing to fully exploit the capacity of the computer.

In order to better assess the physicians needs for specific application of the developed optical elastography system, additional experiments are to be performed with the adaptable and compact version of the system presented in this manuscript. Further testing on anisotropic media and biological tissues should be considered. Investigating a variety of media would be instrumental in assessing the limitations of the current set-up in term of power, sensitivity and stability. Testing on manufactured samples with controlled anisotropic mechanical properties and different sources of noise would also allow to quantify experimentally the accuracy of stiffness estimations with an anisotropic source energy repartition.

At last, discussions with medical professional should be continued on this project to advance towards clinical testings. Discussions have already been held during this thesis over possible clinical applications: monitoring lesions in dermatology, assessing skin grafts and delineating tumors in real-time during tumor removal surgeries. This collaboration should be intensified for the set-up to move towards pre-clinical settings in order to improve the design and customize technical performances for specific applications.

-Appendix A -Résumé étendu en français

Elastographie quantitative en temps réel par holographie numérique

Chapitre 1 -Introduction

La palpation a été utilisée depuis les débuts de la médecine pour évaluer les propriétés mécaniques des tissus biologiques afin de diagnostiquer diverses maladies. Les anomalies de la structure des tissus, comme les tumeurs ou la fibrose présentent la plupart du temps une rigidité plus élevée que dans des tissus sains. Depuis quarante ans, des technologies ont été développées en élastographie pour le biomédical, procurant des mesures de rigidité des tissus biologiques afin d'aider les médecins. Ces techniques ont d'abord été basées sur les ultrasons ou sur l'imagerie par résonance magnétique. Ce n'est qu'à la fin des années 1990 que l'élastographie a été introduite pour la première fois en optique. Les systèmes optiques permettent des mesures très sensibles et sans contact et les récents travaux montrent le potentiel de ces systèmes optiques d'élastographie pour la translation clinique. Ce projet de thèse en élastographie par holographie numérique utilise une approche d'élastographie basée sur les ondes de cisaillement afin d'accéder à une mesure quantitative de la rigidité. Le principe est alors de filmer les ondes de cisaillement dans le champ de vue afin de pouvoir calculer localement leur vitesse. De cette vitesse, on peut alors déduire directement le module de cisaillement du milieu, ce qui caractérise sa rigidité. L'élastographie impulsionnelle classique nécessite le plus souvent une source contrôlée d'onde de cisaillement dans le milieu afin de pouvoir en suivre la propagation, ce qui réduit les applications cliniques possibles. Nguyen et al. ont récemment montré que l'on pouvait s'affranchir de cette limite en adaptant une approche par corrélation de bruit à un système d'élastographie par OCT. Cette méthode est particulièrement prometteuse pour l'application à l'imagerie de tissus biologiques étant donné que le bruit est Chapitre 2 -Imagerie des tissus biologiques Avec l'évolution des technologies et de la science, de nouvelles modalités d'imagerie sont apparues pour les applications biomédicales. Les médecins et les biologistes disposent désormais d'un certain nombre de systèmes fournissant des images avec une gamme de contrastes sur des champs de vue variés. Le second chapitre de cette thèse vise à fournir une introduction basique des systèmes d'imagerie existants afin de placer le travail de cette thèse dans un contexte plus large. Les systèmes d'imagerie médicale présentés combinent souvent un type d'onde pour sonder les tissus (résonance magnétique, ultrasons, lumière visible, rayons X ou gamma) avec un contraste d'imagerie (phase, amplitude, tomographie, fluorescence, Doppler ou imagerie spectrale). Dans la première section, le comportement général des ondes dans les milieux complexes est d'abord introduit afin de poser les bases de la compréhension de la plupart des systèmes d'imagerie. Une brève description de la physique de base qui sous-tend le contraste des systèmes d'imagerie médicale les plus courants tels que les rayons X, les ultrasons, l'imagerie par résonance magnétique (IRM) et l'imagerie nucléaire est ensuite présentée. Les différentes modalités d'imagerie médicale sont présentées ici en fonction du type d'onde utilisé pour l'imagerie des tissus. Les principales caractéristiques, les applications biomédicales et les limites sont discutées pour chaque modalité d'imagerie.

Les principales modalités utilisées dans le diagnostic médical sont présentées dans la

Les sections suivantes se concentrent sur l'imagerie optique. L'optique est principalement utilisée en médecine pour des applications spécifiques telles que les microscopes pour l'anatomopathologie, la tomographie par cohérence optique (OCT) pour l'ophtalmologie ou les caméras pour l'endoscopie. Cependant, l'utilisation de l'optique est récemment apparue comme un outil prometteur pour l'imagerie biomédicale, fournissant de nouveaux contrastes ainsi que des images rapides à haute résolution. Les différentes modalités d'imagerie optique pour les tissus biologiques sont présentées. Les techniques d'imagerie optique basées sur la polarisation, l'absorption ou le contraste de diffusion font l'objet d'une première section. L'imagerie holographique étant au centre de cette thèse, les techniques d'imagerie interférométrique sont spécifiquement présentées dans une autre section. Enfin, la dernière section présente les techniques d'imagerie basées sur la combinaison de l'imagerie acoustique et optique, qu'il s'agisse de combiner les deux ondes ou d'adapter les techniques d'un domaine à l'autre.

Chapitre 3 -Propriétés mécaniques des tissus biologiques

Les propriétés mécaniques des matériaux décrivent leur comportement sous tous types de forces et par conséquent leur stabilité dans des environnements spécifiques. Les propriétés mécaniques sont étroitement liées à la composition et à la structure d'un matériau à différentes échelles. L'évaluation des propriétés mécaniques d'un échantillon peut donc fournir des informations précieuses. Ces informations peuvent être utilisées dans une variété d'applications telles que l'évaluation de nouveaux matériaux pour des applications spécifiques, la cartographie de la structure de la Terre ou l'imagerie biomédicale. Pour les tissus biologiques, une rigidité anormale a été notamment liée à un certain nombre de pathologies. À l'échelle cellulaire, différentes propriétés et mécanismes cellulaires peuvent être mis en évidence. La récupération des informations de rigidité est un outil prometteur pour le diagnostic et le montage développé dans cette thèse vise à développer de telles technologies. Les applications potentielles de l'élastographie en médecine et les motivations de ce travail sont présentées dans une première section.

Après avoir exposé les motivations de ce travail, des grandeurs physiques sont définies afin de qualifier et quantifier les propriétés mécaniques d'un matériau. A partir de la définition de ces grandeurs, des modèles ont été établis pour les tissus biologiques. Ces modèles décrivent le comportement des tissus lorsqu'ils sont soumis à des forces extérieures comme celles illustrées Fig. A.3. La deuxième section de ce chapitre définit les différentes propriétés mécaniques utilisées pour décrire les tissus biologiques et les notions générales nécessaires à la compréhension de l'élastographie, notamment la propagation des ondes dans les tissus et les modèles utilisés pour les tissus biologiques. Sur la base des notions présentées précédemment, différentes approches ont été développées pour déterminer les propriétés élastiques des matériaux. Ces approches d'élastographie reposent toutes sur le même principe. Dans un premier temps, une contrainte mécanique est appliquée sur l'échantillon d'intérêt. Ensuite, la réponse de l'échantillon à cette contrainte est enregistrée. Ensuite, à l'aide des modèles définis précédemment, on peut extraire des propriétés mécaniques telles que la rigidité, l'anisotropie ou la viscosité. Ces différents types d'essais mécaniques et d'imagerie par élastographie sont présentés dans une troisième partie. Leur principe respectif, leurs avantages et leurs limites sont discutés. Grâce aux récents développements de l'élastographie par ultrasons, les médecins disposent désormais de systèmes commerciaux fournissant une image quantitative et non invasive de la rigidité. Si les systèmes à ultrasons sont très utiles pour le diagnostic, leur utilisation nécessite un contact, ce qui est limitant dans un certain nombre de cas, comme pour la chirurgie et les tissus délicats. L'introduction de méthodes optiques pour l'élastographie à la fin des années 1990 a permis des mesures sans contact et précises, particulièrement prometteuses pour des applications en imagerie ophtalmique ou en chirurgie.

Des méthodes d'élastographie plein champ par voie optique ont été étudiées, en utilisant des sources contrôlées de perturbations mécaniques. L'équipe d'Oldenburg [1], [2] a réalisé une élastographie en profondeur plein champ par holographie numérique en utilisant des perturbations harmoniques. A. Nahas et al. [3] ont démontré la faisabilité de l'élastographie impulsionelle plein champ en utilisant l'élastographie par cohérence optique plein champ pour suivre la propagation des ondes de cisaillement en temps réel. L'objectif de cette thèse est de concevoir un système d'imagerie expérimental pour réaliser des images de rigidité plein champ. En choisissant l'imagerie optique, un système d'imagerie rapide, sans contact et à haute résolution est construit. Basé sur l'interférométrie, le dispositif permet de réaliser des images de déplacement à haute sensibilité. Les images résultantes peuvent être exploitées pour réaliser des images d'élastographie par ondes de cisaillement avec une excitation mécanique de faible puissance. Le système d'élastographie optique est présenté et caractérisé ainsi que les différentes approches d'élastographie réalisées dans les premières sections.

Afin d'assurer l'aspect quantitatif des approches développées, les résultats doivent être comparés. Or, les fréquences d'excitation utilisées dans les essais mécaniques (≈ 1 Hz) et les systèmes à ultrasons commerciaux (≈ 50Hz) sont très différentes de celles utilisées ici (≈ 1 kHz). Une étude précédente [4]

Chapitre 7.1 -Conclusion

L'élastographie en imagerie biomédicale a d'abord été développée à partir de systèmes d'imagerie conventionnels comme les ultrasons ou l'IRM. L'élastographie par imagerie optique est pourtant prometteuse pour accéder aux images de rigidité de manière quantitative, avec une haute résolution et sans contact. L'élastographie optique est par conséquent bien adaptée aux applications cliniques de diagnostic et d'imagerie chirurgicale. Le travail présenté ici fait partie de la recherche actuelle portant sur la faisabilité et les perspectives des systèmes d'élastographie optique pour l'imagerie biomédicale. Dans cette thèse, un système d'imagerie optique plein champ basé sur l'holographie numérique a été combiné à différentes approches d'élastographie pour réaliser des mesures de rigidité quantitatives sur des échantillons contrôlés ainsi que des tissus biologiques.

La comparaison des mesures et des images avec d'autres modalités d'élastographie est difficile en raison des différentes fréquences utilisées pour l'excitation mécanique et des différents contrastes d'imagerie entre les ultrasons, l'optique et l'IRM. Le montage optique a donc aussi été conçu pour effectuer des mesures de contrôle basées sur différentes modalités d'élastographie. Le dispositif optique mis au point est très sensible et fournit des images de phase quantitatives plein-champ en une seule prise, compatibles avec des modes d'imagerie lent ou rapide. En utilisant différentes excitations mécaniques, diverses approches d'élastographie ont été testées sur des simulations par différences finies mais aussi expérimentalement : élastographie harmonique et impulsionelle classique pour les mesures de contrôle ainsi que des approches de corrélation du bruit.

Les approches d'élastographie basées sur les champs diffus semblent être de bons candidats pour une application clinique et sont particulièrement bien adaptées au système d'imagerie développé, qui permet d'obtenir des images plein-champ. Par rapport aux impulsions contrôlées utilisées pour l'élastographie impulsionelle classique, un champ diffus est en effet plus facile à obtenir avec des vibrateurs externes ou directement à partir du bruit naturellement présent in-vivo. La preuve de concept de l'élastographie par corrélation de bruit a été réalisée sur le système expérimental d'holographie numérique plein champ. Les images de rigidité ont pu être obtenues avec succès avec différentes sources de bruit et dans des tissus biologiques contrastés.

En exploitant l'imagerie plein-champ offerte par le système experimental, l'algorithme développé pour les champs diffus a été adapté aux sources unidimensionnelles. La méthode a été démontrée avec succès sur des films de propagation issus de l'élastographie impulsionelle classique avec une meilleure qualité d'image dans les images expérimentales par rapport à l'approche conventionnelle du temps de vol, en particulier dans le cas de milieux hétérogènes.

Un champ diffus est une superposition d'ondes mécaniques dont le contenu en fréquence, les phases, les amplitudes et les directions sont aléatoires et non corrélés dans un milieu. Par conséquent, dans le cas d'un champs diffus, l'excitation mécanique sonde l'échantillon d'intérêt dans toutes ses directions. Puisque le système optique développé permet l'imagerie du déplacement dans deux dimensions spatiales, les propriétés anisotropes de la rigidité dans le plan du milieu peuvent également être extraites des acquisitions experimentales. Ces propriétés sont particulièrement intéressantes car la plupart des tissus biologiques présentent un comportement élastique anisotrope, ce qui peut fausser les mesures de rigidité. Des mesures d'anisotropie dans le plan ont été réalisées expérimentalement dans des tissus biologiques avec des résultats prometteurs in-vivo.

Chapitre 7.2 -Perspectives

Les résultats présentés dans ce manuscrit constituent une base de référence pour le développement futur d'un système d'élastographie optique adapté aux essais précliniques. La méthode et le dispositif doivent encore être affinés et développés. Quelques perspectives concernant la méthode d'élastographie, le montage optique et les expériences sont présentées ici.

En ce qui concerne les perspectives relatives à la méthode, une description rigoureuse de la corrélation du bruit spécifique à l'imagerie optique plein champ serait d'une grande aide dans les considérations sur l'anisotropie du milieu et de la source. En effet, les méthodes de corrélation du bruit ont d'abord été développées pour la sismologie, où un espace tridimensionnel est considéré avec de grandes distances entre les récepteurs enregistrant le signal. La théorie de corrélation du bruit pour l'élastographie des tissus mous a été principalement développée pour les systèmes à ultrasons, souvent basés sur des images en profondeur acquises ligne par ligne. Dans le cas de l'optique, l'acquisition d'images est encore différente puisqu'elle consiste dans notre cas en un capteur CMOS acquérant des images plein champ. Des problèmes spécifiques à l'imagerie plein champ sont apparus au cours de cette thèse lors de l'adaptation à notre système d'approches existantes développées pour l'OCT ou les ultrasons. Si une expression analytique de la corrélation peut être obtenue dans le cas d'une détection plein champ dans des milieux anisotropes, il pourrait être possible d'évaluer ou même de quantifier la contribution de la partition de l'énergie de la source dans les estimations de la rigidité.

En ce qui concerne le dispositif optique expérimental, des modifications pourraient être envisagées pour améliorer le signal optique. Le capteur optique de la caméra n'est pas encore saturé avec un temps d'intégration de 200 µ s. Cela rend l'imagerie des tissus biologiques avec absorption plus difficile. Pour résoudre ce problème, la puissance du laser pourrait être directement augmentée ou des alternatives pourraient être trouvées pour les éléments optiques responsables des pertes de puissance tels que les polarisateurs et l'analyseur. Une autre limite expérimentale du montage est le bruit de speckle. Ce bruit crée un grain sur les images de rigidité obtenues et une perte de résolution due au traitement numérique des images. Ce chatoiement est inévitable en imagerie cohérente, mais avec une puissance laser suffisante, l'interférométrie hétérodyne pourrait être étudiée plus en détails pour supprimer optiquement ce motif en induisant un léger décalage de fréquence dans le bras de référence. Le dispositif développé peut actuellement fournir des images de rigidité en 4 à 5 secondes en utilisant l'approche basée sur la corrélation de bruit par dérivation. Une seconde correspond au temps d'acquisition, le reste étant consacré au traitement et à l'affichage des images. Cela ne peut être considéré comme du temps réel, mais le temps de traitement de l'image pourrait être considérablement réduit en améliorant le calcul parallèle pour exploiter pleinement la capacité de l'ordinateur.

Afin de mieux évaluer les besoins des médecins pour une application spécifique du système d'élastographie optique développé, des expériences supplémentaires doivent être réalisées avec la version adaptable et compacte du système présenté dans ce manuscrit. Des tests supplémentaires sur des milieux anisotropes et des tissus biologiques doivent être envisagés. L'étude d'une variété de milieux permettrait d'évaluer les limites de la configuration actuelle en termes de puissance, de sensibilité et de stabilité. Des essais sur des échantillons fabriqués avec des propriétés mécaniques anisotropes contrôlées et différentes sources de bruit permettraient également de quantifier expéri-mentalement la précision des estimations de la rigidité avec une répartition anisotrope des sources d'énergie.

Enfin, les discussions avec les professionnels de la santé doivent être poursuivies sur ce projet pour avancer vers des tests cliniques. Des discussions ont déjà eu lieu au cours de cette thèse sur les applications cliniques possibles : suivi des lésions en dermatologie, évaluation des greffes de peau et délimitation des tumeurs en temps réel lors des opérations d'ablation de tumeurs. Cette collaboration doit s'intensifier pour que le montage évolue vers des contextes précliniques afin d'améliorer le design et de personnaliser les performances techniques pour des applications spécifiques. ) and (d.) for the 4ms-long pulse, (c.) and (e.) for the 1.33 ms-long pulse), chirp function (f.) and noise (g.). 148 6.4 Speed maps obtained from simulated propagation movies (images of 500×500 pixels of 100 µm², 1000 frames at 50 kHz). The medium of propagation was a star-shaped inclusion in a softer background, illustrated by the map (a.). Using time-of flight ((b.) and (c.)) and NCi algorithm ((d.) to (g.)).

Results with different profiles of one-dimensional shear-wave sources are presented: half-sine functions of 4 and 1.33 milliseconds, a chirp function and a noise signal. Profiles of the excitation sources are represented above each speed map: half sine functions ((b.) and (d.) for the 4ms-long pulse, (c.) and (e.) for the 1.33 ms-long pulse), chirp function (f.) and noise (g.).

6.5 Curve representing the average speed of Rayleigh waves measured at the surface of scattering agarose samples of gradually increasing stiffness (from 0.5 to 2.0% agarose mass concentration) using the FWHM approach. Time-of-flight measures are used as a reference and presented on the xaxis. A 1 ms-long half-sine excitation function was used to generate the mechanical waves. The measures were averaged over six acquisitions for each sample. The speed map retrieved from the same set of propagation images using the NCi approach is shown in (c). . . . . . . . . . . . . . . . . . . . . . . . 6.9 Refocused shear-wave fields at t=0 from finite difference simulation images. An isotropic noise field was generated in homogeneous media with different mechanical properties: the shear wave speed in the y-axis was fixed (6 m/s) and different shear-wave speeds were simulated in the xdirection (6, 4, 3 and 2 m/s respectively for a., b., c. and d.). The corresponding shear-wave speed ratio is indicated on the top of each image.

Figure 2 . 1 :

 21 Figure 2.1: Schematic of the acoustic and the electromagnetic spectra with examples of their applications.

Figure 2 . 2 :

 22 Figure 2.2: Table of comparison of the different medical imaging modalities in terms of scale, cost and ionizing radiations.

Figure 2 . 3 :

 23 Figure 2.3: Schematic of the different interactions of waves with matter.

Figure 2 . 4 :

 24 Figure 2.4: Schematic of the refraction and reflection at an optical interface.

Figure 2 . 5 :

 25 Figure 2.5: Illustration of Rayleigh, Mie and geometrical scattering regimes.

Figure 2 . 6 :

 26 Figure 2.6: Illustration of the different scattering regimes.

Figure 2 . 7 :

 27 Figure 2.7: Image known as the first x-ray image, presented in 1895.

Figure 2 . 8 :

 28 Figure 2.8: Mass coefficient of different interactions of x-ray radiations with matter depending on the beam's energy. The area for x-ray medical imaging is highlighted in gray.

Figure 2 . 9 :

 29 Figure 2.9: Schematic of the photoelectric effect and Compton scattering.

Figure 2 .

 2 Figure 2.10: X-ray images of a mouse in-vivo with conventional attenuation (a), differential phase-contrast based on x-ray refraction (b) and dark-field (c) based on x-ray scattering. The white bars correspond to 1 cm. [9]

Figure 2 . 11 :

 211 Figure 2.11: Schematic of different acquisition modes of a medical ultrasound system. z designates the in-depth spatial direction, t the time and x a lateral direction. The gradual shading in the B and M modes indicate how the successively acquired signals can contribute to a final full-field image.

Figure 2 . 13 :

 213 Figure 2.13: Images of a jugular vein and carotid artery from a clinical ultrasounds system in B-mode superimposed with Doppler contrast [13].

Figure 2 . 14 :

 214 Figure 2.14: Schematic of particles magnetic spins in the absence and presence of a magnetic field.

Figure 2 . 16 :

 216 Figure 2.16: Schematic of the radioactive decay leading to the emission of gamma photons exploited for SPECT and PET imaging.

  2.3.4.2 Medical use and accessible contrasts

Figure 2 . 17 :

 217 Figure 2.17: Different approaches in optics theory with illustrations and main equations.

Figure 2 . 18 :

 218 Figure 2.18: Absorption spectra of several contrast agents in biological tissue optical imaging at normal concentrations from[29]

  .19.

Figure 2 . 19 :

 219 Figure 2.19: Illustration of the main components of a cell (a.) [30] and their characteristic sizes in regards to Mie and Rayleigh scattering (b.).

Figure 2 . 21 : 1 =

 2211 Figure 2.21: Illustration of three states of polarization of light from [2].

Figure 2 . 22 :

 222 Figure 2.22: Energy processes after photon absorption by a molecule.

Figure 2 . 23 :

 223 Figure 2.23: Indocyanine Green emission and absorption spectra from[33]

Figure 2 . 24 :

 224 Figure 2.24: Basic optical set-up for fluorescence imaging.

Figure 2 . 25 :

 225 Figure 2.25: Example of some quenching processes

Figure 2 . 26 :

 226 Figure 2.26: Some biomedical optical imaging methods and their typical penetration depth and spatial resolution. Values from [35].

  2.27). Temporal measurements can be applied to various optical imaging methods, including fluorescence imaging.

Figure 2 . 27 :

 227 Figure 2.27: Paths of ballistic, snake and diffusion photons propagating through tissue.For a short pulse of light transmitted through tissue, the ballistic component emerges fist, followed by the diffuse component. The early par of the diffuse component is known as the "snake" component because the path of these photons approximate straight lines in the tissue. Image from from[30]

Figure 2 . 28 :

 228 Figure 2.28: Comparison of wide-field and confocal microscopy on a mouse intestine section from [34].

Figure 2 . 29 :

 229 Figure 2.29: Two-photon and Second Harmonic Generation energy processes

Figure 2 . 30 :

 230 Figure 2.30: Contrast enhancing techniques in optical microscopy. Images from [39].

Figure 2 . 31 :

 231 Figure 2.31: Detailed electromagnetic spectrum from UV to IR ranges.

Figure 2 . 32 :

 232 Figure 2.32: Scattering phenomena used in distributed optical fiber sensors from [40].

Figure 2 . 33 :

 233 Figure 2.33: Schematic of the interaction of two monochromatic plane waves.

2⟨E

Figure 2 . 34 :

 234 Figure 2.34: Schematic of the Young's experiment with two pinholes.

  54) with sinc(x) = sin(x)/x represented in Fig 2.35.

Figure 2 . 35 :

 235 Figure 2.35: Representation of the sinc function.

  62) with k the mean wave number, I(O i ) the average intensity at O i and R i the distance between O i and the central point S of the surface element dσ of the source, as illustrated Fig.2.36.

Figure 2 . 36 :

 236 Figure 2.36: Schematic illustrating the notations of Van-Cittert-Zernike theorem

Figure 2 . 37 :

 237 Figure 2.37: Images and histogram of speckle pattern acquired in phase and amplitude using the digital holography set-up described in Chapter 4 in (a.) amplitude and (b.) phase.

Figure 2 . 38 :

 238 Figure 2.38: Schematic illustrating the configuration considered for characterizing speckle: propagation in free space (a.) and through an optical system (b.)

Figure 2 . 39 :

 239 Figure 2.39: Illustration of a Michelson interferometer in three configurations: (a.) Parallel-faces, (b.)Tilted-faces and (c.) Sample analysis.

Figure 2 . 40 :

 240 Figure 2.40: Illustration of a Mach-Zehnder interferometer.

Figure 2 . 41 :

 241 Figure 2.41: Schematic of a Michelson-based OCT system. The low coherence illumination light beam is divided into two by the beam splitter. The reflected beam goes to illuminate the sample while the transmitted beam is reflected onto a moving mirror.Both beams recombine on the camera to form an interference pattern. This interference pattern varies depending on the optical path difference between both beams. For low coherence light, the interference signal is maximal when this optical path difference is null and is indistinguishable when the path difference is greater than the coherence length of the source.

Figure 2 . 42 :

 242 Figure 2.42: Temporal interference fringe pattern detected as a result of reference mirror motion. Image from [45].

Figure 2 . 43 :

 243 Figure 2.43: Schematic of in-line hologram formation.

Figure 2 . 44 :

 244 Figure 2.44: Schematic of the coordinate system for in-line holography.

Figure 2 . 45 :

 245 Figure 2.45: Schematic of an hologram formation set-up (a.) and the virtual image reconstruction set-up (b.) and the real-image reconstruction set-up (c.).

Figure 2 . 46 :

 246 Figure 2.46: Schematic of the hologram formation in an off-axis configuration. The different component of the hologram are consequently shifted in the spatial frequency domain. They are completely separated if the angle and numerical aperture of the optical system are appropriately chosen.

Figure 2 . 47 :

 247 Figure 2.47: Notations for free space propagation between two planes with transverse coordinates (ξ, η) and (x, y)

Figure 2 . 48 :

 248 Figure 2.48: Illustration of acousto-optic imaging in thick diffuse biological tissues.

Figure 2 . 49 :

 249 Figure 2.49: Non invasive photo-acoustic tomography (PAT) imaging of the superficial layer of the in-vivo rat brain acquired with the skin and skull intact (a) and open-skull photograph of the ray brain surface acquired after the PAT experiment[77]

Figure 2 . 50 :

 250 Figure 2.50: Retinal angiogram using intensity-based Doppler variance OCT, adapted from [82] .

Figure 3 . 1 :

 31 Figure 3.1: Examples of different types of stresses.

Figure 3 . 3 :

 33 Figure 3.3: Schematic representing a classification of mechanical behaviors in solids.

Figure 3 . 4 :

 34 Figure 3.4: Schematic introducing the notations for shear and elastic modulus in a linear elastic solid.

Figure 3 . 5 :

 35 Figure 3.5: Basic models for viscoelastic behavior.

Figure 3 . 6 :

 36 Figure 3.6: Illustration of a compression elastography set-up based on an ultrasound probe, published in [27].

  Commercial ultrasound elastography systems have been first developed by specialized companies (Fibroscan TM by Echosens in 2003, Aixplorer TM by Supersonic Imagine in 2008). In the following years, elastography was introduced to the imaging systems of the dominant manufacturer (Logiq E9 by GE Healthcare in 2014, Virtual Touch TM by Siemens Healthineers in 2013 or ElastPQ TM by Philips Healthcare in 2013).

Figure 3 . 7 :

 37 Figure3.7: Schematic of the eye from[START_REF] Schedin | Pulsed digital holography for deformation measurements on biological tissues[END_REF]

Figure 3 . 8 :

 38 Figure 3.8: Image of transient displacement of a propagating mechanical wave in ex-vivo porcine cornea with 10 mmHg (a.) and 40 mmHg (b.) of IOP using an a phase-sensitive OCT system and non-contact stimulation based on ARF, published in[55]

Figure 3 . 9 :

 39 Figure 3.9: Image of ultrasound transient elastography on the liver of patients different fibrosis grade.The slope on each image represents the propagation velocity and is consequently linked to the liver stiffness, published in[START_REF] Sandrin | Transient elastography: A new noninvasive method for assessment of hepatic fibrosis[END_REF]

Figure 3 . 10 :

 310 Figure 3.10: Image of Magnetic resonance elastography on the liver of patients with different grades of fibrosis and with obesity or ascites, published in[44]

Figure 3 . 11 :

 311 Figure 3.11: Image of ultrasound elastography on a breast with grade I infiltrating ductal carcinoma performed with an Aixplorer (SuperSonic Imagine, France) prototype published in [94].

Figure 3 . 12 :

 312 Figure3.12: Microscopy images of stained breast tumor margin sections with zero (A), 1 mm (B) and 2 mm margins (C), published in[START_REF] Brown | Optical assesssment of tumor resection margins in the breast[END_REF]

Figure 3 . 13 :

 313 Figure 3.13: (a.) Wide-field en-face OCT images of a freshly excised lumpectomy with the elastography image superimposed. (b.) and (e.) Histology acquired in orthogonal plane to OCT from regions indicated by green dashed line and red dashed line, respectively.(c.) and (f.) En-face OCT image corresponding to a zoom of (a.) in regions indicated by green dashed line and red dashed line, respectively. (d.) and (g.) show the overlaid elastography image corresponding to images (c.) and (f.). A, adipose; S, stroma; and T, tumor. Image published in[START_REF] Allen | Wide-field optical coherence microelastography for intraoperative assessment of human breast cancer margins[END_REF]

Figure 3 . 14 :

 314 Figure 3.14: Images of endoscopic view (A) and OCT (B) obtained in a swine bowel in-vivo, published in [107]

Figure 4 . 1 :

 41 Figure 4.1: Schematic of the holographic off-axis set-up and the stroboscopic approach. The optical set-up is represented on the right of the image: the object beam (in red on the image) is back-scattered by the sample before interfering with the reference beam (in purple) (f L1 = 100 mm, f L2 = 150 mm). On the left of the image, the synchronization of the camera and the mechanical actuator for the fast mode is detailed: the impulses represented in green trigger the piezoelectric actuator and the blue signal activates the camera.

Figure 4 . 2 :

 42 Figure 4.2: Processing for the 2D Fourier transform method to retrieve phase and amplitude images from each camera frame.

Figure 4 . 3 :

 43 Figure 4.3: Photo (a) and corresponding schematic (b) of the compact off-axis digital holography set-up with an adaptable objective. The reference beam and the object beam are represented respectively in purple and red. The objective presented on the photo consists of two 100mm focal distance lenses with a diaphragm.

Figure 4 . 4 :

 44 Figure 4.4: a. Frame captured by the CMOS captor with standard experimental condition (integration time of 200µs) without the laser source. The histograms presented in b., c., d. and e. correspond respectively to the top left, top right, bottom left and bottom right quadrant of the sensor. The total histogram of the captor is displayed in f. The average values (Mean) and standard deviation (StDev) are given for each curve.

Figure 4 . 5 :

 45 Figure 4.5: Histogram (a.) and standard deviation (b.) in time averaged over 1440 ×1440 pixels and 500 frames for an typical acquisition in ambient lighting without any object (200µs, 500 frames at 500Hz).

Figure 4 . 6 :

 46 Figure 4.6: Noise measurements performed with the experimental set. The sensor was exposed to gradual intensities. The intensity fluctuations (in gray levels) were traced in function of I/ √ I (in gray levels -1/2). The curve was fitted with an affine line with R 2 = 0, 97.

Figure 4 . 7 :

 47 Figure 4.7: Image auto correlation of a typical acquisition in an agarose sample to determine the speckle's size. (a.) Raw amplitude image, (b.) averaged auto-correlation, (c.) profile of the auto-correlation in blue fitted with the red curve to estimate FWHM.Algorithm by[7].

Figure 4 . 8 :

 48 Figure 4.8: Experimental phase (a.) and phase difference (b.) images of a typical diffuse noise field acquired on agarose. The red square is a zoom on a 100 × 100 pixels zone to better see the speckle pattern. The histogram of the phase difference image is displayed in image (c.), where the central peak represents the useful signal while the secondary peaks around ∓2π are the artifacts due to speckle noise. Image (d.) is the same as image (b.) with a better suited color-map to compare with the filtered images (e.) and (f.) using respectively a median filtering and a threshold and Gaussian filtering. The threshold used is illustrated in (c.) by blue dotted lines.

Figure 4 . 9 :

 49 Figure 4.9: Image of an USAF resolution target taken with the presented optical setup.(a.) is the full-field image and (b.) is a zoom centered on the resolution limit of the presented optical system.

Figure 4 . 10 :

 410 Figure 4.10: Shear-wave speed map calculated from the simulation movie of a sinusoidal impulse 4 ms and 1.3 ms long propagating through a star-shaped soft inclusion. The simulated shear-wave speeds are 3.0 m/s for the inclusion and 6 m/s for the background. In a. and b., only the time-of-flight calculation was used and in c. and d., a propagative filter was added prior to the time-of-flight processing.

Figure 4 . 11 :

 411 Figure 4.11: Acquisitions of harmonic excitation on agarose samples. On (a.), each curve represents a sample with a different agarose concentration. The measured wavelength is computed in function of the excitation frequency. At 1 kHz, images illustrating the different observed frequencies depending on agarose concentration are shown in (b.). These images show the same field of view and correspond to the gray vertical band in the curve. The images were acquired using a digital holography set-up. Gaussian filtering was used on the images.

Figure 4 . 12 :

 412 Figure 4.12: Measured shear-wave speed in different agarose samples calculated from the measured wavelength in function of the excitation frequency and agarose mass concentration.

Figure 4 . 13 :

 413 Figure 4.13: Curve representing the measured shear-wave speeds in homogeneous samples with varying agarose mass concentration. Each point was computed on 5 different acquisitions. The mean square error is displayed using vertical error bars. 1-ms-long mechanical impulses were generated in the samples. To illustrate these results, a. to g. are examples of spatiotemporal profiles from experimental propagation movies acquired in homogeneous agarose samples with gradually varying mass concentration (0.5, 0.75, 1.0, 1.25, 1.50, 1.75 and 2.0 % respectively). The x-axis represents the consecutive frames of the movie and the y-axis is a line on the image along the propagation axis. Directional and Gaussian filtering was used on the images. All the images were acquired using a digital holography set-up and a stroboscopic approach to reach 25 kHz equivalent frame rate.

Figure 4 . 14 :

 414 Figure 4.14: Correspondence between speed measurements on phase and amplitude propagation movies acquired using a digital holography set-up on samples with varying agarose mass concentrations.

  4.15. The change of slope in the image shows the two different speeds recorded in the two different zones of the sample. Speeds of 4.3 and 1.5 m/s are recorded for each side of the sample. These values are coherent with the results on homogeneous samples, even though deviations in values can be attributed to the sample's mechanical properties' susceptibility to temperature or humidity.

Figure 4 . 15 :

 415 Figure 4.15: Speedmap from the propagation of a shear 1-ms-long impulse in an agarose sample where half has 0.5 % and half 1.0 % agarose mass concentration, as illustrated on the left. This image was acquired using a digital holography set-up and a stroboscopic approach to reach 25 kHz equivalent frame rate. Directional and gaussian filtering was used on the image.

Figure 4 . 16 :

 416 Figure 4.16: Images from the propagation of a 1-ms-long mechanical impulse in an inclusion sample. A white light image of the sample's surface in shown in a. and the inclusion is highlighted in red. b. represents frames from the propagation movie when the wave fronts are highlighted in orange when they propagate in the background and in red when they propagate through the inclusion.

Figure 4 .

 4 Figure 4.17: a. and b. presents the recording of the propagation of a 1 ms mechanical impulse through a cooked and raw porcine liver sample recorded with the stroboscopic digital holography experimental set-up at an equivalent frame rate of 25 kHz.

Figure 4 . 18 :

 418 Figure 4.18: Cropped frames from the propagation movie are presented for the nerve in free space under no stress and tensile stress, respectively a. and b. The corresponding speeds are represented in c. and d. These images were acquired using the amplitude acquired using a stroboscopic digital holography set-up at an equivalent 1.25 MHz frame rate.

Figure 4 . 19 :

 419 Figure 4.19: Cropped frames from a propagation movie acquired in amplitude with a digital holography set-up at an equivalent 1.25 MHz frame rate using a sinusoidal source on a murine nerve. The images displayed correspond to differences between successive frames.

Figure 5 . 1 :

 51 Figure 5.1: Schematic of the described time-reversal cavity.

17

 17

  ) with e C the source contribution from C.

Figure 5 . 2 :

 52 Figure 5.2: Illustration of the focal spot size method for noise-correlation image processing on simulation images using a finite difference algorithm. a. represents the modeled medium, which if is a stiff inclusion. Images of the displacement field in the medium in the presence of noise sources are then simulated (b.). The correlation function in one-dimension is presented in c. for a point in the stiff inclusion and in the softer background.

Figure 5 . 3 :

 53 Figure 5.3: Mean wavelengths given by the noise correlation method depending on the pulse speed measured by the control experiment. The red curve is the theoretical curve for a excitation of central frequency 1.1 kHz. These data correspond to agarose samples with different concentrations. Five different one-second acquisitions were made for each sample with a 200 µm exposure time. The mean PSF of the mechanical waves are illustrated for each sample (a, b, c, d and e respectively for 1%, 1.25%, 1.5%, 1.75%, and 2% agarose) in 600 × 600 pixels images

Figure 5 . 4 :

 54 Figure 5.4: Images from a diffuse field acquisition on a 1% inclusion in a 2% agarose sample with 1% titanium dioxide concentration. (a.) Description of the inclusion sample and mechanical noise sources. (b.) Amplitude image of the sample in ambient lighting. the inclusion zone is lined with red dotted lines. (c.) Phase displacement images on 4 consecutive frames. (d.)The wavelength map image is retrieved using the noise correlation method. Averaged images of the correlation functions for a zone in the softer background (e.) and for a zone in the inclusion (f.).

Figure 5 . 6 :

 56 Figure 5.6: Images from a diffuse field acquisition on a 1% inclusion in a 2% agarose sample with 1% titanium dioxide concentration. (a.) Description of the inclusion sample and mechanical noise source. (b.) The wavelength map image is retrieved using the noise correlation method. (c.) The frequency content of the mechanical excitation was computed from the acquisition to extract the central noise-frequency.

c

  shear-wave = F W HM k [C(r, r 0 , t = 0)] F W HM [C(r = r 0 , r 0 , t)] . (6.4) with F W HM [] the FWHM of a 1D function.F W HM k [] corresponds to the FWHM of a 2D function in the propagation direction k.

Figure 6 . 1 :

 61 Figure 6.1: Schematic of the time-of-flight and NCi algorithm for retrieving shear-wave speed estimation from full-field acquisitions.

  4 and 5.7 m/s for a soft star-shaped inclusion. Stiffer star-shaped inclusions were also generated for a 3 m/s background. Similarly to homogeneous simulations, different temporal profiles and duration of mechanical sources were tested. Typical propagation movies are showcased in Fig. 6.2.

Figure 6 . 2 :

 62 Figure 6.2: Frames from simulated propagation movies (500×500 pixels of 100 µm² with 1000 frames at 50 kHz) displayed in Video 1 (MP4, 0.20 MB) and Video 2 (MP4, 0.29 MB). A soft star-shaped inclusion is simulated with 3 m/s shear-wave speed for the inclusion and 6 m/s for the background . (a), (b), (c) and (d) correspond to frames 350, 450, 550 and 650 of Video 1, for which the shear-wave source function is a 4 ms-long half-sine. (e), (f ), (g) and (h) correspond to frames 350, 450, 550 and 650 of Video 2, where the source function is 4 ms-long a noise signal (400-800 Hz).

Figure 6 . 3 :

 63 Figure 6.3: Speed maps obtained from simulated propagation movies (images of 500×500 pixels of 100 µm², 1000 frames at 50 kHz). The medium of propagation was a star-shaped inclusion in a softer background, illustrated by the map (a.). Using time-of flight ((b.) and (c.)) and NCi algorithm ((d.) to (g.)). Results with different profiles of one-dimensional shear-wave sources are presented: half-sine functions of 4 and 1.33 milliseconds, a chirp function and a noise signal. Profiles of the excitation sources are represented above each speed map: half sine functions ((b.) and (d.) for the 4ms-long pulse, (c.) and (e.) for the 1.33 ms-long pulse), chirp function (f.) and noise (g.).

Figure 6 . 4 :

 64 Figure 6.4: Speed maps obtained from simulated propagation movies (images of 500×500 pixels of 100 µm², 1000 frames at 50 kHz). The medium of propagation was a star-shaped inclusion in a softer background, illustrated by the map (a.). Using time-of flight ((b.) and (c.)) and NCi algorithm ((d.) to (g.)). Results with different profiles of one-dimensional shear-wave sources are presented: half-sine functions of 4 and 1.33 milliseconds, a chirp function and a noise signal. Profiles of the excitation sources are represented above each speed map: half sine functions ((b.) and (d.) for the 4ms-long pulse, (c.) and (e.) for the 1.33 ms-long pulse), chirp function (f.) and noise (g.).

Figure 6 . 5 :

 65 Figure 6.5: Curve representing the average speed of Rayleigh waves measured at the surface of scattering agarose samples of gradually increasing stiffness (from 0.5 to 2.0% agarose mass concentration) using the FWHM approach. Time-of-flight measures are used as a reference and presented on the x-axis. A 1 ms-long half-sine excitation function was used to generate the mechanical waves. The measures were averaged over six acquisitions for each sample. The standard deviation is presented and the gray line represents the x=y function.

Figure 6 . 6 :

 66 Figure 6.6: Curves representing the average speed of Rayleigh waves measured at the surface of scattering agarose samples of gradually increasing stiffness (from 0.5 to 2.0% agarose mass concentration) using the NCi derivation-based approach. Time-of-flight measures were used as a reference and presented on the x-axis of each graph. (a) and (b) were obtained respectively using four 1 ms-long half-sine excitation function and using a 2 ms-long chirp excitation function, from 500 to 2 kHz. The measures are averaged over six acquisitions for each sample. The standard deviation is also presented. The gray line represents the x=y function.

Figure 6 . 7 :

 67 Figure 6.7: Images obtained from a 2.0% agarose sample with a 1.0% star-shaped agarose inclusion. (a.) describes the sample and mechanical source. (b.) is an amplitude image of the sample. The inclusion cannot be distinguished by its optical properties in this image and its shape and position is highlighted by a red dotted line. (c.) shows the raw phase difference images obtained with the digital holography set-up in stroboscopic mode (25 kHz equivalent frame rate). Frames 108, 179, 250 and 321 of the 500 are displayed. (d.) is the speed-map reconstruction from the acquired propagation movie achieved with the NCi algorithm.

Figure 6 . 8 :

 68 Figure 6.8: The spatiotemporal representations of a 1 ms-long pulse propagating in a cooked and raw pork liver sample are presented in (a) and (b) respectively. The speed map retrieved from the same set of propagation images using the NCi approach is shown in (c).

Figure 6 . 10 :

 610 Figure 6.10: Images from simulations. a. Images from a movie simulating the displacements generated by a noise field in a dual sided anisotropic medium. f a = cy cx = 2/3 and 3/2 on the upper and bottom side respectively. b. and c. are the focal spots reconstructed for each side. The schematics represent the calculation of the anisotropy direction and ratio from these focal spots. d. Reconstructed anisotropy ratio map obtained with the derivation approach.

  Finger tapping was used to create a unknown low-frequency noise field. Results of this experiment are presented Fig. 6.11. The averaged focal spot profile is presented for this structured sample in d. and for a homogeneous agarose samples in b.

Figure 6 . 11 :

 611 Figure 6.11: Agarose sample description (a.) and (c.) and corresponding reconstructed time-reversed-field at t=0 (b.) and (d.) from experimental data using noise correlation with fingers tapping as the source of the diffuse field. The experimental set-up was a full-field digital holography imaging system. The first sample presented in (a.) is a homogeneous agarose sample. In the second sample (c.), glass rods were introduced to create an anisotropic behavior.

Figure 6 .

 6 Figure 6.12: a. Schematic of the agarose sample structured with rice noodles. b., c.,and d. represent the focal spots calculated from a one-second acquisition using fingers taping as the source of mechanical noise. e., f., and g. represent the focal spots calculated from a one-second acquisition using chirps generated by piezoelectric actuators (300-1100 Hz). The acquisisitons were performed with the noodles' direction horizontal (b. and e.), vertical (c. and f.) and diagonal (d. and g.) on the images.

Figure 6 .

 6 Figure 6.13: a. Schematic of the agarose sample structured with rice noodles with two orientations. b., c., and d. represent the wavelength maps calculated from a one-second acquisition using fingers taping as the source of mechanical noise. b. is the wavelength averaged over all the direction in the image plane, c. and d. are the wavelength maps obtained by considering the x and y-axes separately.

Figure 6 . 14 :

 614 Figure 6.14: Experimental data acquired on the forearm in-vivo using a digital holography set-up. A vibrating phone was used as the source of mechanical vibrations. Two different zones were imaged as displayed in a. The first one is closer to the wrist (b. and d.) and the other closer to the elbow (c. and e.). The estimated shear-wave wavelength maps are presented in b. and c. The corresponding average correlation function from the same acquisitions is presented in (d.) and (e.).

Figure 6 . 15 :

 615 Figure 6.15: Refocused shear-wave fields at t=0 (a., b. and c.) and t= τ (d., e. and f.) obtained from finite difference simulation images. Three configurations are displayed here: equipartition of sources and medium mechanical isotropy (a. and d.), equipartition of sources and medium mechanical anisotropy (b. and e.) and one-sided sources and medium mechanical isotropy (c. and f.).

  Fig. A.1 et comparées dans la Fig. A.2.

Figure A. 3 :

 3 Figure A.3: Exemples de différents types de contraintes sur un solide de forme pavé : forme initiale, étirement, compression, contrainte de cisaillement, torsion et flexion

  avait en effet comparé les mesures d'élasticité issues de tests mécaniques (0,1 à 4 Hz) et de l'élastographie ultrasonore sur le foie. Cette étude avait conclu qu'aucune comparaison directe et quantitative ne pouvait être faite entre les deux techniques de mesure en raison des différentes gammes de fréquences d'excitation mécanique. Des mesures de contrôle ont donc été effectuées avec le dispositif expérimental. Les approches d'élastographie impulsionelle harmonique et classique ont été utilisées ici pour la validation des algorithmes et les mesures de contrôle des échantillons fabriqués. Ce chapitre se concentre sur les mesures de contrôle effectuées à l'aide de l'élastographie harmonique et impulsionelle basée sur le temps de vol. Ces mesures de contrôle ont servi de référence pour les nouvelles approches présentées aux chapitres 5 et 6, basées sur la corrélation du bruit. Le montage expérimental d'holographie numérique hors-axe est décrit Figure A.4. Il a été testé sur des images de simulation, des échantillons de polymère et des tissus biologiques ex-vivo. Le protocole de simulation et celui des échantillons imitant les tissus biologiques sont également décrits. Les résultats des mesures de contrôle sont présentés dans la dernière section de ce chapitre

Figure A. 4 :

 4 Figure A.4: Schéma du montage d'holographie hors-axe : le faisceau objet (en rouge) est diffusé par l'échantillon avant d'interférer avec le faisceau de référence (en violet) (f L1 = 100mm, f L2 = 150mm)

Figure A. 5 :

 5 Figure A.5: A. Résultats expérimentaux de la corrélation de bruit (fréquence d'acquisition de 500 Hz) et du temps-de-vol (fréquence équivalente de 25 kHz). La courbe a été réalisée à partir de mesures sur des échantillons avec des concentrations massiques d'agarose croissantes (1.00, 1.25, 1.50, 1.75 et 2.00 %). Pour chaque concentration d'agarose, la réponse impulsionnelle du milieu a été tracée (a, b, c, d et e). B. Image de la rigidité à la surface d'un échantillon d'agarose à 2.00 % contenant une inclusion moins rigide à 1.00 % à partir de la corrélation de bruit.

Figure A. 7 :

 7 Figure A.7: Refocalisations d'un champs diffus dans un échantillon isotrope (a) et anisotrope (b)

Figure A. 8 :

 8 Figure A.8: Champs d'ondes de cisaillement refocalisés à t=0 (a., b. et c.) et t= τ (d., e. et f.) obtenus à partir d'images de simulation par différences finies. Trois configurations sont représentées ici : équipartition des sources et isotropie mécanique du milieu (a. et d.), équipartition des sources et anisotropie mécanique du milieu (b. et e.) et sources unilatérales et isotropie mécanique du milieu (c. et f.).

4. 2

 2 Processing for the 2D Fourier transform method to retrieve phase and amplitude images from each camera frame. . . . . . . . . . . . . . . . . . . 4.3 Photo (a) and corresponding schematic (b) of the compact off-axis digital holography set-up with an adaptable objective. The reference beam and the object beam are represented respectively in purple and red. The objective presented on the photo consists of two 100mm focal distance lenses with a diaphragm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 a. Frame captured by the CMOS captor with standard experimental condition (integration time of 200µs) without the laser source. The histograms presented in b., c., d. and e. correspond respectively to the top left, top right, bottom left and bottom right quadrant of the sensor. The total histogram of the captor is displayed in f. The average values (Mean) and standard deviation (StDev) are given for each curve. . . . . . . . . . 4.5 Histogram (a.) and standard deviation (b.) in time averaged over 1440 ×1440 pixels and 500 frames for an typical acquisition in ambient lighting without any object (200µs, 500 frames at 500Hz). . . . . . . . . . . . . . . 4.6 Noise measurements performed with the experimental set. The sensor was exposed to gradual intensities. The intensity fluctuations (in gray levels) were traced in function of I/ √ I (in gray levels -1/2 ). The curve was fitted with an affine line with R 2 = 0, 97. . . . . . . . . . . . . . . . . . . . . . . 4.7 Image auto correlation of a typical acquisition in an agarose sample to determine the speckle's size. (a.) Raw amplitude image, (b.) averaged auto-correlation, (c.) profile of the auto-correlation in blue fitted with the red curve to estimate FWHM. Algorithm by [7]. . . . . . . . . . . . . . . . 4.8 Experimental phase (a.) and phase difference (b.) images of a typical diffuse noise field acquired on agarose. The red square is a zoom on a 100 × 100 pixels zone to better see the speckle pattern. The histogram of the phase difference image is displayed in image (c.), where the central peak represents the useful signal while the secondary peaks around ∓2π are the artifacts due to speckle noise. Image (d.) is the same as image (b.) with a better suited color-map to compare with the filtered images (e.) and (f.) using respectively a median filtering and a threshold and Gaussian filtering. The threshold used is illustrated in (c.) by blue dotted lines. . . 4.9 Image of an USAF resolution target taken with the presented optical setup.(a.) is the full-field image and (b.) is a zoom centered on the resolution limit of the presented optical system. . . . . . . . . . . . . . . . . . . . . . 4.10 Shear-wave speed map calculated from the simulation movie of a sinusoidal impulse 4 ms and 1.3 ms long propagating through a star-shaped soft inclusion. The simulated shear-wave speeds are 3.0 m/s for the inclusion and 6 m/s for the background. In a. and b., only the time-of-flight calculation was used and in c. and d., a propagative filter was added prior to the time-of-flight processing. . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Acquisitions of harmonic excitation on agarose samples. On (a.), each curve represents a sample with a different agarose concentration. The measured wavelength is computed in function of the excitation frequency. At 1 kHz, images illustrating the different observed frequencies depending on agarose concentration are shown in (b.). These images show the same field of view and correspond to the gray vertical band in the curve. The images were acquired using a digital holography set-up. Gaussian filtering was used on the images. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Mean wavelengths given by the noise correlation method depending on the pulse speed measured by the control experiment. The red curve is the theoretical curve for a excitation of central frequency 1.1 kHz. These data correspond to agarose samples with different concentrations. Five different one-second acquisitions were made for each sample with a 200 µm exposure time. The mean PSF of the mechanical waves are illustrated for each sample (a, b, c, d and e respectively for 1%, 1.25%, 1.5%, 1.75%, and 2% agarose) in 600 × 600 pixels images . . . . . . . . . . . . . . . . . 133 5.4 Images from a diffuse field acquisition on a 1% inclusion in a 2% agarose sample with 1% titanium dioxide concentration. (a.) Description of the inclusion sample and mechanical noise sources. (b.) Amplitude image of the sample in ambient lighting. the inclusion zone is lined with red dotted lines. (c.) Phase displacement images on 4 consecutive frames. (d.)The wavelength map image is retrieved using the noise correlation method. Averaged images of the correlation functions for a zone in the softer background (e.) and for a zone in the inclusion (f.). . . . . . . . . . 134 5.5 Wavelength maps obtained on a pig liver and on chicken skin using noisecorrelation and digital holography. (a.) presents results on a porcine liver where the upper half was left raw as the other part was cooked. The average full cross-correlation is traced for each red square on the image for a 300 x 300 pixel window : (b.) and (c.) correspond respectively to the raw and cooked liver. (d.) images a circular nylon rod inserted in a chicken thigh. The nylon rod is highlighted with a red dotted line. . . . . 135 5.6 Images from a diffuse field acquisition on a 1% inclusion in a 2% agarose sample with 1% titanium dioxide concentration. (a.) Description of the inclusion sample and mechanical noise source. (b.) The wavelength map image is retrieved using the noise correlation method. (c.) The frequency content of the mechanical excitation was computed from the acquisition to extract the central noise-frequency. . . . . . . . . . . . . . . . . . . . . 136 6.1 Schematic of the time-of-flight and NCi algorithm for retrieving shearwave speed estimation from full-field acquisitions. . . . . . . . . . . . . . . 145 6.2 Frames from simulated propagation movies (500×500 pixels of 100 µm² with 1000 frames at 50 kHz) displayed in Video 1 (MP4, 0.20 MB) and Video 2 (MP4, 0.29 MB). A soft star-shaped inclusion is simulated with 3 m/s shear-wave speed for the inclusion and 6 m/s for the background . (a), (b), (c) and (d) correspond to frames 350, 450, 550 and 650 of Video 1, for which the shear-wave source function is a 4 ms-long half-sine. (e), (f), (g) and (h) correspond to frames 350, 450, 550 and 650 of Video 2, where the source function is 4 ms-long a noise signal (400-800 Hz). . . . . 146 6.3 Speed maps obtained from simulated propagation movies (images of 500×500 pixels of 100 µm², 1000 frames at 50 kHz). The medium of propagation was a star-shaped inclusion in a softer background, illustrated by the map (a.). Using time-of flight ((b.) and (c.)) and NCi algorithm ((d.) to (g.)). Results with different profiles of one-dimensional shear-wave sources are presented: half-sine functions of 4 and 1.33 milliseconds, a chirp function and a noise signal. Profiles of the excitation sources are represented above each speed map: half sine functions ((b.

  The standard deviation is presented and the gray line represents the x=y function. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Curves representing the average speed of Rayleigh waves measured at the surface of scattering agarose samples of gradually increasing stiffness (from 0.5 to 2.0% agarose mass concentration) using the NCi derivation-based approach. Time-of-flight measures were used as a reference and presented on the x-axis of each graph. (a) and (b) were obtained respectively using four 1 ms-long half-sine excitation function and using a 2 ms-long chirp excitation function, from 500 to 2 kHz. The measures are averaged over six acquisitions for each sample. The standard deviation is also presented. The gray line represents the x=y function. . . . . . . . . . . . . . . . . . . 6.7 Images obtained from a 2.0% agarose sample with a 1.0% star-shaped agarose inclusion. (a.) describes the sample and mechanical source. (b.) is an amplitude image of the sample. The inclusion cannot be distinguished by its optical properties in this image and its shape and position is highlighted by a red dotted line. (c.) shows the raw phase difference images obtained with the digital holography set-up in stroboscopic mode (25 kHz equivalent frame rate). Frames 108, 179, 250 and 321 of the 500 are displayed. (d.) is the speed-map reconstruction from the acquired propagation movie achieved with the NCi algorithm. . . . . . . . . . . . . 6.8 The spatiotemporal representations of a 1 ms-long pulse propagating in a cooked and raw pork liver sample are presented in (a) and (b) respectively.

  the delay between the light arriving at M from S passing through O 1 and O 2 . The amplitude contributions of the light from each secondary source O 1 and O 2 are noted s 1,ω (M ) and s 2,ω (M ) respectively.

  3, no. 6, p. 199, Sep. 1998. doi: 10.1364/oe.3. 000199.

Lors de cette thèse, plusieurs approches pour l'élastographie ont été proposées. Tout d'abord, la preuve de concept de l'élastographie par corrélation de bruit sur un montage d'holographie numérique a été effectuée. Pour cela le montage expérimental a été mis en place pour l'élastographie par corrélation de bruit mais aussi pour de l'élastographie impulsionnelle classique, utilisée pour fournir une mesure de référence. Ensuite la méthode de corrélation de bruit a également été adaptée aux sources impulsionnelles. Enfin, le cas des milieux présentant une anisotropie mécanique a été examiné. Les résultats sont présentés en quatre parties : le montage expérimental, la preuve de concept de la méthode passive, son application à des sources spatialement cohérentes et enfin son utilisation dans les milieux anisotropes.

Remerciements

Je tiens à remercier sincèrement les personnes qui ont fait partie de l'expérience qu'a été cette thèse, que ce soit en acceptant de faire partie du jury, en l'encadrant ou en y collaborant. J'espère que les personnes qui y ont principalment contribué se retrouveront dans cette grille de mots croisés aux côtés de quelques mots qui résument ces dernières années de travail.

naturellement présent in-vivo. La corrélation de bruit permettrait donc de s'affranchir de toute source externe d'onde mécanique, rendant notamment l'élastographie adaptable plus facilement au bloc opératoire. L'application d'une méthode de corrélation à un système d'acquisition ligne-par-ligne comme l'OCT est cependant restreinte à une seule dimension spatiale. L'objectif de cette thèse est d'adapter cette méthode passive à un montage interférométrique plein-champ afin de fournir des images de rigidité quantitatives et sans-contact au bloc opératoire. Le système développé présente également l'avantage de pouvoir fournir des mesures sur des tissus biologiques présentant des propriétés mécaniques anisotropes.

Chapitre 5 -Elastographie par corrélation de bruit

Les approches d'élastographie impulsionnelle pour les tissus biologiques ont d'abord été introduites sur la base de la propagation d'une impulsion mécanique unique. Pour l'optique, la limitation majeure de cette approche classique est la génération d'une onde mécanique contrôlée et son suivi sur des champs de vue plus petits. En effet, l'élastographie impulsionnelle classique nécessite le suivi d'ondes de cisaillement se propageant à 1-10 m/s dans les tissus mous. Ceci implique l'utilisation d'une caméra rapide ou d'une synchronisation par verrouillage de phase entre la caméra et la source mécanique comme montré dans le chapitre précédent de cette thèse. Une autre limitation de l'approche classique de l'élastographie impulsionnelle est qu'elle est basée sur la propagation balistique des ondes mécaniques. Comme nous l'avons montré dans le chapitre précédent sur des images de simulation et des échantillons ex-vivo, la présence de la diffusion et des interférences provenant des inhomogénéités ou des limites des tissus peut entraver le suivi des ondes de cisaillement. De même, la présence d'autres ondes mécaniques provenant de vibrations naturelles (battements cardiaques et autres activités musculaires) peut limiter l'application clinique de cette méthode. De nouvelles méthodes d'élastographie passive ont donc été développées pour travailler avec des champs élastiques complexes, qu'ils soient générés de l'extérieur [5], [6] ou présents naturellement in-vivo. [7], [8].

Les approches utilisant la corrélation du bruit en sismologie [9] et le retournement temporel (TR) [10] ont été récemment adaptées à l'élastographie des tissus biologiques mous [11]. [12] pour surmonter la principale contrainte des méthodes classiques en exploitant les champs d'ondes de cisaillement diffus. Ces approches sont basées sur une refocalisation numérique du champ d'ondes de cisaillement en utilisant des algorithmes de corrélation du bruit. À chaque pixel, le champ refocalisé donne accès à la longueur d'onde locale des ondes de cisaillement qui est directement liée à la rigidité locale.

Cette méthode a d'abord été développée pour les systèmes à ultrasons mais en 2016, Nguyen et al. [13] ont adapté la corrélation du bruit à l'élastographie optique. La preuve de concept de l'approche de corrélation du bruit sur les images optiques a été faite sur des images OCT. En raison de la nature des images OCT, la corrélation a été calculée ligne par ligne. Le montage plein champ présenté dans cette thèse offre une adaptation de cette méthode de corrélation de bruit aux images plein champ où la corrélation peut être faite sur deux dimensions. La méthode d'élastographie présentée ici est basée sur un champ de bruit mécanique pour imager l'élasticité, ce qui élimine le besoin d'un générateur d'impulsions contrôlé.

L'approche de corrélation du bruit est présentée dans ce chapitre, ainsi que les résultats de sa combinaison avec le dispositif d'holographie numérique plein champ présenté dans le chapitre précédent. La preuve du concept est présentée ainsi que les résultats sur les échantillons d'inclusion et les tissus ex-vivo. [14].

La Figure A.5 A. présente la preuve de concept de l'élastographie par corrélation de bruit. Des mesures de rigidité sont évaluées par élastographie impulsionnelle classique (axe des abscisses) et par élastographie par corrélation de bruit (axe des ordonnées). La courbe rouge correspond à l'équivalence des deux méthodes. On peut voir sur cette figure une forte corrélation des deux types de mesures. La Figure A.5 B. présente quand à elle un exemple d'une carte de rigidité obtenue sur un échantillon aux propriétés mécaniques Chapitre 6.1 -Méthode alternative de traitement des images de propagation issues de sources spatialement cohérentes Les approches de l'élastographie présentées précédemment ont décrit deux méthodes principales pour réaliser l'élastographie impulsionnelle : l'une utilisant des sources d'excitation mécanique unidimensionnelles pour les mesures de contrôle (chapitre 4) et l'autre basée sur les champs de bruit (chapitre 5). Ces deux méthodes n'ont été envisagées que pour des milieux isotropes. Cette section vise à combler le fossé entre ces deux approches en examinant l'influence de l'anisotropie de la source et du milieu dans les reconstructions de cartes de vitesse.

En effet, contrairement aux systèmes à ultrasons ou OCT utilisés précédemment pour l'élastographie, l'approche plein champ de l'holographie numérique permet des corrélations bidimensionnelles. Cela permet d'obtenir des informations directionnelles sur la réponse de l'échantillon. À partir de ces données, il est possible d'extraire des informations sur la source ou l'anisotropie du milieu.

La première partie de ce chapitre 7 se concentre sur l'utilisation d'un algorithme basé sur la corrélation du bruit pour les sources non bruyantes telles que les sources d'ondes de cisaillement spatialement cohérentes utilisées en élastographie impulsionelle classique. En exploitant la corrélation bidimensionnelle, la réponse du milieu à l'impulsion peut être traitée pour extraire la direction de propagation locale et la vitesse de l'onde de cisaillement. Les résultats obtenus sur des images plein champ se révèlent plus robustes aux phénomènes de réflexion et de réfraction que l'approche par temps de vol [15].

Des résultats sur des images de simulation sont présentés

4.13

Curve representing the measured shear-wave speeds in homogeneous samples with varying agarose mass concentration. Each point was computed on 5 different acquisitions. The mean square error is displayed using vertical error bars. 1-ms-long mechanical impulses were generated in the samples. To illustrate these results, a. to g. are examples of spatiotemporal profiles from experimental propagation movies acquired in homogeneous agarose samples with gradually varying mass concentration (0.5, 0.75, 1.0, 1.25, 1.50, 1.75 and 2.0 % respectively). The x-axis represents the consecutive frames of the movie and the y-axis is a line on the image along the propagation axis. Directional and Gaussian filtering was used on the images. All the images were acquired using a digital holography set-up and a stroboscopic approach to reach 25 kHz equivalent frame rate. . . . .