Shengfu Par

Li

Didier Vray

CREATIS Professeur

Silva Da

Anabela

Kouame Denis

R Gourdon

M Jean

Marc Lancelin

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

M Obadia

Contribution to quantitative photoacoustic reconstruction: forward models and inversion schemes Délivrée par

Keywords: Photoacoustic imaging, biomedical imaging, optical fluence, analytical model, numerical model, inverse problem, multigrid methods, quantitative reconstruction, tomography Imagerie photoacoustique, imagerie biomédicale, fluence optique, modèle analytique, modèle numérique, problème inverse, méthodes multigrilles, reconstruction quantitative, tomographie

Photoacoustic imaging (PAI) of biological tissues tries to combine the advantages of optical and acoustical imaging. The main endogenous contrast for PAI is derived from blood vessels due to the strong absorption of hemoglobin compared to the background tissues. Furthermore, blood vessels are roughly cylindrical and hemoglobin concentration can be assumed to be uniform inside the vessel. Therefore, the blood vessels can be considered as "cylindrical inhomogeneities". As a first contribution, we have developed in this thesis an analytical model of optical fluence for multiple parallel cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium.

Analytical models only exist for simple cases. To deal with more complex situations like biological tissues, numerical methods are required. The second contribution of this thesis is to develop a multigrid solver of optical diffusion equation and therefore to propose an efficient numerical method to resolve the optical fluence.

Finally, our third contribution is concerned with quantitative PA tomography (QPAT) reconstruction. Based on the efficient models presented in the first and second contributions, we have proposed an analytic-based reconstruction method for simple cases and a multigridbased inversion scheme for more realistic cases. The advantages of multigrid-based inversion scheme are shown in both computation and convergence speed. An experimental validation is presented in the last chapter of this thesis, proving the validity and analyzing the performances of the developed methods.

Notations

and (e) are 2.1 cm and 2. 3cm, depths of the c2 center in (a) and (e) are 1.5cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (e)): Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1 and c2. Third row ((c) and (g)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different radius (radius of c1 in (a) and (e) is 0.5cm and 0.7cm, respectively, radius of c2 in (a) and (e) is 0.3cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) and (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row ((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm . 49 Fig. 3.7 Comparison between the analytical model and NIRFAST for two inhomogeneities with different absorption. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different absorption coefficient (µ a of c1 in (a) and (e) is 0.4cm -1 and 0.8 cm -1 , respectively, µ a of c2 in (a) and (e) is 0.8 cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row ((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] × [-1.2cm, 1.2cm . 3.9 Comparison between the analytical model and NIRFAST. (a)Horizontal cross section of the first numerical phantom with five cylindrical inhomogeneities (µ a =0.8cm -1 and µ' s =10cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c1 and c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the position of the inhomogeneity c5. (d) MRE corresponding to the area [-1.2cm, 1.2cm] × [-1.2cm, 1.2cm]. .. 53 Fig. 3.10 Comparison between the analytical model and NIRFAST. (a) Horizontal cross section of the second numerical phantom with three cylindrical inhomogeneities (see Table 3.1 for optical properties) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). ( b Fig. 3.12 Multiple updates of x (0) . The grid levels are noted by (0) to (3). Red broken line represents multigrid algorithm recursively moves forth and back. The dotted lines represents grid with different spacing. The position of each update of x (0) is noted by black points along the green line with an arrow. The number of updates is noted by 1,2,…. ............................... 57 Fig. 3.13 Illustrations of grid in finite difference. red circles represent the positions of space points. (x i ,y j ) represents coordinate values of an arbitrary point in space................................ 59 Fig. 3.14 Validation of home-developed multigrid solver of diffusion equation, (A) geometry of numerical phantom consisting of three absorber with absorption coefficient of 0.3cm -1 (larger absorber with radius of 0.5cm),0.2cm -1 (middle absorber with radius of 0.4cm) and 0.15cm -1 (small absorber with radius of 0.3cm) embedded in a homogeneous medium with absorption coefficient of 0.1cm -1 . µ s =10cm -1 . (B) Optical fluence distribution calculated with multigrid method. Profile through the optical fluence distribution calculated with 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm (absorber c2). (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm (absorber c3). (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm (absorber c1). . 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm (absorber c2). (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm (absorber c3). (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm (absorber c1). .

Introduction Motivation

Photoacoustic imaging (PAI) uniquely combines the merits of optical imaging and acoustical imaging. PAI is based on PA effect which refers to the generation of acoustic waves from an object being illuminated by pulsed optical waves. PAI has the advantages of pure optical imaging by using non-ionizing and non-invasive sources, providing sensitive absorption contrast related to molecular constitution in soft tissues [START_REF] Cox | Quantitative spectroscopic photoacoustic imaging: a review[END_REF][START_REF] Favazza | In vivo functional photoacoustic microscopy of cutaneous microvasculature in human skin[END_REF][START_REF] Niederhauser | Combined ultrasound and optoacoustic system for real-time high-contrast vascular imaging in vivo[END_REF][START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF]. Moreover, scattering limits the penetration of ballistic photons to a few mm in biological tissues while PAI takes advantage of the penetration of scattered photons to detect PA signals from deeper region [START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF]. Therefore, PAI has shown fast developments during the past decade in both theoretical methods and innovative applications for biomedicine. This thesis is trying to bring a personal contribution to this field.

PAI applications in human and animal imaging

In the last decade, biomedical applications in human and animal PAI have been investigated widely. The following paragraphs focus on the more significant works.

PAI of small animal models: pre-clinic application

The small animals such as mice have been widely used in research platform to study structural and functional imaging with PAI system. The penetration depth is the principal limit for PAI applications in biological tissues. Thus, the dimension of the small animals, a few centimeters, provides the possibility to image inner structure. Fig. 0.1 shows an example of whole-body three-dimensional PAI for small animals [START_REF] Brecht | Whole-body three-dimensional optoacoustic tomography system for small animals[END_REF].

As well as providing structure information, PAI has been used to extract the functional changes in biological tissues. The brain of small animal models has been widely studied. The attractiveness of small animal brain lies in the sufficiently thin skull compared with larger animals such as indeed humans. Hard tissue such as bone strongly scatters ultrasound waves, thus a thin skull means that it is potential to obtain high image quality of the mouse or rat brain. PAI can be used to image the vasculature and to extract quantitative functional information by measuring of blood oxygen saturation and hemoglobin concentration. Fig. 0.2 shows an example of PAI of small animal brain [START_REF] Yang | Photoacoustic tomography of small animal brain with a[END_REF]. contrast on a PA image. As usually, tumor development is associated with neovascularization around the tumor, PAI can be applied to aid cancer detection and diagnosis, some researches of PAI of breast [START_REF] Ermilov | Laser optoacoustic imaging system for detection of breast cancer[END_REF][START_REF] Lee | Optical mammography: Diffuse optical imaging of breast cancer[END_REF] and prostate cancer detection [START_REF] Yaseen | Optoacoustic imaging of the prostate: development toward image-guided biopsy[END_REF] can be found in published works.

Skin

PAI of the skin is relatively easy to implement because the penetration depth requirement is reduced to a few mm achievable under considering safety standard. PAI instruments have been developed with the required ability for imaging the vasculature and other features of skin [START_REF] Favazza | In vivo functional photoacoustic microscopy of cutaneous microvasculature in human skin[END_REF][START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF]Zhang, Laufer, Pedley, & Beard, 2009).

Image formation in PAI

In PAI, the recorded signal represents an acoustical pressure wave that depends on the absorption of optical illumination. By detection of the pressure wave with acoustic transducers on the sample surface, a 2-D or even 3-D image related to the absorption of optical radiation can be produced by using an appropriate image formation method.

Two major forms of PAI systems have been implemented [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF]. One is based on a scanning-focused acoustic transducer, so-called PA microscopy (PAM). Another form is called PA tomography (PAT), in which an array of unfocused acoustic transducers is used to record PA waves and an image can be obtained by applying appropriate image formation method.

In PAM, the time-resolved acoustic waves recorded at each scanning location are converted into a 1-D depth-resolved image based on the acoustic velocity of sample. Then, raster scanning in a plane with a certain step size produces a 3-D image.

In PAT, unlike PAM, a reconstruction algorithm is required, which consists of producing the optical properties of sample from the boundary measurements of the acoustic waves. Typically, two steps are involved: (1) the first step concerns the reconstruction of the deposited optical energy, which is called conventional PAT reconstruction in literature; (2) the second step, so-called quantitative PAT (QPAT) consists of reconstructing optical parameter maps, particularly the absorption coefficient, from the deposited optical energy. Reconstruction methods will be more detailed in chapter 2 of this document.

Objectives

As a general objective, we aim at reconstructing the optical properties of the illuminated sample. First, we will develop a forward analytical model of light propagation for simple geometry and also a multigrid numerical model for more complex 3D configuration. Secondly, we will propose an efficient inversion scheme for QPAT reconstruction based on the forward models developed previously. This inversion scheme will be validated with experimental data.

Plan of this thesis

This thesis includes five chapters. The first two chapters will give basic principles, from physical principle of PA effect to PAI systems. We will mainly focus on the state of the art about existing optical models and reconstruction methods. Chapter 1 presents the principle of the PA effect in biological tissues, firstly illustrating the concept of optical absorbers, and then presenting the basic models to describe the propagation of light, including Monte Carlo, finite element method (FEM) and analytical model. From the aforementioned concepts of PAI, it can be known that a complete PAI system consists of optical illumination, acoustic detection and image formation. In chapter 2, the three elements are to be discussed, including the parameters of illumination source, ultrasound detectors and imaging formation method (reconstruction algorithm). In chapter 3, we will present our contribution of optical forward models, including an analytical model of optical fluence for multiple cylinders embedded into an otherwise homogeneous medium, and also an efficient multigrid numerical model of (2) the number of parameters to be recovered; (3) optimization method. Chapter 5 presents an experimental validation of our QPAT reconstruction algorithms. The experiment has been realized in collaboration with Dr. Zhen Yuan, University of Macau. We firstly present the experimental system and the design of a realistic phantom. The discussion focuses on the comparison between the phantom parameters and recovered results with our algorithms.

The manuscript is concluded with the perspectives of this work.

Photoacoustic effect

The photoacoustic (PA) effect, also called the optoacoustic effect, concerns the generation of acoustic waves from an object being illuminated by pulsed or modulated laser source [START_REF] Razansky | Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion[END_REF][START_REF] Xi | Evaluation of breast tumor margins in vivo with intraoperative photoacoustic imaging[END_REF][START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF]. This phenomenon was firstly discovered by Alexander Graham Bell in 1880 [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF]. During his investigations of "photophone" (see in Fig. 1.1), he noticed that a sound wave could be produced directly by rapidly interrupted beam of sunlight focused on a solid sample. He observed that the produced acoustic signals were dependent on the composition of the sample. He conjectured that these signals were caused by absorption of the incident light. He explained the phenomenon as following: the light pulse was absorbed by the sample of matter, then converted to equivalents energy, and then partially radiated as heat (generation of wave). Consequently a pressure wave can be detected by an acoustic sensor. (http://photonics.cusat.edu/). The optical properties of biological tissue consist of the absorption coefficient, µ a (cm -1 ), the scattering coefficient µ s (cm -1 ) , the scattering function p(θ) (sr -1 ) where θ is the deflection angle, and the refractive index, n' (Wang & Jacques, 1992). The strong scattering of biological tissue averages the deflection angle and then an average parameter, g = <cos θ>, called the anisotropy of scatter, is used to describe tissue scattering in terms of the relative forward versus backward direction of scatter [START_REF] Boas | The Migration of Diffuse Photon Density Waves through Highly Scattering Media[END_REF]. Reduced scattering coefficient ' s  is often used to describe the scattering of biological tissues, which is defined as, '

(1 ) ss g

 

. Optical scattering can be modeled as scattering by particles. A mixture of spheres of different sizes is used to describe the optical scattering behavior of biological tissue (Jacques, 2013).

Absorber in biological tissues

In PAI of biological tissues, a fraction of incident light energy will be absorbed along the optical path length within the medium. The absorption coefficient is generally noted by µ a , with a unit of cm -1 . Absorption coefficient can be measured by the transmission through a nonscattering medium containing chromophores as  [START_REF] Reif | Analytical model of light reflectance for extraction of the optical properties in small volumes of turbid media[END_REF][START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF], where C is the concentration (mol•L -1 ) and ε is the extinction coefficient (L•mol -1 •cm -1 ) of the chromophores. Absorption coefficient of biological tissues is equal to the sum of contributions from all absorbing chromophores within the tissue,   ln 10

a i i i C    .
Studies of tissue optical properties usually use values of a tissue's average absorption coefficient, since the molecular composition of the biological tissues is not well known (Jacques, 2013).

To illustrate the concept of absorption coefficient of biological tissues, we take the blood as an example. The hemoglobin molecule provides the main endogenous chromophore, including oxyhemoglobin (HbO 2 ) and deoxygenated hemoglobin (Hb). Thus the absorption coefficient of blood is equal to the sum of contributions from Hb and HbO 2 , as following equation [START_REF] Yang | Photoacoustic tomography of small animal brain with a[END_REF], Another way to calculate the optical absorption coefficient of biological tissues is based on the volume fraction of a tissue f v.i and the absorption coefficient of that pure component μ a.i , with a unit of cm -1 . Using this approach, the absorption coefficient can be expressed as,

        2 2 , ( ) ( )
BS B S W F M C C                 
Therefore, it is important to consider the absorption spectra of all kinds of absorbing chromophores to obtain the total absorption coefficient. Figure 1.2 shows absorption spectrum of fully oxygenated blood, fully deoxygenated blood, fat, water and aorta. These curves in Fig. 1.2 are based on data from published work. Note that reliable data beyond 1000 nm wavelength are difficult to find in the literature.

Multi-Physics Coupling: optics and acoustics

When incident light propagates through biological tissues, part of energy is absorbed and produces an absorbed energy distribution, noted by H(r,t) in published papers [START_REF] Xu | Photoacoustic tomography of water in phantoms and tissue[END_REF]. A fraction of absorbed energy converts into heat and produces a rise of temperature, which then produces a pressure distribution, so called initial pressure, noted by p 0 in papers. The pressure distribution p(r,t) produced by the heat source H(r,t) obeys to the following equation [START_REF] Xu | Photoacoustic tomography of water in phantoms and tissue[END_REF]: The initial pressure of the absorber at location r after absorbing laser energy can be expressed as [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF] From the aforementioned concept, it can be seen that a model of optical fluence and a method to resolve Eq. (1.2) are required to completely determine PA wave. In this chapter, 1.2 will present the model of optical fluence that is used to describe optical fluence and 1.3 will present the method to resolve Eq. (1.2).

      2 2

Optical propagation model in biological tissues

As mentioned before, in PAI of biological tissues, PA waves are produced by absorbed energy which is the product of optical fluence and absorption coefficient. In this section, the existing models that are used to describe optical fluence distribution are presented. The optical fluence has been obtained by Monte Carlo simulation or by resolving the differential equations that describe the optical flux distribution.

Statistical method: Monte Carlo

Once a single photon packet enters biological tissues, it may be scattered, absorbed, propagated undisturbed, internally reflected or transmitted out of the tissue. The typical trajectory of a single photon packet in a homogeneous medium is given in Fig. 1.3. (Prahl, Keijzer, Jacques, & Welch, 1989). Upon launching, the photon is moved a distance ∆s, in which the photon may be scattered, absorbed, propagated undisturbed, internally reflected, or transmitted out of the tissue. The photon is repeatedly moved until it either escapes from or its energy is totally absorbed by the tissue. If the photon escapes from the tissue, the reflection or transmission of the photon occurs due to the change of properties of medium at the interface. If the photon is absorbed, the position of the absorption is recorded. This process is repeated until achieving the desired precision (determining by the number of launched photons). For a tissue with the specified optical properties, the recorded reflection, transmission, and absorption distribution will approach true values as the number of launched photons approaches infinity. ∆s must be less than the average mean free path length of a photon in the medium. The mean free path length is defined as the reciprocal of the total attenuation coefficient,

1 as s  D 
The probability that the photon may be reflected by interface is determined by the Fresnel reflection coefficient. Figure 1.5 represents optical fluence distribution produced by a MCbased package, which was developed by Wang, etc. (Wang & Jacques, 1992). In Fig. 1.5, the medium is a homogeneous turbid slab,

Differential equation (1): Radiative transfer equation (RTE)

The general model of light transport that is expressed as differential equation is radiative transfer equation (RTE), which has been studied and resolved with numerical methods [START_REF] Tarvainen | Gauss-Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation[END_REF]. RTE, also referring to as Boltzmann equation, conveys that a light beam loses energy by divergence and extinction as passing through the medium [START_REF] Hielscher | Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues[END_REF]. In RTE, coherence, polarization and nonlinearity are ignored. Because biological tissues are high scattering medium and photons are multiply scattered in tissues, thereby phase and polarization are quickly randomized, and play little role in light energy transport (Chapter 2, Boas, 1996). In RTE, optical properties such as refractive index n', absorption coefficient μ a , scattering coefficient μ s and scattering anisotropy g are assumed as time-invariant but spatial-variant. Generally, only elastic scattering (changing the direction of the photon, not the frequency) is considered.

RTE describes the behavior of the radiance   ,, r t s  with units of Wm -2 sr -1 which represents the energy traveling along the direction normal to ŝ , per unit solid angle, per unit time, and per unit area. The final equation for the time-dependent case is [START_REF] Haskell | Boundary conditions for the diffusion equation in radiative transfer[END_REF][START_REF] Lee | Scattering of evanescent wave by multiple parallel infinite cylinders near a surface[END_REF],

            '' 4 1 ˆˆˆˆˆˆˆˆ, , , , , , , , , , , 4 t
                (1.6)
Where v is the velocity of light inside the medium determined by v=c/n' , c is the speed of light in vacuum and n' is the refractive index of the medium. µ t =µ a +µ s is the total extinction coefficient where µ a and µ s are the scattering and absorption coefficients inside the medium, respectively.   ,, S r t s refers to the spatial and angular distribution of the source with units of 31 Wm sr  ,   ˆ,' p s s refers to the scattering function or the phase function which defines the probability of a photon traveling in the direction ' ŝ to be scattered into direction ŝ .

The two first terms on the left-hand side of RTE represents the temporal and spatial variations of radiance along the direction ŝ . The third term represents the radiance loss of photons due to extinction. The first term on the right-hand side of RTE represents the photons being scattered from all directions ' ŝ into direction ŝ (Chapter 2, [START_REF] Boas | The Migration of Diffuse Photon Density Waves through Highly Scattering Media[END_REF].   ,, S r t s is the source term. Therefore, RTE requires that the total space and time variation of the radiance along certain direction in an elementary volume be equal to the variation of radiance due to scattering and absorption inside the volume.

Analytic solutions of the RTE typically exist only for some special geometry, such as planar geometries illumination with a plane wave. It is resolved by numerical methods and needs large amounts of computational power because both the angular space and spatial space need to be discretized.

Differential equation (2): Diffusion equation

The diffusion equation is expressed as the relation between the optical fluence rate and optical parameters. It is deduced from RTE based on some reasonable approximations. Here, we briefly review the process from RTE to diffusion equation. The optical fluence rate is given by [START_REF] Haskell | Boundary conditions for the diffusion equation in radiative transfer[END_REF][START_REF] Hielscher | Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues[END_REF],

  ˆ( , )
, , r t ds r s t

  

Photon flux is given by,

  ˆˆ( , )
, , F r t ds r s t s   

Fluence and flux have same units, Wm -2 .

From RTE to diffusion equation, the variables in RTE need to be expanded with respect to spherical harmonics (Chapter 2, [START_REF] Boas | The Migration of Diffuse Photon Density Waves through Highly Scattering Media[END_REF][START_REF] Cox | Quantitative spectroscopic photoacoustic imaging: a review[END_REF], usually noted by Y l,m (l is the order, m is the degree m, l is a positive integer and m meets -l≤m≤l ) (Jackson, 1995)   and the scattering is not highly anisotropic (e.g. g<0.9, here g is anisotropy factor which has been defined before).

(2) Source term refers to isotropic source. Fortunately, the two assumptions above are suitable for biological tissues, since most biological tissues meet ' a s

  , and laser will lose the direction information after one transport mean free path in biological tissues due to the fact of highly scattering, although laser beam has high coherence in free space [START_REF] Patterson | Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties[END_REF]Farrel,1992).

Given these assumptions, we arrive at the diffusion equation of   , rt  for homogeneous turbid medium [START_REF] Li | Fluorescent diffuse photon density waves in homogeneous and heterogeneous turbid media: analytic solutions and applications[END_REF],

Shengfu LI Photoacoustic Imaging INSA Lyon 15         2 0 , , , , a rt D r t r t S r t vt           (1.7)
Note that in the frequency-domain the diffusion equation can be rewritten as the Helmholtz equation [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF],

      0 22 AC AC Sr kr D     
Where the wavenumber is complex, i.e., In PAI, the acoustic propagation occurs on a timescale several orders of magnitude longer than the heat deposition, therefore the time-integrated absorbed power density (i.e., the absorbed energy density) is the quantity of interest. The time-independent diffusion equation has the form [START_REF] Cox | Quantitative spectroscopic photoacoustic imaging: a review[END_REF]

, 22 () ( ) ( ) Sr r k r D       and a k i D   Where 2 1 i  , ' s  is the reduced scattering coefficient, a  is the absorption coefficient,   1 ' 3 as D      
 is the diffusion coefficient and () Sr is the source term.

Methods resolving the diffusion equation

The optical fluence distribution can be obtained by solving the diffusion equation. The solutions of a diffusion equation can be categorized as numerical solutions and analytical solutions. In PAI applications, numerical solutions are usually adopted due to their simplicity to implement. Meanwhile, we will also focus on situations where analytical solutions have a role to play.

Analytical model

Analytical models exist only for some simple geometry, like spherical wave incidents on infinite homogeneous medium or a spherical (cylindrical) inhomogeneity embedded an otherwise homogeneous medium.

Homogeneous medium

In an infinite homogeneous medium the solution to diffusion equation for a point source at s r is given by other studies [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF][START_REF] Patterson | Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties[END_REF]. Optical fluence at an arbitrary position ( , )

     then,     2' 2' '2 0 inc inc dr kr dr     (1.9)
The general solution of Eq. (1.9) is,

Inhomogeneous medium

In presence of a spherical (cylindrical) heterogeneity, () r  is found by constructing general solutions to diffusion equation outside and inside spherical (cylindrical) heterogeneity. Then, the appropriate boundary conditions are applied to resolve the unknown coefficients in general solutions [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF][START_REF] Li | Fluorescent diffuse photon density waves in homogeneous and heterogeneous turbid media: analytic solutions and applications[END_REF]. Generally, the problem is analyzed in spherical (cylindrical) coordinates whose origin coincides with the center of the spherical (cylindrical) heterogeneity. Boundary conditions require that: (1) the flux normal to the boundary of the heterogeneity must be continuous;

(2) the optical fluence must be continuous across the boundary of the heterogeneity. The general solution outside the heterogeneity is in the form of a superposition of incident (inc) and scattered (scatt) waves i.e.,

Numerical model

In the research about PAI, finite element method (FEM) was widely adopted [START_REF] Laufer | Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration[END_REF][START_REF] Tarvainen | Gauss-Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation[END_REF]Z. G. Wang, Sun, Fajardo, & Wang, 2009;[START_REF] Yao | Finite-element-based photoacoustic tomography in time domain[END_REF]. Here, FEM method is simply reviewed from the point of implementing standpoint. In FEM scheme, the domain of the problem, Ω, is divided into a collection of subdomains, e.g., mesh (2D) or tetrahedron (3D). The fluence at a given point, Φ(r), is approximated by choosing a set of basis functions that are simple polynomials and then the exact solution is approximated by a combination of those polynomial functions,

    1 V h ii r U u r    
. Ω h is a finite dimensional subspace spanned by a set of basis functions,  

( ); 1,2, i u r i V  .
Then, diffusion equation is converted to a discrete equation of the form, KU =F. Here, K is a matrix of size V×V, F is a vector of size V×1 and U is a V×1 vector, V is the number of nodes in the mesh. Therefore, to resolve the diffusion equation, three steps are necessary: (1) The domain of the problem is divided into a collection of subdomains. ( 2 The K,U and F have entries given by (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF][START_REF] Bhowmick | Numerical Modeling of Near-Infrared Light Propagation in Tissue[END_REF],

      , n i j i j a i j K D r u u r u u d r            n ii F u r S r d r   
Here, a simple example of optical fluence produced by NIRFAST (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF] with steps above is given in Fig. 1.7, in which a homogeneous turbid medium is used with



. The light source is placed at the point (x=1.9 cm, y=0 cm).

Discussion

MC-based model of photon migration in biological tissues has been proven to be an efficient method. Due to the generality and few hypothetical conditions, this method has been used as the gold-standard in this area [START_REF] Fang | Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units[END_REF]. However, this kind of model is time-consuming due to its low computational efficiency. For a typical domain size in the human head, it takes hours of computational time, therefore this method is low efficiency in comparison to a few seconds for solving diffusion equation or minutes for solving the RTE. Although several optimization techniques have been developed to improve the typical MC- based model of photon migration [START_REF] Fang | Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units[END_REF][START_REF] Fang | Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates[END_REF], it remains far from the efficiency needed for day-to-day use .

There are two partial differential equations to describe optical propagation in turbid tissues, RTE and the diffusion equation. RTE is usually solved by numerical method due to the fact that analytical solutions exist only for a few simple geometries. Also, it needs spatial and angular discretization to solve the RTE by numerical method. Therefore, fast way of accurately calculating the light distribution in large heterogeneous scattering media is based on the diffusion equation. And it will be used in this thesis.

Methods resolving PA wave equation

PA wave equation, like diffusion equation, can be resolved by two ways: analytical and numerical methods.

Analytical solution of PA waves emitted by spherical source

Analytical model of PA wave from a spherical source with uniform initial pressure distribution can be found in published papers [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF], it is expressed as,

     0 , 2 ss A U a R v t R v t p R t R     (1.12)
Here, p(R,t) represents acoustic pressure, at R and time t, emitted from this spherical source with a radius of a, R is the distance between the center of sphere and the position transducer. U(x) is the step function, with U(x) = 1 when x>0, and U(x) = 0 when x<0. A 0 is the amplitude of initial pressure the sphere. Some related parameters can be found in Fig. 1.8, in which PA signal is recorded by a transducer positioned at location R from the center of spherical source. assumes an acoustic speed of 1500 m•s -1 and a unit initial pressure. The pressure wave generally includes both positive and negative values even A 0 is uniform and constant inside the sphere. This signal is so called 'N-shape' signal.

A simple and intuitive interpretation of this signal shape is as follows: the initial pressure produces two waves [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF][START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF]. One propagates outward as a diverging spherical compression wave, producing the positive amplitude. The other propagates inward as a converging spherical compression wave. As long as reaching the center of spherical source, the converging spherical wave becomes a diverging spherical rarefaction wave, generating the negative pressure. This interpretation is easy to understand for the onedimensional situation due to the fact of undisturbed amplitudes in propagation. .10 shows that PA signal amplitude decays with the distance between center of sphere and detector increasing, even if the detected signal shapes have similar shape. The detected pressure amplitude is much weaker than its initial pressure amplitude. From Eq (1.12), it can be concluded that the amplitude of PA signal at the distance far from the spherical source is approximately proportional to the source size, but inversely proportional to the distance from the source object, R.

Fig. 1.11 shows that the signal shape varies with the size of spherical source. The temporal width of the PA signal from a spherical source is proportional to its diameter. The corresponding PA signal Fourier analysis shows a wide band spectrum. Moreover, the smaller the object is, the higher the frequency components in the generated PA signal spectrum are (see lower part of Fig. 1.11). Therefore, the ultrasound detector should be selected according to the size of imaged absorber.

Note that all these simulations start from a homogeneous distribution of initial pressure inside the source. In fact, it is directly related to the absorbed energy which determined by light source properties, medium properties and optical absorber properties. Therefore, the distribution of initial pressure is not uniform in practice and that could be taken into account in propagation models. This point will be discussed further in this thesis. K-wave toolbox is based on k-space solution of PA wave equation. This software is easy to use and provides rich examples to understand the principle of PA simulation and reconstruction. Here, the method to obtain k-space solution of PA wave equation is simply reviewed.

In a lossless medium, the equation of motion, equation of continuity, and equation of state can be expressed as first order form [START_REF] Treeby | k-Wave : MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[END_REF],

0 1 u p t       0 u t         2 pc  
where the initial pressure source is given by,

0 a p     0 0 p t   
Here, u is the acoustic particle velocity, 0  is the ambient density,  is the acoustic density, c is the thermodynamic sound speed, p is the acoustic pressure, 0 ( 0) p p t  is the initial photoacoustic pressure distribution,  is the Grüneisen parameter (the proportionality constant between the absorbed light and the initial pressure), a  is the optical absorption coefficient, and  is the optical fluence, all of which may be spatially varying. In k-wave toolbox, the pressure distribution is determined as the following process: the pressure distribution within the medium is used to calculate the spatial derivatives, which are used to update the corresponding velocity terms by use of a first-order FD. Then, the spatial derivatives of the velocity for each Cartesian direction are computed, which are used to update the values of the acoustic density within the domain by use of a first-order FD. Last, the pressure is updated by use of the appropriate equation of state. The process can be simply summarized as following figure, in which subscripts x and z represents the different Cartesian directions.

p p x p z     x z u u x z u x u z     x z  

Conclusion

In this chapter, the physical basis and related models of photoacoustic effect have been presented. The start point of this chapter is the general photoacoustic effect and then introduces two key points to completely understand PA effect in biological tissues: absorbers in biological tissues and the coupling between optical illumination and acoustic propagation. Related models governing the propagation of incident light and PA wave have been discussed. The advantages and disadvantages of related models have been stated. Up to now, two major forms of system have been developed [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF]: photoacoustic microscopy (PAM) and photoacoustic tomography (PAT). This chapter presents the principles of the two major forms of system from flowing three aspects: ultrasound detectors, laser source and image formation method, including conventional PA reconstruction algorithm and quantitative PA reconstruction algorithm.

PAI systems

PAM is based on a scanning-focused ultrasonic transducer to achieve a high-resolution volumetric image. Figure 2.1 shows PAM system developed by Lihong WANG's group [START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF]. Here the principle of this system is simply reviewed. A tunable dye laser produces a short-pulsed laser (the pulse width is ~10ns), which is delivered through an optical fiber to a conical lens, then forming a ring-shaped illumination pattern. It is then weakly focused into the tissue with the focal region coaxially overlapping the ultrasonic focus inside the tissue, see right part of Fig. 2.1. The sound waves acquired by ultrasonic transducer at each location can be converted into a 1-D depth-resolved image assuming constant sound velocity in soft tissue (~1540 m/s). Then, a raster scanning of the dual optical ultrasonic focuses in the horizontal (x-y) plane produces a 3-D image. [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF] In such this system, if the focal spot formed by optical system is smaller than the size of focal spot of ultrasound transducer, the spatial resolution is determined by the size of optical focal spot, so called OR-PAM (optical resolution photoacoustic microscopy) [START_REF] Beard | Biomedical photoacoustic imaging[END_REF][START_REF] Maslov | Optical-resolution photoacoustic microscopy for in vivo imaging of single capillaries[END_REF]. If the focal spot formed by optical system is larger than the size of focal spot of ultrasound transducer, lateral resolution of PAM is determined by the size of focal spot of ultrasound transducer, so called AR-PAM (acoustic resolution photoacoustic microscopy) [START_REF] Beard | Biomedical photoacoustic imaging[END_REF]. The axial resolution is inversely related to the bandwidth of the ultrasonic transducer (this point is to be explained below).

Shengfu LI Photoacoustic Imaging INSA Lyon 26

Photoacoustic tomography (PAT) is another major form of implementation of PAI system. PAT uses an unfocused ultrasonic transducer scanning around or an array of unfocused ultrasonic transducers to record photoacoustic waves. A representative of PAT system [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF][START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF] is shown Fig. 2.2. The principle of this system is simply reviewed hereafter. A laser provides a short-pulsed laser (the pulse width is ~10ns), which is expanded and homogenized to produce relatively uniform incident optical field. Once laser excitation has produced photoacoustic effect in the whole illuminated tissue, photoacoustic waves are generated. In this system, instead of a focused ultrasonic transducer, an unfocused ultrasonic transducer with a central frequency of ~10 MHz, scans around the objects in test with a radius of several centimeters and an appropriate step size to collect PA signals. After amplification and digitization, the signals are processed for image reconstruction. In Fig. 2.2., water is used as coupling medium to detect photoacoustic waves. [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF] From aforementioned descriptions, it can be seen that a general PAI system consists of three basic elements: optical excitation, acoustic detector and image formation method. This stands for both PAM and PAT. In the next sections, the three elements of PAI system will be detailed.

System axial resolution

Both in PAM and PAT system, the axial resolution is determined by acoustic transducer. In definition, axial resolution can be expressed as,  is the pulse width of time impulse response of ultrasound transducer. If the axial distance between two objects is less than the axial resolution, the two objects cannot be distinguished. Generally, the pulse width of time impulse response of ultrasound transducer can be expressed as a function of the bandwidth of transducer B, ~1 B  . So, axial resolution can be expressed also by,

1 2 s ax v R B 
System lateral resolution PAM system lateral resolution is determined by the focal spot formed by optical system, this point has been explained before. PAT system lateral resolution, R la , is determined by diffraction effect and is generally expressed by, , f L is the focus length and D is diameter of the transducer. Consider, for instance, a focused ultrasound transducer with a diameter of 2 cm and a focus length at 3cm. Suppose the center frequency at 6MHz with a 6MHz bandwidth, then then axial and lateral resolution of such a system are as follows. Lateral resolution could be improved with an image formation algorithm taking into account the scanning.

1 1500 125 2 2 6 10 s ax v R m B       2 62 1500 3 10 375 6 10 2 10 la focuslength R fnumber m activeaperture               .

Ultrasound detector

Ultrasound detectors of PAM system

Ultrasound detector parameters should be selected according to different objectives and applications. For microscopy, a focused transducer with high center frequency (several tens of MHz) is usually used to record PA signals to achieve high resolution. For example, a ultrasonic detector with 50 MHz center frequency and a 35 MHz nominal bandwidth was adopted in the system of Lihong WANG's group [START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF][START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF]). An axial resolution of ~15 μm, a lateral resolution of ~45 μm and an imaging depth of 3 mm were achieved.

Ultrasound detectors of PAT system

For PAT, there are two major forms of ultrasound detectors: (1) a focused or unfocused transducer, with a central frequency of several MHz, scanning at a closed shape (such as circle, spherical surface etc.). ( 2) an array of transducers with different shape (linear array, 2-D plan array (see Fig. 2.4), curved array (see Fig. 2.3)) [START_REF] Ermilov | Laser optoacoustic imaging system for detection of breast cancer[END_REF][START_REF] Vaithilingam | Three-dimensional photoacoustic imaging using a two-dimensional CMUT array[END_REF].

New sensors for PAI system

In published papers, new schemes detecting PA waves can be found. All these techniques try to resolve some typical problems of conventional acoustic detector, like limited bandwidth, physical dimension, laser incident position and axial resolution etc. Some groups detect PA waves by use of optical interferometry techniques, such as Beard et al (Zhang et al., 2009). A planar Fabry-Perot (FP) polymer film sensing interferometer is the center element of the system. The photoacoustic waves modulate the optical thickness of the FP interferometer (FPI), changing then locally its reflectivity, as shown in Fig. 2.5. By scanning a continuous wave focused laser beam (at a wavelength at which the FPI mirrors are highly reflective) across the surface of the FPI, the time-varying reflected optical power can be acquired as a 2D array. Other form of interferometer also can be used to record PA waves, such as Mach-Zehnder interferometer [START_REF] Paltauf | Photoacoustic tomography using a Mach -Zehnder interferometer as an acoustic line detector[END_REF]. Some other new sensors can also be found in published works, such as micro ring resonator [START_REF] Xie | Pure optical photoacoustic microscopy[END_REF], as shown in Fig. 2.6.

Laser source

Laser beam

Laser beam size should be selected according to different objectives and applications. Generally, a focused light beam is usually employed for PAM system to achieve high lateral resolution. On the opposite, a broad area illumination with a diameter of few centimeters is usually employed for PAT system to maximize the depth of exposure and acquire acoustic signals from deeper structures [START_REF] Razansky | Sensitivity of molecular target detection by multispectral optoacoustic tomography[END_REF].

Pulse width

Shengfu LI Photoacoustic Imaging INSA Lyon 30

To generate PA signals efficiently, laser pulse width  needs to satisfy the stress confinements and thermal confinement condition [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF][START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF]. The thermal confinement can be approximated by . Hence,  must be less than 100 ns to guarantee the stress and thermal confinement conditions.

Pulse energy

To guarantee safety of PA application in living humans, the pulse energy of the laser illumination must meet the safety criterion, given by American National Standards Institute (ANSI) (Laser Institute of America 2000) [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF][START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF]. Table 2.1 shows the maximum permissible exposure (MPE) values for the skin in the wavelength region of 400-1400 nm. According to the ANSI laser safety standards (Laser Institute of America 2000), MPE depends on the optical wavelength, pulse duration, exposure duration and exposure aperture.

) MPE a (mJ•cm -2 ) MPE b (mJ•cm -2 ) MPE b (mW•cm -2 )
400-700 700-1050 1050-1400 20 20×10 2(λ-700)/1000 100 1100•t 0.25 1100 •t 0.25 •10 2(λ-700)/1000 5500 •t 0.25 200 20×10 2(λ-700)/1000 1000 a: 10 -9 s≤t≤10 -7 s b: 10 -7 s≤t≤10s c: t>10s t: exposure duration.

Laser wavelength

In general, a wavelength region, from 650 nm to 900 nm (so called diagnostic window) (E. [START_REF] Zhang | Backward-mode multiwavelength photoacoustic scanner using a planar Fabry-Perot polymer film ultrasound sensor for high-resolution three-dimensional imaging of biological tissues[END_REF], is used for PAI because this optical range allows large imaging depths in biological tissue, due to the lack of strong absorption. It is noteworthy that Nd:YAG laser with a wavelength of 1064 nm can be often found in published papers due to the relative low cost and simplicity of use for biological tissues.

Wavelength-dependent behavior of scattering and absorption coefficient of biological tissues produce signal amplitude variation with adopted laser wavelength. In multiple wavelengths system, some key characters of absorption coefficient are essential. For example,

Shengfu LI Photoacoustic Imaging INSA Lyon 31
looking at the hemoglobin spectra (see in Fig. 1.2), the isosbestic point near the 800nm where oxy-and deoxy-hemoglobin show identical absorption (Jacques, 2013). Hence, the PA signal is proportional to the changes in hemoglobin concentration. The range including evident declining slope in the absorption spectra is often used to measure oxygen saturation (SO 2 ) changes of blood volume.

Reconstruction

The forward model in PAI refers to obtain PA signal from optical parameters of tested medium. Conversely, the inverse problem refers to obtain optical parameters from recorded PA signals by acoustic sensors. In published work, the reconstruction methods fall into two categories depending on the physical quantity to be reconstructed: (1) Initial pressure is reconstructed. Methods are so-called "conventional reconstruction"; (2) Optical parameters are reconstructed. Methods are so-called "quantitative reconstruction". DAS algorithm is to be explained with Fig. 2.7. Each element of a transducer array receives the signal from the same source at different time, ∆t i varies with i. In DAS algorithm, the In addition to the delay and sum process, the image quality can be further improved by utilizing the signal coherence information. The coherence of the delayed signals is estimated by the coherence factor (CF) [START_REF] Deng | Two-dimensional synthetic-aperture focusing technique in photoacoustic microscopy[END_REF], 2012;[START_REF] Liao | Optoacoustic imaging with synthetic aperture focusing and coherence weighting[END_REF][START_REF] Park | Adaptive beamforming for photoacoustic imaging using linear array transducer[END_REF][START_REF] Yoon | Enhancement of photoacoustic image quality by sound speed correction : ex vivo evaluation[END_REF][START_REF] Zemp | Photoacoustic imaging of the microvasculature with a high-frequency ultrasound array transducer[END_REF], which is defined as,

Conventional reconstruction algorithm

Delay and sum (DAS) algorithm

2 0 2 0 () () () e e N ii i N e i i i p t t CF t N p t t   D    D   (2.1)
The numerator in Eq. (2.1) represents the energy of the delayed signal, and the denominator denotes the total incoherent energy of the delayed signals. CF is maximal when the delayed signals are identical across the synthesized aperture. A high CF indicates that the image intensity should be maintained because the sources are in the direction of the synthesized beam. Otherwise, a low CF should be used to reduce the image intensity because of the presence of significant focusing errors. 

Fourier reconstruction



Then the wave equation can be transformed to an equation about P(k x ,z,ω), if combining with boundary value P(k x ,0,ω) ( Fourier transform of p(x,z,0)), one can obtain the solution of equation of P(k x ,z,ω). Once P(k x ,z,ω) has been determined, the last task is to resolve p(x,z,0) by use of inverse Fourier transform of P(k x ,z,ω).

The essence of Fourier reconstruction is based on the assumptions that the propagating PA wave transients p(x,y,t) can be described by a superposition of plane harmonic waves. As it propagates, each wave is an independent solution of the acoustic wave equation. The propagation corresponds to a phase modulation of the amplitude with the frequency  corresponding to the acoustic wavelength, while the modulus of the amplitude remains constant. These harmonic components superpose at time t=0 with complex pressure amplitudes 0,k P to form the original pressure distribution. And the detail about Fourier reconstruction can be found in references [START_REF] Burgholzer | Fourier Domain Reconstruction Methods in Laser Ultrasonics and Photoacoustic Imaging. 1st International Symposium on Laser Ultrasonics[END_REF][START_REF] Jaeger | Fourier reconstruction in optoacoustic imaging using truncated regularized inverse k -space interpolation[END_REF]. In time reversal reconstruction method, the recorded measurements are used as a time varying Dirichlet boundary condition at the position of the detectors, and then solving the PA wave equation [START_REF] Bal | Multi-source quantitative PAT in diffusive regime[END_REF]B T Cox, Kara, Arridge, & Beard, 2007),

Time reversal (TR) image reconstruction

2 2 2 2 2 2 2 1 ( , ) ( , ) ( , ) p r t p r t p r t x z v t        ( , ) ( , ) d p

r t p r t r S 

Therefore, in time reversal reconstruction method, the time evolution of the wave field in the medium is calculated from the measurements by using forward propagation model with zero initial conditions. The reconstruction result is then produced as the acoustic pressure distribution within the medium after time t. As the time reversal reconstruction is based on the forward propagation model, it is possible to include the medium parameters into the reconstruction simply by choosing the appropriate model parameters.

Filtered radial back-projection(FBP) algorithm

Filtered radial back-projection (FBP) algorithm can be used to reconstruct initial pressure. The algorithm can be expressed by (M. [START_REF] Xu | Time-domain reconstruction for thermoacoustic tomography in a spherical geometry[END_REF], ( , ) 1 ( , 0)

Discussion

The results of beamforming method is not initial pressure distribution, it only shows the position of PA wave source, and therefore cannot used to research quantitative initial pressures. The output of the other three algorithms is initial pressure (p 0 ) or absorbed energy (A(r)) (see in section 1.1.2).

Filtered radial back-projection (RBP) algorithm is an algorithm in time domain, the artifact introduced by RBP algorithm usually cannot be neglected [START_REF] Paltauf | Iterative reconstruction algorithm for optoacoustic imaging[END_REF]. Fourier reconstruction is often used for linear (planar) array and therefore only takes into account a part of waves travelling to transducer array. Some works (Cox & Beard, 2007) have been done to improve the quality of PA image recovered from linear array, still rather far from practical use in producing deposited energy distribution. ) and the geometry of transducer array (black points) that has a closed-shape geometry (evenly distributed across a circle). p(r,0) is taken as 1 and 0 when r locates inside and outside the absorber, respectively. The initial pressure distribution pattern is obtained from k-wave toolbox. In this geometry, the waves produced by initial pressure traveled outward and were recorded by transducer array. The recorded waves were used as input for reconstruction algorithm. d) shows that DAS introduces non-negligible artifacts around the absorbers and cannot reflect the true initial distribution. However, TR results can be used to recover the initial pressure distribution, and the artifacts can be neglected. In this thesis, a closed-shape array will be selected in our experiment (see in chapter 5).

Quantitative reconstruction

Several groups focus on the quantitative PAT, for example: (1) Department of Biomedical Engineering, University of Florida [START_REF] Yao | Finite-element-based photoacoustic tomography in time domain[END_REF], 2011;[START_REF] Yuan | Quantitative photoacoustic tomography: Recovery of optical absorption coefficient maps of heterogeneous media[END_REF][START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF][START_REF] Yuan | Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method[END_REF]; (2) Department of Medical Physics and Bioengineering, University College London [START_REF] Cox | Quantitative photoacoustic imaging : fitting a model of light transport to the initial pressure distribution[END_REF][START_REF] Cox | Quantitative Photoacoustic Image Reconstruction using Fluence Dependent Chromophores[END_REF][START_REF] Cox | Quantitative spectroscopic photoacoustic imaging: a review[END_REF]Cox, Arridge, & Beard, 2007;[START_REF] Laufer | Quantitative spatially resolved measurement of tissue chromophore concentrations using photoacoustic spectroscopy: application to the measurement of blood oxygenation and haemoglobin concentration[END_REF]; (3) Department of Mathematics, University of California (Gao, Zhao, & Osher, 2010).

Shengfu LI Photoacoustic Imaging INSA Lyon 36

Quantitative PA reconstruction consists in producing the spatial distribution of the optical properties of the tissue, including absorption coefficient µ a or scattering coefficient ' s  , from an absorbed energy distribution H m (r) that can be obtained by using conventional PAT. The process can be written as, Therefore, once H m (r) is measured, the task consists in minimizing the errors between the measured data and the modeled data, thus yielding the final absorption and scattering coefficients distribution.

        2 ' ' (

Conclusion

Two major forms of implementations of PAI system are presented in this chapter. The principal parameters have been detailed. Three characteristic elements of PAI system have been discussed: light source, ultrasound sensors and image reconstruction. The main existing conventional and quantitative reconstruction methods have been reviewed. This thesis will focus on quantitative reconstruction, thus this chapter can provide theory for following chapters and theory guide for experimental part. As mentioned in chapter 1, PA forward models involve multi-physics coupling between optics and acoustics. The generation of PA waves can be summarized as follows: when incident light propagates through a medium, part of energy is absorbed and produces an absorbed energy distribution, a fraction of absorbed energy converts into heat and produces a rise of temperature, which produces a pressure distribution, so called initial pressure, which produces acoustic waves travelling outwards.

This chapter is mainly divided into 2 parts. In the first part, two models for optical propagation in biological tissues are to be presented: (1) an analytical model of optical fluence for multiple cylinders embedded in an otherwise homogeneous turbid medium is to be presented in the section 3.1. This model has been published recently with the following reference: S. Li, B. Montcel, W. Y. Liu, and D. Vray, "Analytical model of optical fluence inside multiple cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium for quantitative photoacoustic imaging", Optics Express, vol. 22, no. 17, 2014. (2) a new efficient multigrid numerical method resolving the diffusion equation is developed and presented in section 3.2. These models will be used for quantitative reconstruction in the following chapters.

The second part of this chapter presents in section 3.3 two PA wave simulation models which are based on the optical models developed in section 3.1 and 3.2, respectively.

Analytical optical models in PAI

In the domain of diffusion optical tomography (DOT), analytical models of optical fluence for homogeneous and inhomogeneous media have been presented and validated with experiments and simulations. Patterson et al. developed analytical models for semi-infinite and finite homogeneous tissue slab [START_REF] Kienle | Improved solutions of the steady-state and the time-resolved diffusion equations for reflectance from a semi-infinite turbid medium[END_REF][START_REF] Patterson | Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties[END_REF]. Boas et al. developed analytical models for an inhomogeneity (sphere or cylinder) embedded in an otherwise homogeneous turbid medium. They introduced boundary conditions that take into account the index of refraction mismatch between the embedded inhomogeneity and the background medium [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF][START_REF] Walker | Photon density waves scattered from cylindrical inhomogeneities: theory and experiments[END_REF]. Their models have been validated experimentally [START_REF] Rocco | CW laser transilluminance in turbid media with cylindrical inclusions[END_REF]. Some models based on the Born approximation [START_REF] Boas | Three dimensional Monte Carlo code for photon migration through complex heterogeneous media including the adult human head[END_REF] were also developed, which are limited to small fluctuations of the optical properties. Ripoll et al. developed analytical models based on the Kirchhoff approximation [START_REF] Ripoll | Kirchhoff approximation for diffusive waves[END_REF][START_REF] Ripoll | Experimental demonstration of a fast analytical method for modeling photon propagation in diffusive media with arbitrary geometry[END_REF] by assuming that the total intensity at a certain point in the medium is equal to the sum of the incident field and the wave reflected from the plane tangent to the interface. Some works incorporate an extrapolated boundary into analytical models, which focus on the points located on the surface [START_REF] Arridge | The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis[END_REF]Zhang, Piao, & Bunting, 2012). These studies focused on DOT, and therefore the optical fluence is mainly investigated at the detector position. Relatively little works have been done on the analytical model of optical fluence for multiple inhomogeneities embedded in an otherwise homogeneous turbid medium, probably because this model can hardly be integrated in a DOT reconstruction process for biological tissues. Indeed, optical properties of biological tissues are very inhomogeneous. Furthermore, for DOT, the optical detectors are positioned onto the tissue surface, and therefore surface boundary conditions are mandatory. However, we believe that this model could be relevant for forward problems and inverse problems of PA imaging, since PA waves are only generated by absorbers inside tissues. The absorbers in tissues are mainly localized in anatomical structures. In particular, the main endogenous contrast for PA imaging is the blood vessels due to the strong absorption of hemoglobin. Hence, considering a model of multiple inhomogeneities of absorption could be more relevant in PA imaging than in DOT. Blood vessels are roughly cylindrical and have a much stronger absorption coefficient than the background tissues. Furthermore, hemoglobin concentration can be assumed to be uniform inside the vessel. Therefore, the blood vessels can be considered as "cylindrical inhomogeneities." In this section we present an analytical model for multiple parallel cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium to analyze the factors that influence optical fluence distribution.

Analytical models for such configuration are a well-known problem in acoustics, electromagnetics and the light scattering of particles [START_REF] Felbacq | Scattering by a random set of parallel cylinders[END_REF][START_REF] Lee | Scattering of evanescent wave by multiple parallel infinite cylinders near a surface[END_REF][START_REF] Lin | Sound scattering by a group of oscillatory cylinders[END_REF][START_REF] Linton | Multiple scattering by random configurations of circular cylinders: Secondorder corrections for the effective wavenumber[END_REF][START_REF] Mishchenko | Scattering of light by polydisperse, randomly oriented, finite circular cylinders[END_REF][START_REF] Wu | Multiple Scattering by Parallel Plasma Cylinders[END_REF][START_REF] Young | Multiple scattering by two cylinders[END_REF]. However, all these related works focused only on the analytical models corresponding to the physical quantity outside the "inhomogeneities". It was also assumed that a plane wave was incident on the "inhomogeneities" in acoustics, electromagnetics [START_REF] Felbacq | Scattering by a random set of parallel cylinders[END_REF][START_REF] Lee | Scattering of evanescent wave by multiple parallel infinite cylinders near a surface[END_REF][START_REF] Lin | Sound scattering by a group of oscillatory cylinders[END_REF][START_REF] Linton | Multiple scattering by random configurations of circular cylinders: Secondorder corrections for the effective wavenumber[END_REF][START_REF] Wu | Multiple Scattering by Parallel Plasma Cylinders[END_REF][START_REF] Young | Multiple scattering by two cylinders[END_REF]. In the research of the light scattering of cylindrical particles, only single scattering has been considered based on the assumption that the interaction of particles can be neglected [START_REF] Mishchenko | Scattering of light by polydisperse, randomly oriented, finite circular cylinders[END_REF]. However, in biological tissues, PA waves are generated by the local absorption contrast introduced by "inhomogeneities" in biological tissues, thus the optical fluence inside inhomogeneities is a key point to analyze the PA waves. The objective of this section is to present a model that focuses on the optical fluence both outside and inside the inhomogeneities embedded in turbid medium.

Our model arises out of the solution of the diffusion equation. The model represents the optical fluence distribution in tissues as a function of the reduced scattering and absorption coefficients. The migration of light through biological tissues in this model can be treated as a wave that is called a diffuse photon density wave (DPDW) in literature [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF]. In tissues with homogeneous optical properties, the incident DPDW of light source propagates unobstructed. However, the presence of "inhomogeneities" in the optical properties leads to a distortion, modeled as a scattered DPDW [START_REF] Walker | Photon density waves scattered from cylindrical inhomogeneities: theory and experiments[END_REF]. Light that diffuses through a medium including "inhomogeneities" can be considered as a superposition of the incident wave and scattered waves [START_REF] Cheng | Diffuse optical reflection tomography using continuous wave illumination[END_REF]. In view of high scattering properties of biological tissues, this incident wave cannot be modeled as a plane wave and the multiple scattering of "inhomogeneities" needs to be considered. The interplay between cylindrical inhomogeneities is incorporated by using new boundary conditions that take into account the multiple scattering of inhomogeneities. This model is validated by comparing with the numerical solution obtained by use of NIRFAST software (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF].

Shengfu LI Photoacoustic Imaging INSA Lyon 41

Geometry

The geometry of the problem is illustrated in Fig. 3.1. A point light source is incident on cylindrical inhomogeneities embedded in an otherwise homogeneous infinite turbid medium. The inhomogeneities have different optical parameter from background medium.

Incident model

In PA imaging the acoustic propagation occurs on a timescale several orders of magnitude longer than the heat deposition. Therefore the time-integrated absorbed power density (i.e., the absorbed energy density) is the quantity of interest [START_REF] Cox | Quantitative spectroscopic photoacoustic imaging: a review[END_REF]. In most biological tissues the optical fluence obeys the diffusion equation [START_REF] Walker | Photon density waves scattered from cylindrical inhomogeneities: theory and experiments[END_REF]. In an infinite homogeneous medium, the time-independent diffusion equation has the form [START_REF] Hielscher | Comparison of finite-difference transport and diffusion calculations for photon migration in homogeneous and heterogeneous tissues[END_REF], , in an infinite homogeneous medium for a point source at s r has been given by other studies [START_REF] Boas | Scattering of diffuse photon density waves by spherical inhomogeneities within turbid media: analytic solution and applications[END_REF][START_REF] Patterson | Time resolved reflectance and transmittance for the noninvasive measurement of tissue optical properties[END_REF]. In this study we use the incident optical fluence at an arbitrary position r inside the medium,

  ( , ) exp 4 inc s s s S r r jk r r D r r      (3.3)
In the presence of cylindrical inhomogeneities, it is natural to analyze the problem in cylindrical coordinates that are defined with respect to each center of the cylinders, as in Fig. 3.1. The z-axis of each cylindrical coordinate is the axis of each cylinder. The polar axis of all coordinates is parallel to the x-axis. Important notations are explained in Fig. 3.1, in particular i is index of i th cylinder and ' i represents an arbitrary remaining cylinders, with respect to the i th cylinder.

Based on the published works [START_REF] Arridge | The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis[END_REF]Jackson,1995), the cylindrical expansion of ( , ) inc s rr  with respect to the i th inhomogeneity center can be written as,

    2 0 exp ( ) cos( ) 2 i i i i i inc s n n n out S dp j n p z I r K r D                (3.4)
Where, out D is the diffusion coefficient of the background. out k is the wave number of the background. i  and i s  are the azimuth angles of a given position in the medium and of the light source, respectively, with respect to the i th inhomogeneity center.

Scattering model

If the individual cylindrical inhomogeneities are sufficiently far from each other, the interplay between cylindrical inhomogeneities can be neglected and the boundary conditions given by Li et al. [START_REF] Li | Fluorescent diffuse photon density waves in homogeneous and heterogeneous turbid media: analytic solutions and applications[END_REF] are valid.

If the cylindrical inhomogeneities are closer, the physical interpretation of their interplay is shown in Fig. 3.2. inc  , the optical fluence incident on the i th inhomogeneity, produces the first-order scattered wave,

  1 scat i 
(see left one in Fig. 3.2). In the same manner all the remaining cylinders produce first-order scattered waves,

  ' ' 1 scat i i  
, which are incident on the i th cylinder and produce the second-order scattered wave, (2) i scat  . Proceeding in this way, one can obtain the l-order scattered waves, () il scat  , with l from 1 to infinity. According to this physical interpretation, the second order of scattering,

(2) i scat  , can be interpreted as arising out of the consecutive scattering of the incident wave by two different inhomogeneities.

() il scat  arises out of the scattering of the incident wave by l consecutive scattering process. Therefore, () il scat  approaches 0 when l approaches infinity. This enables the truncation by considering only the first t orders of scattering to generate a solution, where t is a positive integer.

On the basis of the general solution given by Walker et al. [START_REF] Walker | Photon density waves scattered from cylindrical inhomogeneities: theory and experiments[END_REF] and of Eq. (3.4), the l-order scattered wave from the i th cylinder can be expressed as follows, (

i l i i l i scat n n n dp j n p z B p K x              (3.5) ) ( ) 0 exp( ) cos( ) ( ) ( )
and the sum of the first t orders of scattering waves from the i th cylinder i SC t , can be expressed as follows:

0 exp( ) cos( ) ( ) ( ) i t i i t i nn n SC dp j n p z B p K x             (3.6) With () 1 ( ) ( ) t i t i l nn l B p B p    .
Correspondingly, the fluence inside the i th cylinder has the form, 0 exp( ) cos( ) ( ) ( )

i t i i t i in n n n dp j n p z C p I y              (3.7) Here, 2 2 1 2 () ii out x p k     , 2 2 1 2 () i i i in y p k     , i in k
is the wave number inside the i th cylinder, ()

it n Bp , () () il n
Bp and ()

it n
Cp are unknown coefficients that are determined by use of boundary conditions.

Boundary conditions

Shengfu LI Photoacoustic Imaging INSA Lyon 44

As explained by Fig. 3.2, we can see that

i inc  produces   1 scat i  and that   ' ' 1 scat i i   produces   2 scat i  . Thus,   ' ' 1 i inc scat i i     produces   2 2 1 l i scat l i SC   
. Proceeding in this way, and considering the first t orders of scattered waves, we can see that ' '

1 i i t inc i SC    produces () 1 t i t i l scat l SC   
. Thus the fluence outside the i th cylinder can be written as:

' ' 1 i t i i t i t out inc i SC SC        (3.8)
It is noteworthy that it out  include three terms that are incident waves, the sum of the first t orders of scattering waves from the cylinder in question and the sum of the first t-1 orders of scattering waves from the remaining cylinders, respectively. Boundary conditions require that: (1) the flux normal to the boundary of the i th cylinder must be continuous; (2) the optical fluence must be continuous across the boundary of the i th cylinder. They can be written as:

i i i i i i i i i t i i t out out in in ii aa i t i t out in aa DD                                 
(3.9)

Solving the system of linear equations we have,

                      ' ' ' 2 ' ' 2 i out n n n b b in n n n b b it n i out out n b n b b in n b n b b D S a SC I y y D S a SC I y x S Bp D D K x I y x D K x I y y                   (3.10) and                       ' ' ' 2 ' ' 2 out n n n b b out n n n b b it n i out out n b n b b in n b n b b D S a SC K x x D S a SC K x x S Cp D D K x I y x D K x I y y                   (3.11) Where 22 i b out x a p k    , 22 i b in y a p k    and 22 i b s out z p k     .               '' ' 1 exp( ( ) ) 1 exp( ( ) ) ( ) i n n b s n m s m b m n ii i t n n b n m m i m S a I x j n m I z K z SC I x j n m K z B p                                            '' ' '' ' ' 1 exp( ( ) ) 1 exp( ( ) ) ( ) i n n b s n m s m b m n ii i t n n b n m m i m S a I x j n m I z K z SC I x j n m K z B p                              Where, 22 i s s out z p k     , ' ' 22 ii out i z r p k    , '
ii r and '

ii  are, respectively, the radial distance and the azimuth angle of the center of the i th inhomogeneity with respect to the center of the 'th i inhomogeneity.

()

it n
Bp and ()

it n
Cp are derived under the adoption of the addition theorems of Bessel functions (Jackson,1995). The details of the addition theorems are presented in Appendix A and Appendix B.

Therefore, the unknown coefficients (3.12) Considering Eq. (3.12), the stop condition of the recurrence process can be written as:

    1 1 1, 2,..., it n it n Bp iN Bp      (3.13)
 is an arbitrarily small quantity that depends on the required precision.

Validations

Our aim was to investigate the validity of this new model when used as a model to analyze the factors that influence optical fluence distribution in medium. The quantity of interest is the optical fluence inside the cylindrical absorption inhomogeneities, which represent the blood vessels. Indeed the PA signal is directly linked to the absorbed energy into the blood vessels, which is proportional to the product of the absorption coefficient and the optical fluence. We used NIRFAST software (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF] as the gold standard for the characterization of the optical fluence into cylindrical inhomogeneities.

Several numerical phantoms were defined so as to validate the model. They were all based on a cube ( 3 444 cm  ), and a point light source was positioned 0.1 cm ( ' 1 s  ) from the border. The distance between nodes was taken as 0.05 cm ( ' 0.5 s  ) for NIRFAST simulations. Optical properties of phantoms were close to the typical parameters of biological tissues (M. [START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF] in the near-infrared region. The reduced scattering coefficient was constant in the whole medium (background and inhomogeneities), '1 10

s cm    .
The first three sets of phantoms include two inhomogeneities to investigate the influence of different parameters of inhomogeneities (depth, size, and absorption) on the interaction between inhomogeneities. The last two phantoms include multiple inhomogeneities to validate our model for more complex cases.

For each phantom we simulated the optical fluence distribution with NIRFAST. This solution is denoted "NIRFAST" and is used as the gold standard. We then simulated the optical fluence distribution with the analytical model in each phantom, considered in this case as an infinite cube.  was set to 3 10  , which fixed a value of t max by use of Eq. (3.13). We evaluated the influence of the orders of scattering wave by simulating different analytical solutions with increasing values of t. The optical fluence distribution corresponding to, t=1, was denoted "1-order." We calculated the optical fluence distribution as the value of t increased from 2 to t max , these solutions were denoted "t-order". n in our analytical model arises from the separation of variables and should be infinity in theory. In the following simulation results, the value of n ranges from -20 to 20. The upper limit of n can be fixed by a method similar to Eq. (3.13). In this section, it was calculated by following equation, In order to evaluate quantitatively the error of the proposed model compared with the chosen standard, we calculated the mean relative error (MRE) with the following equation: 3.3(d). The corresponding depths of the c1 center are 2.1 cm and 2.3 cm, the depth of the c2 center in Fig. 3.3(a) and Fig. 3.3(e) is 1.5 cm. The absorption coefficient of all inhomogeneities is 0.8 cm -1 . The radius of all inhomogeneities is 0.3 cm. The corresponding fluence distribution along the line y=0 cm is shown in the second rows of Fig. 3.3(b) and Fig. 3.3(e). The curves of MRE versus the order of scattering are given in Fig. 3.3(c) and Fig. 3.3(f). The resolved t max values are 4 and 3, the difference of t max indicates that a higher order of scattering needs to be considered if the two inhomogeneities are closer. The t max values were determined as c1 was moved toward the -x direction along the dotted line in Fig. 3.3(a). Figure 3.4 shows the curve of t max versus the depth of the c1 center. This curve shows that a higher order of scattering needs to be considered if the depth of the c1 center decreases. This is due to the fact that the scattering wave is an outgoing wave, which is reduced as the distance of propagation increases.

Two inhomogeneities with different size

This set of simulation examples shows the interaction of two inhomogeneities with different radius. The position of the inhomogeneities (denoted by c1 and c2) and light source are shown in Fig. 3.5(a) and Fig. 3.5(e). The corresponding radius of c1 is 0.5 cm and 0.7 cm, the radius of c2 in Fig. 3.5(a) and Fig. 3.5(e) is 0.3 cm. The absorption coefficient of all inhomogeneities is 0.8 cm -1 . The corresponding fluence distribution along the line y=-0.5cm is shown in Fig. 3.5(b) and Fig. 3.5(f). The fluence distribution along the line y=0.5cm is shown in Fig. 3.5(c) and Fig. 3.5(g). The resolved t max values are 3 and 5. The curve of MRE versus the order of scattering is given in Fig. 3.5(d) and Fig. 3.5(h). The t max values were determined as the radius of c1 was changed from 0.1 cm to 0.7 cm. Figure 3.6 shows the curve of t max versus the radius of c1. This curve shows that the order of scattering to be considered increases with the increase in the radius of c1. This is due to the fact that the interaction of c1 and c2 is improved as the radius of c1 increases. In this set of simulations the center of two inhomogeneities is fixed because the source of scattering waves can be considered as the center of inhomogeneities according to Eq. (3.5). These curves of optical fluence indicate the higher the absorption coefficient of c1 is, the faster the optical fluence decays. This is due to the fact that the rise of absorption results in an increase of wavenumber. The resolved t max value is 4. The curve of MRE versus the order of scattering is given in Fig. 3.7(d) and Fig. 3.7(h). The t max values were determined as the absorption of c1 was changed from 0.1 cm -1 to 0.9 cm -1 . Figure 3.8 shows the curve of t max versus the absorption coefficient of c1. This curve shows that t max is slightly influenced by the absorption value. We did not note any significant influence of the absorption coefficient properties variation.

Multiple inhomogeneities

In this set of simulations, the first phantom had five identical cylindrical inhomogeneities (denoted c1 to c5 in Fig. 3.9(a)) embedded within the cube. They had a radius of 0.3 cm and an absorption coefficient µ a of 0.8 cm -1 . The second phantom had three different cylindrical inhomogeneities embedded within the cube. Their radius were set in the range 0.3 cm to 0.5 cm and their absorption coefficient in the range 0.4 cm -1 to 0.8 cm -1 (see Table 3.1). Figure 3.9(a) and Fig. 3.10(a) show, respectively, the horizontal cross sections of the first and second phantoms. Figure 3.9(a) shows the horizontal cross sections of the first phantom. The t max value was found to be 5. The configuration of the inhomogeneities was symmetrical about the line, y=0. Figure 3.9(b) and Fig. 3.9(c) show the optical fluence distribution through the horizontal cross sections corresponding to lines y=-0.7cm and y=0cm, respectively. We can clearly see that the analytical solutions obtained by only considering first order and first two orders of scattering waves had large errors. This means that the higher-order (larger than second order) scattering waves need to be considered in this case.

The MRE corresponding to the fluence distribution in the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ] is given in Fig. 3.9(d). The two remaining inhomogeneities had the same behavior because of the symmetry of the phantom. The curve of MRE versus the order of scattering used in the simulation confirms that a higher order of scattering has to be considered.

In the second phantom we considered a geometry including inhomogeneities with different parameters. The t max value was found to be 5. Figure 3.10(b), Fig. 3.10(c),and Fig. 3.10(d) show the optical fluence distribution through the lines y=0.7 cm, y=0 cm, and y=-0.7 cm, respectively. The analytical solutions obtained by only considering first order and first two orders of scattering waves still had larger errors. Simulation results showed that the analytical solution inside and outside the inhomogeneities was in good agreement with the NIRFAST solution when t=t max . The curves of the MRE corresponding to the area [-1.2cm, 1.2cm] × [-1.2cm, 1.2cm] versus the number of order, Fig. 3.10(e), confirm the observation. These results also demonstrate the relevance of the new analytical model for the investigation of real vasculature, which presents numerous blood vessels with different sizes, depths, and absorption coefficients. Indeed oxy-hemoglobin and deoxy-hemoglobin have very different absorption coefficients, and thus the oxygen saturation of blood vessels greatly influences its absorption coefficient. From the MRE curves we can see that the increasing number of the order of scattering considered, t, optimizes the analytical solution. MRE falls sharply when t goes from 1 to 2, and it falls more slightly when t goes from 2 to t max . This means that the scattered waves are weakened with increasing orders. As we mentioned previously, this property enables the convergence of the recurrence process.

To compare our analytical model in terms of computation speed, we implemented the algorithms with Matlab (R2011b) and they were run on an ordinary computer (Intel(R) core(TM) i7-2760QM CPU @ 2.40GHz). It took 642 s to run the first example by using NIRFAST software compared with 54 s with the analytical model. This is a 12-fold improvement in computation speed, which means that the analytical model is more likely to be used for real-time imaging when implemented in a faster computing environment. In all the results, the fluence was given within the interval [-1.2 cm, 1.2 cm] not within the interval [-2 cm, 2 cm], and the fluence distribution was normalized by its maximum value. In the area close to the border, there were larger errors due to different boundary conditions.

The index of refraction is considered as a constant in the whole media. If one wants to consider the variation of the index of refraction, the boundary conditions at the surface of a certain cylinder need to be modified to incorporate Fresnel reflection [START_REF] Walker | Photon density waves scattered from cylindrical inhomogeneities: theory and experiments[END_REF]. If the point source is replaced with a light beam, the source term can be considered as being composed of many point sources, which means that the analytical solution of a beam source can be obtained by an integral of the solution of the point source. It is noteworthy that our model can also be used to research the influence of scattering coefficient on optical fluence distribution, though the scattering coefficient in all numerical phantoms used in simulations has been assumed as a constant in whole medium.

It is noteworthy that the analytical model in this section is confined to the geometry of multiple infinite cylindrical inhomogeneities embedded into a homogeneous medium. However, the idea dealing with multiple inhomogeneities can be directly extended to other cases, for example multiple spherical inhomogeneities and multiple discal inhomogeneities.

Multigrid solver of the diffusion equation

As mentioned before, it is not possible to obtain an analytical solution for an arbitrary medium. Analytical models only exist for simple cases, such as cylindrical and spherical objects. Therefore, analytical model is not a flexible model. In most cases, researchers analyze the optical fluence distribution or absorbed energy distribution with numerical method.

Usually, finite element method (FEM) is used to resolve the forward solution required in inverse problem [START_REF] Cox | Estimating chromophore distributions from multiwavelength photoacoustic images[END_REF][START_REF] Cox | Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method[END_REF][START_REF] Razansky | Hybrid photoacoustic fluorescence molecular tomography using finite-element-based inversion[END_REF][START_REF] Yuan | Quantitative photoacoustic tomography[END_REF]. However, there are several practical challenges. Firstly, the optical fluence is described by the solution of diffusion equation (PDE) which is computationally demanding to solve. Second, the unknowns to be reconstructed are formed by the coefficients of the PDE, so the forward model is highly nonlinear, even when the PDE is itself linear. Finally, these problems typically are inherently three-dimensional (3-D) due to the 3-D propagation of energy in the scattering media being modeled. Meanwhile, the finite element method involves the matrix inversion while resolving the solution of diffusion equation. If the measurement region is large, this constitutes a large-scale inverse problem and poses some practical difficulties. Furthermore, all these works are implemented on a fixed finer grid. Although a finer fixed grid reduces error of numerical method and enhances the resolution, resolving problems at finer resolution also leads to slow convergence.

In this thesis, we will propose a multigrid-based inversion scheme for quantitative photoacoustic tomography reconstruction in which both the forward model and optimization method is to be expressed at multigrid (see in chapter 4). In this chapter, multigrid solver of diffusion equation is presented hereafter.

This section firstly presents basic idea of multigrid method and then explains the multigrid solver of diffusion equation.

Multigrid idea briefly sketched

Multigrid methods are usually used as solvers for linear equation, as Eq. (3.15), usually representing a discretization form of a differential equation.

Ax b  (3.15) The principle of multigrid can be found in associated references [START_REF] Garg | Multigrid Approach for Solving Elliptic Type Partial Differential Equations[END_REF][START_REF] Sharan | A Numerical Solution of Burgers' Equation Based on Multigrid Method Full Discretization of Equation[END_REF]Sundar & Stadler, 2014). In this scheme, the problem is resolved in multiple grids with different spacing. Let x (0) denote the solution of Eq. (3.15) at the finest grid and let x (q) be the solution at an arbitrary coarser grid, here q is a positive integer. Grid spacing corresponding to (q) th grid is 2 q times the spacing of the finest grid. The map of the related quantity from (q) th to (q+1) th (from (q+1) th to (q) th ) is based on the linear decimation matrix, noted by ( 1) () q q I  (linear interpolation matrix, noted by () ( 1)

q q I  ).
To have a better explanation of ( 1) () q q I  , let's take 1D case for an example. If x (q) consists of ( ) ( ) ( ) 1 2 7 , , , q q q x x x   , ( 1) () q q I  can be chosen as a form,

Shengfu LI Photoacoustic Imaging INSA Lyon 56     1 1 2 1 1 1 2 1 4 1 2 1 q q I       
Then the map from x (q) to x (q+1) has the form,

() 1 () 2 () ( 1) 3 1 ( 1) ( 1) ( ) 24 ( 1) () 3 5 () 6 () 7 1 2 1 1 1 2 1 4 1 2 1 q q q q q q q q q q q x x x x x x x x x x x                         
Likewise, in this case, one can understand the map from x (q+1) to x (q) ,

    1 1 2 11 1 2 2 11 2 1 q q I             

Two grids algorithm

Multigrid resolution of Eq. (3.15) can be understood as multiple recursive call of the two grid algorithm. Therefore, the basic principle of the two grid algorithm is the key point to understand multigrid method. The two grid algorithm is to be explained hereafter.

Suppose that an approximation x (q) , has been found for (q) th grid, that is

( ) ( ) ( ) q q q
A x b  (3.16) Then the residual corresponding to x (q) can be expressed as, ( ) ( ) ( ) ( ) q q q q r b A x  . x (q) is to be corrected by the defect correction, () q e , which can be obtained by resolving Eq. (3.17) at the (q+1) th grid and using linear interpolation operation,

    ( ) ( 1) 1 q qq q e I e    .   ( 1) ( 1) ( 1) ( ) () q q q q q A e I r     (3.17)
And thus x (q) is to be improved by Eq. (3.18),

( ) ( ) ( ) ( 1) ( 1) q q q q q x x I e    (3.18)
In this thesis, we use the recursion known as V-cycle multigrid [START_REF] Garg | Multigrid Approach for Solving Elliptic Type Partial Differential Equations[END_REF][START_REF] Sharan | A Numerical Solution of Burgers' Equation Based on Multigrid Method Full Discretization of Equation[END_REF]Sundar & Stadler, 2014), as shown in Fig. 3.11, in which a four-level grid is used as an example to illustrate the procedure of one update of x (0) . To implement one update of x (0) , V-cycle multigrid algorithm moves from fine to coarse grids (from (0) th to (3) th grid) by recursive call Eq. (3.16) and Eq. (3.17), and then backtracks to coarse grid (from (3) th to (0) th grid) by recursive call Eq. (3.18). Fig. 3.11 shows one update of x (0) . In practice, multiple updates which move forth and back fine to coarse grids are required to obtain a solution with certain precision, as shown in Fig. 3.12, in which the number of updates is noted by 1, 2, . The maximum number of updates is usually determined by required precision. Fig. 3.12 Multiple updates of x (0) . The grid levels are noted by (0) to (3). Red broken line represents multigrid algorithm recursively moves forth and back. The dotted lines represents grid with different spacing. The position of each update of x ( 0) is noted by black points along the green line with an arrow.

The number of updates is noted by 1,2,….

The root of multigrid methods

The very root of multigrid methods relies on the higher effectiveness of relaxation methods in the high frequency (oscillation components) dampening effect. If the error has high frequency components for certain grid spacing, then the relaxation attenuates it quickly. For the low frequency component (smooth components), the relaxation does not work as well. In multigrid scheme, the problem is to be re-defined on a coarser grid. The relaxation methods are to be more effective on the smooth components of error at the coarse grid spacing. Further explanation of the effectiveness of relaxation is consigned to Appendix C.

Multigrid solution of diffusion equation

According to the aforementioned statement, to use multigrid to resolve diffusion equation, two questions need to be answered: (1) How to prove the feasibility of multigrid in resolving the optical fluence? (2) How to choose a set of multigrid? These questions are to be answered in the following sections.

Feasibility

According to the last section, the migration of light through biological tissues can be treated as a wave that is called a diffuse photon density wave (DPDW) in literature. This point has been explained in last section. If the physical quantity can be considered as a wave in a finite medium, it also can be considered as consisting of a set of "modes" that includes different frequency, so-called smooth and oscillation components. Thus it provides the feasibility to resolve the physical quantity in multigrid which can be used as following process, oscillation components is to be resolved on fine grid, and smooth components is to be resolved on coarse grid. In mathematics, the diffusion equation has a form of an elliptic equation. Therefore, relying on the difference scheme (this point is to be explained below), the diffusion equation can be discretized in to a form as Eq. (3.15).

The choice of grids

In published papers, one has resolved the optical fluence in a fixed grid that meets the following condition, , to preserve the same oscillation components as in published papers. 2) According to theory of multigrid, to have a faster convergence speed, the grid spacing of other grids must be larger than the fixed grid adopted.

Multigrid solver

According to the multigrid scheme, the diffusion equation needs to be discretized as a form as Eq. (3.15). In this section, the finite difference method (FDM) is adopted to obtain the discretization form of the diffusion equation. FDM is one of the simplest and of the oldest methods to solve differential equations.

In the finite difference scheme for dealing with the time-independent case, the research medium is divided into a mesh and approximations of the solution are computed at the space points with an interval of h, as shown in Fig. 3 The approximations of the diffusion equations at the space points can be expressed as,

,, 2 , 1 , 1, 1, , 1 , 1 , 22 , 44 i j i j ij k k k k k i j i j i j i j i j i j a a i j D h S h D h D                    (3.19)
Here, k is iteration number. Assembling the approximations of the diffusion equations of all grid points in medium, one can obtain a linear equation as Eq. (3.15).

In numerical methods, one cannot consider infinite medium. It is necessary to adopt appropriate boundary conditions. Generally, there are three kinds of boundary condition: (1) Dirichlet boundary condition, its physical meaning is that the solution is known on the boundary, this means that detectors are placed at the boundary. ( 2 In our model, Dirichlet boundary condition is used, the fluence on an extrapolated boundary (at a distance of 2D from the medium boundary) is set as zeros. According to the works of Farrell (Farrell, 1992), this boundary condition is equivalent to adopting the third type of boundary condition on the true boundary.

Implementations of Multigrid solver

The solutions of diffusion equation on the q th level of grid are noted by Φ (q) (r). As shown in Fig. 3.11 and Fig. 3.12, multigrid solver of diffusion equation can be stated by two steps: (1) the diffusion equation is resolved in each grid; (2) the solutions obtained in first step is corrected with the order from the coarsest grid to the finest grid. To obtain precise solution, the grid spacing at the finest grid should meet the condition, It is worthy to be mentioned that the final objective of multigrid solver is to obtain the solution on the finest grid. Other levels of grid are adopted only to increase the convergence speed. The multigrid method can obtain a solution that has the same precision as the solution obtained by use of normal scheme with sufficient iteration.

By using finite difference approximation, diffusion equation can be expressed as MΦ=S, where M consists of medium properties and grid spacing, S is source term. Thus the first step can be expressed as follows,

                              1 11 n q q q n q q q q nq q q q q MS r S M M I r       
The relation between this set of equations can be explained as follows. One can obtain the solution     n q  on fine grid by resolving     () q qq MS 

. And then the residue on the fine grid   q r can be obtained by

          n q q q q r S M  , then     1 n q 
is resolved on coarse grid by using

            1 11 nq q q q q M I r    
, where     1  q q I refers to a linear decimation matrix.

As mentioned above, a high frequency error component, as compared with the grid spacing, will decay quickly. For the low frequency error component, the relaxation methods will not work as well as for high frequency. At the multigrid methods, the solution will be refined on a coarser grid. Once the coarse grid solution has been reasonably obtained, this solution is interpolated to a finer grid. Therefore, the residue   q r of the fine grid mainly includes the lower frequency parts relative to fine grid spacing, and the solution will be refined with coarse grid (q+1) in second step.

The second step is to correct the solution with coarse grid.

                1 11 n n q n q q q q I       Here,     1  q q I
denotes the corresponding linear interpolation matrix.

Validation

Multigrid solver of diffusion equation is validated by simulation example on a numerical phantom which consisted of three cylindrical inhomogeneities with parameters given in Table 3.2, embedded into a cube ( 3 444 cm  ) with absorption of 0.1cm -1 , the cross section of the numerical phantom is shown in Fig. 3.14(A). Scattering coefficient was constant in whole medium, 10cm -1 .

Figure 3.14(C) and Fig. 3.14(D) show profiles through the optical fluence distribution calculated by NIRFSAT and multigrid scheme. Close agreement of these curves indicates that the multigrid scheme can be used to solve the forward model of optical absorbed energy as finite element method. The calculation time of NIRFAST was 643 s. And multigrid scheme yielded a calculation time of 213 s, which is 33% of calculation time of NIRFAST. Furthermore, the agreement of these results indicates that the multigrid scheme can be used to obtain an accurate solution of the diffusion equation. The difference of optical fluence distribution near the border is caused by different boundary condition. In NIRFAST, third type boundary condition is adopted. In our multigrid scheme, we adopted the first type boundary condition, as mentioned before. For same geometry, the analytical model proposed in last section yields a calculation time of 57 s, which is 11-fold improvement in computation speed compared with NIRFAST.

Table 3.2 Parameter of the numerical phantom Inhomogeneity Position (x,y) Radius(cm) Absorption(cm -1 ) C1 (0.7, 0.7) 0.5 0.3 C2 (0.7, 0.7) 0.4 0.2 C3 (-0.7,0) 0.3 0.15 Fig. 3.14 Validation of home-developed multigrid solver of diffusion equation, (A) geometry of numerical phantom consisting of three absorber with absorption coefficient of 0.3cm -1 (larger absorber with radius of 0.5cm),0.2cm -1 (middle absorber with radius of 0.4cm) and 0.15cm -1 (small absorber with radius of 0.3cm) embedded in a homogeneous medium with absorption coefficient of 0.1cm -1 . µ s =10cm -1 . (B) Optical fluence distribution calculated with multigrid method. Profile through the optical fluence distribution calculated with NIRFAST(green), normal finite difference scheme (black) and multigrid (red) through the line y=0.7cm (C) and x=0.7cm (D).

Shengfu LI Photoacoustic Imaging INSA Lyon 62

PA signal model

According to the PA effect, a complete PA model needs to couple optical propagation and ultrasound propagation. This can be achieved based on the fact that the acoustic propagation occurs on a timescale several orders of magnitude longer than the heat deposition. Thus, one can consider the process as the acoustic signal is produced by a steady state optical absorbed energy distribution, as mentioned before. Numerical model of PA signals has been proposed widely. In this section, firstly, a pure analytical PA signal model for multiple cylinders is proposed based on our analytical model of optical fluence. Secondly, pure numerical simulation of PA wave based on the proposed multigrid solver of diffusion equation will be also presented.

Pure analytical model of PA waves for multiple cylinders

An analytical model of optical fluence has been given in section 3.1. The absorbed energy can be obtained by multiplying the optical fluence distribution by absorption coefficient map. Then one can obtain the pure analytical model of PA signal for multiple cylinders by incorporating absorbed energy model into following acoustic propagation equation, with the notations presented in chapter 1. shows PA signals corresponding to the three absorbers. The interactions between absorbers are not evident in this geometry. This is due to the fact that the size of absorbers is relatively small and thus the scattering of absorbers can be neglected compared with incident light fluence.

Figure 3.16(a) shows the three absorbers with different parameters that are given in Table 3 .3. Figure 3.16(b) shows signals corresponding to the three absorbers recorded by ultrasound transducers. The signal is deviated from "N-shape". This indicates that the interactions between absorbers are evident. This is due to the fact that the sizes of the absorbers are relatively large and thus the scattering of absorbers cannot be neglected. These examples state that it is necessary to consider optical part to simulate PA signals.

Furthermore, this simulation method can produce data which can be used to research quantitative photoacoustic reconstruction and furthermore extended into producing data which can be used to multispectral photoacoustic imaging by varying absorption coefficient in optical fluence model. This model has been published in the conference, Biomedical Optics, Miami, Florida United States, April 26-30, 2014. Once the optical distribution has been obtained, PA signal simulation can be implemented by incorporate () Ar into k-wave toolbox (see in chapter 1), which is a numerical simulation tool of PA signals.

Here, an example simulated with the method mentioned above is shown in Fig. 3.17. The signal is from a spherical source with a diameter of 2 mm, as shown in lower part of Fig. 3.17. This signal is produced by a geometry as shown in upper part of Fig. 3.17, in which each pixel has a volume of 0.2mm×0.2mm×0.2mm, the absorber is spherical volume in figure and the transducer is placed on the border of medium (linear shape in figure

Conclusion

In this chapter, an analytical model of optical fluence has been proposed for multiple cylinders embedded in an otherwise homogeneous turbid medium. The model is based on the diffusion equation and represents the optical fluence distribution inside and outside inhomogeneities as a series of modified Bessel functions. We take into account the interplay between cylindrical inhomogeneities by introducing new boundary conditions on the surface of inhomogeneities. The model is compared with the numerical solution obtained with NIRFAST software. The close agreement between the two methods used to simulate the optical fluence distribution into absorption inhomogeneities permits the use of our model as a forward model needed in quantitative PA imaging.

To develop efficient numerical solver for complex medium, multigrid method has been proposed in this chapter. We firstly proved the feasibility of multigrid solver of the diffusion equation from physics and mathematics consideration. Then a multigrid solver of the diffusion equation has been proposed. Simulation results show that the proposed multigrid solver is an efficient way to resolve the diffusion equation. As mentioned in chapter 2 of this document, the quantitative PA reconstruction can be expressed as,

        2 '' ( , ) ˆ, , , arg min as a s m a s uu H r H r      Here,   ' ,,
as Hr  is the forward model of absorbed energy, r refers to an arbitrary point in the medium, ˆa  and ' ˆs  are the recovered absorption and scattering coefficient map, respectively. Therefore, once H m (r) is measured, the task consists in minimizing the errors between the measured data and the modeled data, thus yielding the final absorption and scattering coefficients distribution.

In this chapter, H m (r) refers to the absorbed energy corresponding to single light source position. This leads to nonuniqueness problem for recovering µ a and ' s  simultaneously. To overcome the nonuniqueness problem, it is assumed that the scattering coefficient has been known as proposed in some other works [START_REF] Cox | Quantitative Photoacoustic Image Reconstruction using Fluence Dependent Chromophores[END_REF][START_REF] Cox | Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method[END_REF][START_REF] Yuan | Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton method[END_REF], thus only absorption coefficient needs to be reconstructed. Thus the reconstruction process can be written as,

      2 ˆ, arg min a a m a u H r H r  
Generally, there are two points to be discussed for such inverse problem: (1) Forward model; (2) Optimization method.

Two optimization methods are often used for the quantitative PA reconstruction: (1) Gradient-based methods are to optimize a cost function by a Newton-like method; (2) Fixedpoint iteration methods are to optimize a cost function by a simple iteration method. Such reconstruction schemes are based on a fine fixed grid. The finer the mesh is, the more accurate the solution is with however a larger number of variables to be reconstructed. Thus high computation complexity and slow convergence are inevitable problems for a fine fixed mesh. To overcome these problems, a multigrid-based inversion scheme will be proposed in this chapter.

This chapter is organized as follows. In section 4.1 and section 4.2, the forward models proposed in chapter 3 are used to combine with the two optimization methods often used for the quantitative PA reconstruction, gradient-based algorithm and the fixed-point iteration algorithm, respectively. In the gradient-based algorithm, to reduce the number of variables to be reconstructed, the medium is segmented into several regions, and then the optical parameters are assumed constant in each region. In section 4.3, a multigrid-based inversion scheme will be proposed.

Quantitative reconstruction: gradient-based

In this section, typical gradient-based algorithm is optimized by using the analytical model of previous chapter and by segmenting the medium into several regions. A simulation case including multiple cylindrical inhomogeneities is taken as an example to validate our approach.

Shengfu LI Photoacoustic Imaging INSA Lyon 68

Geometry of the numerical phantom

Figure 4.1 shows the geometry of the numerical phantom to be used in this section. A 2D rectangular background region (5×5 cm 2 ) includes three circular objects with different parameters (see in Table 4.1) representing a horizontal slice of a 3D phantom. A point light source is situated in the medium, noted by S in Fig. 4.1(a). This figure also shows the relative positions of the inhomogeneities, noted by c1, c2, c3. The absorption coefficient is 0.1 cm -1 in the background. The reduced scattering coefficient is a constant (10 cm -1 ) in whole medium. Figure 4.1(b) shows the absorbed energy distribution that was produced by NIRFAST software (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF], in which the distance between nodes was taken as 0.05 cm. To remove the effect of boundary conditions, the absorbed energy distribution corresponding to the area [-2 cm, 2 cm ] ×[-2 cm, 2 cm ] is used as testing data, as shown in Fig. 4.1.

Image segmentation

The positon and size of absorbers being required for analytic-based algorithm, an image segmentation algorithm is applied to the absorbed energy distribution map in order to clearly separate absorbers and background. The absorbed energy distribution is not uniform inside the absorbers neither in the background, because of the nonlinear optical propagation in turbid medium, as biological tissues. So a simple threshold is not able to segment the absorbers. Considering the absorption contrast between absorbers and surrounding medium, we have

Analytic-based algorithm

In this section, an analytical model is used to obtain the modeled data and a Newton-based method is used to solve the optimization problem for quantitative reconstruction of the optical parameters.

The start point of this algorithm is to guess an initial value of parameters to be reconstructed, here noted by X 0 . Actually, X 0 consists of variables to be reconstructed. The optimization can be expressed as,

  0 , ( ) m X OptimizationMethod H H X D (4.1) 00 X X X   D
The task of OptimizationMethod is to find the best ∆X around X 0 to minimize the difference between measurements m H and the modeled data 0 () HX . Then updating X 0 with ∆X and repeating the process to search new ∆X.

Here, a Newton-based method is selected to find ∆X, which can be expressed as, ( ) ( )

TT m J J I X J H H   D  
Here, with N the number of measurements. Thus N and NI satisfy the following condition, N≥NI. This ensures the existence and uniqueness of solution. Considering the large dynamic range of absorbed energy, here the logarithmic of absorbed energy distribution is used as follows,

                             

Simulation validation

The simulated absorbed energy is shown in Fig. 4.1(B). With the proposed algorithm above and the segmented results, the absorption coefficient is recovered from simulated absorbed energy. Table 4.2 shows the reconstruction results. The error is less than 7%. Even if the simulated measurement H m does not include noise, error is introduced by image segmentation process that divides the medium into several regions. This can be further validated by the recovered results shown in Table . 4.3. These results have been obtained by assuming that the position of absorbers is known exactly. If the position of the absorbers can be detected exactly, the recovered results have higher precision, as shown in Table . 4.3.

Although analytic-based algorithms are limited to some simple cases, there is an obvious advantage related to the fact that it is not necessary to divide the medium into mesh. It is also noteworthy that the typical gradient-based algorithm can also optimized by using multigrid method of previous chapter and by segmenting the medium into several regions. The manner of optimizing the typical gradient-based algorithm by using multigrid method and by segmenting is consigned to Appendix D.

In our simulation, the medium includes three regions, thus the computation of Jacobian matrix is relatively cheap. If the number of regions is large, the computation of Jacobian matrix is expensive. In this case, Quasi-Newton method can be used to further improve the efficiency of the proposed algorithm.

Quantitative reconstruction: Fixed point iteration

In this section, the multigrid solver proposed for the diffusion equation in chapter 3 is to be combined with a fixed point optimization to achieve quantitative PA reconstruction.

Algorithm

As mentioned before, absorbed energy can be written as, Here,  is an arbitrary small real value.

Simulation validation

Several simulation examples are provided to validate the proposed algorithm for different cases, recovering absorption coefficient from noisy and noise-free absorbed energy. The first simulation is to recover the absorption coefficient from noise-free absorbed energy shown in Fig. 4.1(B). With the proposed algorithm, the absorption distribution is reconstructed. Horizontal profiles of the absorption distribution through the lines y=-0.7cm, y=0cm and y=0.7cm show that the absorption has been recovered exactly, as shown in Fig. 4.3. Next example is to demonstrate the performance of our algorithm for noisy data. The noisy absorbed energy is produced with noise-free data added with normally distributed noise. Figure 4.4 shows reconstruction results from absorbed energy including noise. Calculated signal-to-noise ratio (SNR) is ~50 dB. The horizontal profiles across the reconstructed map show that all the absorption inside the absorbers has been well recovered. However, the recovered absorption near the border which locates far from the light source has larger errors, this is because that the optical fluence decreases rapidly with the distance from the light source. These results show that the region near light source is less severely affected by noise. This indicates that a wide beam may provide better results.

Wide beam

In the following examples, a wide beam is used as light source to produce the absorbed energy. The point light source in last simulation example is replaced with a wide light beam. Wide light beam is approximated by putting a point source at each node, represented as crosses (×) in Fig. 4.5(A). It is expected that a wide beam can provide higher noise tolerance than a point source. To validate this point, noisy absorbed energy with a same SNR value than with the point source situation is produced with noise free data in Fig. 4.5(B) added with normally distributed noise. Figure 4.6 shows the reconstruction results. Horizontal profiles of the absorption distribution through the lines y=-0.7cm, y=0cm and y=0.7cm show that the absorption has been well reconstructed.

To further confirm our expectation mentioned above, a simulation example with a lower SNR than last example is designed. Noisy absorbed energy, with SNR of 25 dB, is produced also with noise-free data in Fig. 4.5(B) added with simulated normally distributed noise. Figure 4.7 shows a reconstruction results from absorbed energy including noise. These recovered results show that the geometry with a wide light beam is less affected by noise compared the point source situation. This indicates that a wide beam can provide better results.

Quantitative reconstruction: multigrid-based

Multigrid methods have been explained in the previous chapter and applied to resolve the diffusion equation. The basic principles and advantages have been presented. Relatively little work has been done in multigrid-based inversion scheme. Thomas Dreyer et al have presented the application of multigrid methods to optimization [START_REF] Dreyer | Multigrid optimization in applications[END_REF]; they resolved linear-like problems with multigrid methods. However, the inverse problem in quantitative PA reconstruction is a nonlinear problem, thus their works cannot be extended into this area. Jong Chul Ye et al have proposed a multigrid-based algorithm for diffusion optical tomography (Ye, Bouman, Webb, Member, & Millane, 2001). Their algorithm is an effective method in term of convergence speed. However, they made a linear approximation of nonlinear cost function by use of Newton-like method and then used multigrid methods to resolve a linear system.

In our proposed pure multigrid-based inversion scheme, both the forward and inverse problems are resolved at different grid resolutions. The computation is substantially reduced by reducing the required number of iterations for evaluating the forward model and the optimization scheme at multiple resolutions.

Algorithm

For the sake of convenience, we will omit the variable r, which refers to an arbitrary point in the medium. The modeled data of absorbed energy is noted by () a H  , the measured absorbed energy is noted by m H , then the essence of reconstruction is to find the solution of the following equation, Let (0) a  denote the solution of inverse problem at the finest grid and let () q a  be the solution at an arbitrary coarser grid, here q is a positive integer. Grid spacing corresponding to (q) th grid is 2 q times the spacing of the finest grid.

Suppose that an approximation () q a  has been found for (q) th grid, that is

( ) ( ) ( ) () () q q q am HH   (4.5)
Then the residual corresponding to () q a  can be expressed as,

  ( ) ( ) ( ) ( ) q q q q ma r H H   .
According to multigrid idea, () q a  is improved or corrected by use the solution at the (q+1) th grid in order to achieve a faster convergence.

And let ( 1) ( 1) qq a e    be the exact solution at (q+1) th grid, then,

    ( 1) ( 1) ( 1) ( 1) ( 1) ( 1) q q q q q q aa H e H r           (4.6) Shengfu LI Photoacoustic Imaging INSA Lyon 77
If ( 1) q r  is chosen as the decimation of the residual of (q) th grid,

    1 ( 1) ( ) q qq q r I r    , then     ( 1) ( 1) ( ) ( ) ( ) () q q q q q q m a r I H H    , then one can obtain,         ( 1) ( 1) ( ) ( 1) ( 1) ( 1) ( ) ( 1) ( ) ( ) ( ) ( ) ( ) ( ) ( )
q q q q q q q q q q q q a q a q m a H

I e H I I H H              (4.7)
Then, letting ( 1) ( 1) ( ) ( 1) () q q q q a q a Ie    



, thus, ( 1) ( 1) ( 1) ( ) () q q q q a q a eI    



, and

  ( ) ( ) ( ) ( 1)
( 1) ( ) ( 1) ( )

q q q q q q a a q a q a II          (4.8)
Multigrid algorithm improves () q a  with ( 1) q a   by use of Eq. (4.8). Multigrid optimization is implemented by recursively applying the equations from Eq. ( 4.6) to Eq. ( 4.8), the recursive process has been explained in chapter 3.

In following simulations and experiments, Equation (4.9) is used to measure the residual of recovered solution.

()

ma m HH E H    (4.9)

Grid spacing

To have a better understanding for the grid spacing, an example of 2 -level multigrid is given hereafter, noted by (q) and (q+1). It is clear that the number of nodes of the finer one is 11×11=121, the number of nodes of the coarser one is 6×6=36. Thus the node spacing h has the following relation, h (q+1) =2×h (q) . And the advantages of multigrid are summarized in the following figure. More details about multigrid principles have been given in chapter 3. According to simulations in last sections, the geometry with a large beam can produce better results than with the point source. Therefore, in this section, our simulations will be based on a wide beam. The proposed algorithm described above is validated with simulated data. The geometry is shown in Fig. 4.5(A) where three circular inhomogeneities were embedded in an otherwise homogeneous medium (5×5 cm 2 ), and a wide light beam was adopted.

Simulation validation

Figure 4.5(A) shows the relative positions of inhomogeneities that are noted by c1, c2 and c3. The inhomogeneities parameters are given in Table 4.1. Figure 4.5(B) shows the absorbed energy distribution that was produced by FEM (Dehghani et al., 2010;[START_REF] Jermyn | Fast segmentation and high-quality three-dimensional volume mesh creation from medical images for diffuse optical tomography[END_REF].

To explain the advantages of multigrid-based scheme in term of convergence speed, the number of levels of grid was taken as 4 (the number of grid levels is discussed below). The absorption coefficient is recovered with multigrid-based scheme and fixed grid method. The recovered results are compared, as shown in Fig. 4.8 where (A) shows the recovered reconstruction map. Close agreement between true and recovered absorption coefficient in Fig. 4.8(B) -(D) indicates that multigrid-based algorithm can be used to reconstruct the absorption coefficient as normal scheme (based on the fixed grid). Figure 4.8(E) and Figure 4.8(F) show that multigrid-based algorithm converges faster than fixed grid algorithm. In this simulation, Equation (4.10) was used as stop condition of iteration process in multigrid and fixed grid algorithm. Here, E n is the error, defined by Eq. ( 4.9), at n th iteration. The per-iteration times of multigrid-based method is longer than fixed grid algorithm, the faster convergence of the multigrid algorithm is due to the substantially fewer required iterations. Theoretically, the larger the number of level of grid is, the faster the convergence is, as shown in Fig. 4.9. These curves are obtained by use of the absorbed energy distribution shown in Fig. 4.5.

Next simulation is to test the proposed algorithm with noisy data that is produced by the noise-free absorbed energy added with normally distributed noise. In the following examples, a 4 levels grid is used to reconstruct absorption coefficient. Figure 4.10 shows reconstruction results from absorbed energy including noise, SNR is ~25 dB in whole medium.

Conclusion

In this chapter, different QPAT algorithms have been presented. The proposed algorithms aim at optimizing the reconstruction from the following three aspects: (1) forward model; (2) number of variables to be reconstructed; (3) optimization scheme. Following several trials, an efficient multigrid-based inversion scheme has been presented, in which both the forward model and inverse problems are resolved by a multigrid method. This scheme improves the efficiency of algorithm due to the fact that substantially fewer iterations are required compared with the algorithm based on a single fixed grid. Simulation and experimental results show that multigrid-based inversion scheme has a better efficiency than compared with normal scheme (based on a fixed grid) and can obtain a reliable result with the same precision as fixed grid algorithm. . In this chapter, analytical-based and multigrid-based algorithms proposed in chapter 4 of this document are validated with experimental data. The experiment has been realized in collaboration with Dr. Zhen Yuan, University of Macau. In section 5.1, the experimental system will be presented; Section 5.2 will present the experimental phantom and the acquired PA signals. Section 5.3 will present the conventional PA reconstruction result which represents the initial pressure or absorbed energy. Section 5.4 will validate the quantitative reconstruction algorithms based on the absorbed energy which conduct to the absorption map of the medium.

Experiment system

Beamforming result

Firstly, the acquired PA signals are processed to build a PA image. This is completed by using delay and sum algorithm for circle shape array. Result is shown Fig. 5.4.

Absorbed energy distribution

To calculate the absorbed energy distribution within the phantom, we used the time reversal algorithm presented in chapter 2. The results recovered with the time reversal algorithm represents the initial pressure (absorbed energy) distribution. The acoustic speed is taken as 1495 m•s -1 (average sound speed measured for the phantom). The grid spacing in x and y direction is ~0.03 cm. After time reversal reconstruction, one can obtain initial pressure distribution (absorbed energy), as shown in Fig. 5.5. The next step is to determine a calibration factor taking into account the sensibility of US transducers.

Quantitative PA reconstruction

Quantitative PA reconstruction is based on the absorbed energy distribution obtained in last section. A calibration step is necessary before quantitative PA reconstruction.

Calibration

It is assumed that the reconstructed results by time reversal algorithm in last section is proportional to true absorbed energy, thus the quantitative PA reconstruction process can be rewritten as,

      2 ˆ, arg min a a m a u H r C H r    
with C a calibration factor that is assumed as a constant to calibrate the modeled data. This is related to the sensibility of the acoustic transducer and the gain of the acquisition system. Then the calibration factor can be defined as,

m H C H  .
Here, H m represents the recovered absorbed energy shown in Fig. 5.5. H is modeled data which is obtained by calculating the absorbed energy with the forward models for experimental phantom shown in Fig. 5.2. Its parameters have been given in calibration factor is to be determined by dividing the mean value of H m with the mean value of H.

Analytic-based reconstruction

Forward model

To adapt to experimental geometry, a wide beam is used here in the analytical forward model. In chapter 3, a point light incident source has been chosen. The differences between the two geometries are simply reviewed hereafter.

The incident beam is greater than the size of embedded objects. Therefore, the incident beam is assumed as a planar wave, this point is further explained in Appendix E. In our experimental phantom, the light beam is incident on the phantom from the top. The distortion of incident planar wave is determined by the part of inhomogeneities located on the wave front, which have discal shape. Thus the problem can be dealt with 2D models. As mentioned previously, this incident light is scattered by the embedded objects, then the general solution of scattered waves, expressed in Eq. (3.5) for a point light source, need to be modified as, ( )

( ) exp( ) ( ) i l i l i i scat n n n B j n H x          
(5.1) Correspondingly, the general solution of optical fluence inside the object, expressed in Eq. (3.7) for a point light source, needs to be written as, (

i l i i l i in n n n j n C j y           ) ( ) exp( ) ( )

Results

As mentioned before, image segmentation is necessary to divide the medium into several regions due to the constraints of the analytical approach. The results of segmentation of the absorbed energy distribution are shown in Fig. 5.6. Pure scattering target (c1) is not detected by the segmentation algorithm. It is also assumed that optical parameters in each closed region are constant. The reduced scattering coefficient is assumed as a constant, , s µ = 10 cm -1 . The calibration factor is 30.9. After applying the proposed analytical-based algorithm in section 4.1 of this document, the absorption of each region can be recovered. The reconstructed results are listed in Table 5.2. The error of recovered absorption for c1 and c2 is less than 15%, which means that analyticalbased algorithm proposed in section 4.1 in this document could be used to quantitative PA reconstruction. These results also indirectly provide experimental validations of the analytical-based forward models proposed in chapter 3.

Multigrid-based algorithm

In previous chapter, we have shown that multigrid can be used as a solution of forward model and also as an efficient optimization method for quantitative PA reconstruction.

The proposed algorithm is now applied to experimental data. The calibration factor is 31. Close agreement of the results between multigrid-based algorithm and fixed grid based algorithm in Fig. 5.7(B) -(D) indicates that multigrid-based algorithm can achieve same precision with the fixed grid based algorithm. Both of them achieve a precision with an average error of 0.034; which is less than 10%. Figure 5.7(E) and Figure 5.7(F) show that multigrid-based algorithm converges faster than fixed grid algorithm. The multigrid algorithm reduces the number of iterations required substantially, the number of iterations of multigrid algorithm is about ~2 times fewer than the fixed grid algorithm under adopting 2 grid levels. The per-iteration times of multigrid method is longer than fixed grid algorithm. Faster convergence of the multigrid algorithm is due to substantially less iteration required. Finally, the multigrid algorithm takes ~1.6 times less than the fixed algorithm to meet the condition expressed in Eq. (4.11). These results also indirectly provide an experimental validation of the multigridbased forward models proposed in chapter 3.

Conclusion

In this chapter, different QPAT algorithms have been validated with experimental data. Experimental results show that multigrid-based inversion scheme has a faster convergence speed compared with fixed grid based algorithm. Moreover, multigrid algorithm can provide a reliable result that has the same precision as the fixed grid algorithm. Experimental results also shows that the analytic-based algorithm proposed in chapter 4 can be used in simple case as shown experimental phantom shown in this chapter.

The difference between DAS (delay and sum) and TR (time reversal) has been explained in chapter 2 of this document. In this chapter, their difference has been further validated with experimental data. Meanwhile, the experimental results in this chapter also indirectly provide experimental validations of the forward models proposed in chapter 3. This thesis focuses on quantitative PA reconstruction. To achieve quantitative information about medium, a forward model of optical fluence and an optimization method are required. The general objectives of this thesis are to present efficient forward models and develop new optimization schemes.

In the first chapter, existing models of optical fluence have been simply reviewed and compared. Published works have proven that Monte Carlo-based model of photon migration in biological tissues is an efficient model. Due to its generality and few of hypothetical conditions, this method has been used as a gold-standard in this area. However, this kind of model is time-consuming due to its low computational efficiency in comparison to a few seconds for solving the diffusion equation or minutes for solving RTE. RTE is usually solved by numerical method. However, it needs both spatial and angular discretization to solve the RTE by numerical methods. Considering these factors, a time-independent diffusion equation modeling is selected in this thesis.

In the second chapter, two forms of PAI system, PAM and PAT, are presented. Both PAM and PAT system consist of three elements: light source, acoustic sensor and image formation. Existing image formation methods of PAT include conventional PA reconstruction and quantitative PA reconstruction. Usually, the finite element method (FEM) is used to obtain the modeled data required in inverse problem of quantitative PA reconstruction. However, if the measurement region is large, the implementation of FEM causes some practical difficulties. In the second part of the second chapter, existing models of PA signals have been simply reviewed. Most of these models start from a homogeneous distribution of initial pressure. In fact, initial pressure is directly related to absorbed energy which determined by the light source and medium properties. So, the distribution of initial pressure inside absorbers and in the surrounding medium is usually not homogeneous in practice. To overcome the difficulties mentioned above, this thesis proceeds in following ways.

We developed an analytical model for multiple parallel cylindrical inhomogeneities embedded in an otherwise homogeneous turbid medium to analyze the factors that influence optical fluence distribution. We consider that this model should be useful for biological applications when blood vessels are modeled as piece-wise cylinders. Our model arises out of the solution of the diffusion equation. The model represents the optical fluence distribution in tissues as a function of the reduced scattering and absorption coefficients. The migration of light through biological tissues in this model can be treated as a wave that is called a diffuse photon density wave (DPDW) in literature. In tissues with homogeneous optical properties, the incident DPDW of light source propagates unobstructed. However, the presence of "inhomogeneities" in the optical properties leads to a distortion, modeled as a scattered DPDW. Light that diffuses through a medium including "inhomogeneities" can be considered as a superposition of the incident wave and scattered waves. In view of high scattering properties of biological tissues, multiple scattering of "inhomogeneities" needs to be considered. The interplay between cylindrical inhomogeneities is incorporated by using new boundary conditions that take into account the multiple scattering of inhomogeneities. This model has been validated by comparing with the numerical solution obtained by use of NIRFAST software. Once analytical model of optical flounce has been given, the absorbed energy can be obtained by multiplying the optical fluence distribution by absorption coefficient map. Then one can obtain the pure analytical model of PA signal for multiple cylinders by incorporating the analytical optical energy model into PA wave equation. Simulation results of pure analytical model of PA signal show that the obtained signal is deviated from "N-shape". This is due to the fact that the scattering of absorbers has evident influence on PA signals. Therefore it is important to consider optical propagation to simulate PA signals. Furthermore, these simulation methods can produce data which can be used to research quantitative photoacoustic reconstruction with possible introduction of multispectral illumination by varying absorption coefficient in optical fluence model.

Our analytical model has been applied to quantitative PA reconstruction and satisfactory results have been presented after segmentation of the acquired image in several regions where optical parameters are supposed constant. However, complex media cannot be segmented into several regions. In such cases, optical parameters for each point need to be recovered, and only numeric-based algorithms can be used to recover optical parameters. According the numerical theory, high computation complexity and slow convergence are inevitable problems. To overcome these problems, this thesis presents a multigrid-based inversion scheme.

Multigrid methods have been applied successfully to resolve partial differential equations (PDE). Then, a pure multigrid -based inversion scheme has been proposed in this thesis. In our proposed pure multigrid-based inversion scheme, both the forward and inverse problems are resolved at different grid resolutions. The computation is substantially reduced by reducing the required number of iterations. Simulation and experimental results show that multigrid-based inversion scheme has a better efficiency compared with fixed grid algorithm and can obtain a reliable results that have the same precision than with fixed grid algorithm but with significant reduced computation time.

Perspectives of this work could concern more general analytical model for absorbers like generalized cylinder model to fit with real 3D vascular tree. This will conduct to find more general analytical solution of the diffusion equation adapted to biomedical PA imaging. The authors will also continue to research new reliable segmentation algorithms of PA imaging which take into account a priori information concerning absorber shapes, geometry and topology of the medium. Multigrid QPAT scheme will be also validated with various experimental situations. In this context, to find the best adapted light illumination of the medium will be a challenge to be able to image and reconstruct biological tissues in 3D. Another point is to improve our multigrid-based inversion scheme to adapt to the data corresponding to multiple optical wavelengths.

B. Addition theorems of Bessel function

This annex simply explains the addition theorems of Bessel functions. It can be found in many books. Here, the main reference is WATSON's book named by "A treatise on the theory of Bessel functions". The generating function of the Bessel coefficients is

 

exp ( 1 ) 2 z t t  It can be expanded into a series of powers of t. The coefficient of n t in the expansion is called the Bessel coefficient of argument z and order n, denoted by the symbol () n Jz .

 

exp ( 1 ) / 2 ( )

n n n z t t t J z       (1)
Multiplying the expansion Eq. ( 1) by

 

(2)

The integrals on the right all vanish except the one for which m = n; and so we obtain the formula.

C. Effectiveness of relaxation schemes for resolving differential equations

Many relaxation schemes for resolving differential equations have the smoothing property, where oscillatory modes of the error are eliminated effectively (the convergence factor λ approaches 0), but smooth modes are damped very slowly (the convergence factor λ approaches 1), as shown in Fig. 1. Thus the residual norm is damped slowly after certain number of iterations, as shown in Fig. 2.

D. Gradient-based reconstruction optimized with multigrid

In this algorithm, the numerical multigrid solver proposed for solving the diffusion equation in chapter 3 is used to obtain the modeled data and a Newton-based method is used to solve the optimization problem for quantitative reconstruction of the optical parameters. It can be written as,

E. Broaden beam in experimental setup

If a broaden beam is incident on a homogeneous medium, the optical fluence distribution can be expressed by Lambert law like equation, as following figure shows (LI and [START_REF] Wang | Multiscale photoacoustic microscopy and computed tomography[END_REF]. From following equations, broaden beam produces a plane wave which propagates in z-direction.

As mentioned before, the incident wave adopted in experimental setup described in chapter 5 can be considered as plane wave. In our experimental phantom, the light beam is incident on the phantom from the top. Considering the cylindrical targets in phantom are parallel to each other, the distortion of plane wave is determined by the discal inhomogeneities located on the wave front. Therefore, the problem can be dealt with 2D model.

  Fig. 0.1 PAI of whole body of small animal (Brecht et al., 2010) ......................................... 2 Fig. 0.2 PAI of small animal brain showing the changes of optical absorption (Yang & Wang, 2009) .............................................................................................................................. 2 Fig. 0.3 An example of PAI application to breast cancer (Ermilov et al., 2010) ................... 3 Fig. 0.4 An example of PAI of the skin (E. Z. Zhang et al., 2009) ........................................ 3 Fig. 1.1 Illustration of Bell's experiment. A represents the transmitter diaphragm. B represents the lampblack receiver with hearing tube. See in the link: http://digital.library.cornell.edu/ ................................................................................................ 8 Fig. 1.2 Absorption spectra of different components of biological tissues. From: omlc.ogi.edu/spectra .................................................................................................................. 9 Fig. 1.3 The photon travels in a random trajectory. From: http://oilab.seas.wustl.edu/ ..... 11 Fig. 1.4 A general flow chart of a Monte Carlo program .................................................... 12 Fig. 1.5 Optical fluence distribution produced by a MC-based package. Photons are incident from the top-center of the medium. The medium is a homogeneous turbid slab, 1 ' 1 0.2 , 2 as cm cm   

Fig.

  Fig. 1.6 Optical fluence at point P is the sum of incident and scattered fluence by embedded inhomogeneity C ..................................................................................................................... 17Fig. 1.7 Optical fluence distribution obtained by NIRFAST for a homogeneous turbid medium with 1 ' 1 0.1 , 10 as cm cm   

  Fig. 2.8 A linear array is positioned at z=0. x-axis represents the linear array. z-axis represents depth. The recorded signals can be expressed as p(x, 0, t). .................................... 32 Fig. 2.9 Illustration of notations to be used in TR algorithm. Transducer array has a closed geometry (evenly distributed across a circle). d i represents the position of i th element of transducer array. A represents PA wave source. Purple circle represents surface of medium, noted by S. ............................................................................................................................... 33 Fig. 2.10 Illustration the difference between DAS and TR method. (a) Initial pressure distribution and the geometry of transducer array (black points). the recovered results of TR (b) and DAS (c). (d) Horizontal cross slice through the line, x=0 mm, of initial pressure, recovered results of TR and DAS. ........................................................................................... 35 Fig. 3.1 Horizontal cross section of the geometry. Cylindrical inhomogeneities of absorption coefficient embedded in an otherwise homogeneous infinite turbid medium. i is index of i th cylinder and i' represents an arbitrary remaining cylinders. Cylindrical coordinate systems have their origins at the center of each inhomogeneity. The source is noted by S. P indicates an arbitrary point within the medium, and its coordinates values with respect to each cylinder coordinate system are indicated. The z-axis is along the axis of each inhomogeneity and comes out of the cross section. ......................................................................................... 41 Fig. 3.2 Physical interpretation of the different orders of scattering waves by the i th inhomogeneity. Left: the incident wave produces the first-order scattered wave by the i th inhomogeneity; Right: the sum of first-order scattered waves from the remaining cylinders produces the second-order scattered wave by the i th inhomogeneity ...................................... 43 Fig. 3.3 . Comparison between the analytical model and NIRFAST for two inhomogeneities with different depth. Top row ((a) and (d)): Horizontal cross section of the first numerical phantom of two cylindrical inhomogeneities with different depth (depths of the c1 center in (a) and (e) are 2.1 cm and 2. 3cm, depths of the c2 center in (a) and (e) are 1.5cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (e)): Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1 and c2. Third row ((c) and (g)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. .......................................................................... 47 Fig. 3.4 Curve of t max versus the depth of the c1 center. ...................................................... 48 Fig. 3.5 Comparison between the analytical model and NIRFAST for two inhomogeneities with different radius. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different radius (radius of c1 in (a) and (e) is 0.5cm and 0.7cm, respectively, radius of c2 in (a) and (e) is 0.3cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) and (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row ((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. ................................................................................................................... 48

  Fig. 2.8 A linear array is positioned at z=0. x-axis represents the linear array. z-axis represents depth. The recorded signals can be expressed as p(x, 0, t). .................................... 32 Fig. 2.9 Illustration of notations to be used in TR algorithm. Transducer array has a closed geometry (evenly distributed across a circle). d i represents the position of i th element of transducer array. A represents PA wave source. Purple circle represents surface of medium, noted by S. ............................................................................................................................... 33 Fig. 2.10 Illustration the difference between DAS and TR method. (a) Initial pressure distribution and the geometry of transducer array (black points). the recovered results of TR (b) and DAS (c). (d) Horizontal cross slice through the line, x=0 mm, of initial pressure, recovered results of TR and DAS. ........................................................................................... 35 Fig. 3.1 Horizontal cross section of the geometry. Cylindrical inhomogeneities of absorption coefficient embedded in an otherwise homogeneous infinite turbid medium. i is index of i th cylinder and i' represents an arbitrary remaining cylinders. Cylindrical coordinate systems have their origins at the center of each inhomogeneity. The source is noted by S. P indicates an arbitrary point within the medium, and its coordinates values with respect to each cylinder coordinate system are indicated. The z-axis is along the axis of each inhomogeneity and comes out of the cross section. ......................................................................................... 41 Fig. 3.2 Physical interpretation of the different orders of scattering waves by the i th inhomogeneity. Left: the incident wave produces the first-order scattered wave by the i th inhomogeneity; Right: the sum of first-order scattered waves from the remaining cylinders produces the second-order scattered wave by the i th inhomogeneity ...................................... 43 Fig. 3.3 . Comparison between the analytical model and NIRFAST for two inhomogeneities with different depth. Top row ((a) and (d)): Horizontal cross section of the first numerical phantom of two cylindrical inhomogeneities with different depth (depths of the c1 center in (a) and (e) are 2.1 cm and 2. 3cm, depths of the c2 center in (a) and (e) are 1.5cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (e)): Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1 and c2. Third row ((c) and (g)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. .......................................................................... 47 Fig. 3.4 Curve of t max versus the depth of the c1 center. ...................................................... 48 Fig. 3.5 Comparison between the analytical model and NIRFAST for two inhomogeneities with different radius. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different radius (radius of c1 in (a) and (e) is 0.5cm and 0.7cm, respectively, radius of c2 in (a) and (e) is 0.3cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) and (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row ((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. ................................................................................................................... 48

  Fig. 3.6 Curve of t max versus the radius of c1 ....................................................................... 49 Fig. 3.7 Comparison between the analytical model and NIRFAST for two inhomogeneities with different absorption. Top row((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different absorption coefficient (µ a of c1 in (a) and (e) is 0.4cm -1 and 0.8 cm -1 , respectively, µ a of c2 in (a) and (e) is 0.8 cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. ................................................................................................................... 50 Fig. 3.8 Curve of t max versus absorption coefficient of c1 ................................................... 51 Fig. 3.9 Comparison between the analytical model and NIRFAST. (a)Horizontal cross section of the first numerical phantom with five cylindrical inhomogeneities (µ a =0.8cm -1 and µ' s =10cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c1 and c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the position of the inhomogeneity c5. (d) MRE corresponding to the area [-1.2cm, 1.2cm] ×[-1.2cm, 1.2cm]. .. 53 Fig. 3.10 Comparison between the analytical model and NIRFAST. (a) Horizontal cross section of the second numerical phantom with three cylindrical inhomogeneities (see Table 3.1 for optical properties) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1. (d) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (e) MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. .............................................. 54 Fig. 3.11 One update of x (0) by use of V-cycle. The grid levels are noted by (0) to (3). Red broken line represents the path which multigrid algorithm traverses. The dotted lines represents grid with different spacing. .................................................................................... 57 Fig. 3.12 Multiple updates of x (0) . The grid levels are noted by (0) to (3). Red broken line represents multigrid algorithm recursively moves forth and back. The dotted lines represents grid with different spacing. The position of each update of x (0) is noted by black points along the green line with an arrow. The number of updates is noted by 1,2,…. ............................... 57 Fig. 3.13 Illustrations of grid in finite difference. red circles represent the positions of space points. (x i ,y j ) represents coordinate values of an arbitrary point in space................................ 59 Fig. 3.14 Validation of home-developed multigrid solver of diffusion equation, (A) geometry of numerical phantom consisting of three absorber with absorption coefficient of 0.3cm -1 (larger absorber with radius of 0.5cm),0.2cm -1 (middle absorber with radius of 0.4cm) and 0.15cm -1 (small absorber with radius of 0.3cm) embedded in a homogeneous medium with absorption coefficient of 0.1cm -1 . µ s =10cm -1 . (B) Optical fluence distribution calculated with multigrid method. Profile through the optical fluence distribution calculated with

  Fig. 3.6 Curve of t max versus the radius of c1 ....................................................................... 49 Fig. 3.7 Comparison between the analytical model and NIRFAST for two inhomogeneities with different absorption. Top row((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different absorption coefficient (µ a of c1 in (a) and (e) is 0.4cm -1 and 0.8 cm -1 , respectively, µ a of c2 in (a) and (e) is 0.8 cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. ................................................................................................................... 50 Fig. 3.8 Curve of t max versus absorption coefficient of c1 ................................................... 51 Fig. 3.9 Comparison between the analytical model and NIRFAST. (a)Horizontal cross section of the first numerical phantom with five cylindrical inhomogeneities (µ a =0.8cm -1 and µ' s =10cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c1 and c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the position of the inhomogeneity c5. (d) MRE corresponding to the area [-1.2cm, 1.2cm] ×[-1.2cm, 1.2cm]. .. 53 Fig. 3.10 Comparison between the analytical model and NIRFAST. (a) Horizontal cross section of the second numerical phantom with three cylindrical inhomogeneities (see Table 3.1 for optical properties) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1. (d) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c3. (e) MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ]. .............................................. 54 Fig. 3.11 One update of x (0) by use of V-cycle. The grid levels are noted by (0) to (3). Red broken line represents the path which multigrid algorithm traverses. The dotted lines represents grid with different spacing. .................................................................................... 57 Fig. 3.12 Multiple updates of x (0) . The grid levels are noted by (0) to (3). Red broken line represents multigrid algorithm recursively moves forth and back. The dotted lines represents grid with different spacing. The position of each update of x (0) is noted by black points along the green line with an arrow. The number of updates is noted by 1,2,…. ............................... 57 Fig. 3.13 Illustrations of grid in finite difference. red circles represent the positions of space points. (x i ,y j ) represents coordinate values of an arbitrary point in space................................ 59 Fig. 3.14 Validation of home-developed multigrid solver of diffusion equation, (A) geometry of numerical phantom consisting of three absorber with absorption coefficient of 0.3cm -1 (larger absorber with radius of 0.5cm),0.2cm -1 (middle absorber with radius of 0.4cm) and 0.15cm -1 (small absorber with radius of 0.3cm) embedded in a homogeneous medium with absorption coefficient of 0.1cm -1 . µ s =10cm -1 . (B) Optical fluence distribution calculated with multigrid method. Profile through the optical fluence distribution calculated with

  ), normal finite difference scheme (black) and multigrid (red) through the line y=0.7cm (C) and x=0.7cm (D). ........................................................................................ 61 Fig. 3.15 Simulation of PA signals with pure analytical model. (A) Horizontal slice of schematic diagram of the simulation domain. Black squares represent ultrasound (US) transducers. Red circle represents point light source. Red discs represent inhomogeneities embedded in an homogeneous turbid medium. (B) Simulated PA signals recorded by the transducer. ............................................................................................................................... 63 Fig. 3.16 Simulations of PA signal with pure analytical model. (a) Horizontal slice of schematic diagram of the simulation domain. Black squares represent ultrasound (US) transducers, the position of upper, middle and lower transducer are (x=2,y=-0.5), (x=2,y=0.2) and (x=2,y=0.9) , respectively. The point light source position is (x=2, y=0). c1 c2 and c3 are inhomogeneities embedded in an homogeneous turbid medium. (B) Simulated PA signals recorded by three transducers by three US transducers. .......................................................... 63 Fig. 3.17 An example simulating PA signals with pure numerical method. Upper part represents configuration in which each pixel has a volume of 0.2mm×0.2mm×0.2mm. The absorber is spherical volume in figure and the transducer is placed on the border of medium (linear shape in figure). Lower part represents the signal recorded by one of elements in linear array. ........................................................................................................................................ 64 Fig. 4.1 Geometry and absorbed energy distribution produced with NIRFAST. (A) geometry of the numerical phantom, c1, c2 and c3 represent three inhomogeneities with different parameters (see Table 4.1). Red circle represents the position of the point light source. (B) absorbed energy distribution produced by FEM. .................................................. 68 Fig. 4.2 Segmentation of the absorbed energy distribution shown in Fig 4.1(B). (A) detected borders of the absorbers. (B) results obtained after boundary-filling method. (C) labeled regions. ........................................................................................................................ 69 Fig. 4.3 Quantitative reconstruction of absorption coefficient from noise-free absorbed energy shown in Fig. 4.1(B). (A) true absorption coefficient distribution (see Table

  ........................................................................................................... 72 Fig. 4.4 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise, signal-to-noise ratio (SNR) is ~50 dB. (A) true absorption coefficient distribution (see Table

  ............................................................ 73 Fig. 4.5 Geometry of the numerical phantom and absorbed energy distribution produced with FEM for a light source with a wide beam. (A) geometry used to produce the simulated absorbed energy, c1, c2 and c3 represent three inhomogeneities with different parameters (see

  .1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm. (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm. (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm. ................................................................................................................................................. 75 Fig. 4.7 Quantitative reconstruction of absorption coefficient from the absorbed energy that includes noise with SNR of ~25 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm. (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm. (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm. ................................................................................................................................................. 75 Fig. 4.8 Simulation tests of quantitative reconstruction of absorption coefficient by use of fixed grid and multigrid algorithm. (A) recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the line x= -0.5 cm (B) and along the line x= 0.25 cm (C). The curves of relative errors versus iteration number (D) and CPU time (E) for multigrid-based and fixed grid algorithm. ..................................... 78 Fig. 4.9 Comparison of the convergence speed of multigtrid-based algorithms with different number of grid levels ............................................................................................................... 79 Fig. 4.10 Simulation tests of quantitative reconstruction of absorption coefficient from the absorbed energy that includes noise with SNR of 25 dB. (A) recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the line x= -0.5 cm (B) and along the line x= 0.25 cm (C). The curves of relative errors versus iteration number (D) and CPU time (E) for multigrid-based and fixed grid algorithm 80 Fig. 5.1 Experimental system for PA signals acquisition .................................................... 82 Fig. 5.2 Top-view photograph of the phantom used for experimental test .......................... 83 Fig. 5.3 Recorded PA signals ............................................................................................... 83 Fig. 5.4 PA image formed by delay and sum algorithm ...................................................... 84 Fig. 5.5 Image formed by time reversal algorithm .............................................................. 85 Fig. 5.6 Segmentation of the absorbed energy distribution shown in Fig 5.5. (A) detected borders of the absorbers. (B) results obtained after boundary-filling method. (C) labeled regions. .................................................................................................................................... 86 Fig. 5.7 Quantitative reconstruction of absorption from experimental data by use of fixed grid and multigrid algorithm. (A) Recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the line x= -0.5 cm (B) and along the line x= 0.25 cm (C). The curves of relative errors versus iteration number (D) and CPU time (E) for Multigrid-based and fixed grid algorithm. ..................................... 87

Fig. 0

 0 Fig. 0.1 PAI of whole body of small animal[START_REF] Brecht | Whole-body three-dimensional optoacoustic tomography system for small animals[END_REF]

  Fig. 0.3 presents an example of PAI of breast cancer obtained by Ermilov team (Ermilov et al., 2010).

Fig. 0 . 3

 03 Fig. 0.3 An example of PAI application to breast cancer[START_REF] Ermilov | Laser optoacoustic imaging system for detection of breast cancer[END_REF]

  Fig. 0.4 gives an example of PAI of skin in human palm where the microvasculature is clearly visible. Although most studies have demonstrated the ability of PAI to visualize the skin microvasculature, the characterization of specific pathologies remains a challenge.

Fig. 0 . 4

 04 Fig. 0.4 An example of PAI of the skin (E. Z.Zhang et al., 2009)

4 . 2

 42 Numerical model ..................................................................................... 1.2.5 Discussion ........................................................................................................ solution of PA waves emitted by spherical source ......................... 1.3.2 Numerical solution of PA equation .................................................................

Fig. 1

 1 Fig. 1.1 Illustration of Bell's experiment. A represents the transmitter diaphragm. B represents the lampblack receiver with hearing tube. See in the link: http://digital.library.cornell.edu/ Therefore, a complete understanding of PA effect in biological tissues consists of a multiphysics coupling of optical illumination and acoustic propagation produced by local absorption character of the medium.The optical properties of biological tissue consist of the absorption coefficient, µ a (cm -1 ), the scattering coefficient µ s (cm -1 ) , the scattering function p(θ) (sr -1 ) where θ is the deflection angle, and the refractive index, n'(Wang & Jacques, 1992). The strong scattering of biological tissue averages the deflection angle and then an average parameter, g = <cos θ>, called the anisotropy of scatter, is used to describe tissue scattering in terms of the relative

  concentrations of two kinds of hemoglobin.

Fig. 1

 1 Fig. 1.2 Absorption spectra of different components of biological tissues.From: omlc.ogi.edu/spectra

Fig. 1

 1 Fig. 1.3 The photon travels in a random trajectory.From: http://oilab.seas.wustl.edu/

Figure 1 .

 1 Figure 1.4 is a flow chart of a Monte Carlo program(Prahl, Keijzer, Jacques, & Welch, 1989). Upon launching, the photon is moved a distance ∆s, in which the photon may be

.

  Photons are incident from the top-center of the medium.

  unknown constants B, substitute Eq. (1.10) into the condition in Eq.

  spherical (cylindrical) heterogeneity.

Fig. 1 . 6

 16 Fig. 1.6 Optical fluence at point P is the sum of incident and scattered fluence by embedded inhomogeneity C

  the vector F. (3) Solving the equation KU=F to determine the best approximation to Φ(r).

.

  The light source is placed at the point (x=1.9 cm, y=0 cm).

Fig. 1

 1 Fig. 1.7 Optical fluence distribution obtained by NIRFAST for a homogeneous turbid medium with 1 ' 1 0.1 , 10 as cm cm  

  Fig. 1.8 The PA waves detected by a point detector at a certain time t coming from sources over a spherical shell centered at the detector position and a radius of v s t Fig.1.9 shows a typical example of the PA wave from a spherical source, detected at a distance of 2 cm from a uniform spherical source with a radius of 0.2 cm. The calculation

Fig. 1

 1 Fig.1.9 An example of the pressure generated by a uniform spherical PA source Fig.1.10 shows that PA signal amplitude decays with the distance between center of sphere and detector increasing, even if the detected signal shapes have similar shape. The detected pressure amplitude is much weaker than its initial pressure amplitude. From Eq (1.12), it can be concluded that the amplitude of PA signal at the distance far from the spherical source is approximately proportional to the source size, but inversely proportional to the distance from the source object, R.Fig.1.11 shows that the signal shape varies with the size of spherical source. The temporal width of the PA signal from a spherical source is proportional to its diameter. The corresponding PA signal Fourier analysis shows a wide band spectrum. Moreover, the smaller the object is, the higher the frequency components in the generated PA signal spectrum are (see lower part of Fig.1.11). Therefore, the ultrasound detector should be selected according to the size of imaged absorber.Note that all these simulations start from a homogeneous distribution of initial pressure inside the source. In fact, it is directly related to the absorbed energy which determined by light source properties, medium properties and optical absorber properties. Therefore, the distribution of initial pressure is not uniform in practice and that could be taken into account in propagation models. This point will be discussed further in this thesis.

Fig. 1 .

 1 Fig. 1.10 An example of the pressure generated by a uniform spherical PA source placed at different distance from acoustic detector

Figure 1 .

 1 Figure 1.12 shows an example simulated by using k-wave toolbox. Left part in Fig. 1.12 shows the geometry, initial pressure distribution and ultrasound detectors (black points). Middle part shows the detected signals by ultrasound detectors. Right part shows the reconstructed image from signals shown in middle one. This simulation starts from "initial pressure", one can incorporate absorbed energy produced by incident light source to have a complete understanding real PA signals.

Fig. 2

 2 Fig. 2.1 PAM system developed by Lihong WANG's group[START_REF] Wang | Tutorial on Photoacoustic Microscopy and Computed Tomography[END_REF]

Fig. 2

 2 Fig.2.2 A representative of PAT (photoacoustic tomography) system[START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF]

  Where,  is the PA signal wavelength,s vf  with f the frequency of the PA signal, f fnumber L D 

Fig. 2

 2 Fig. 2.3 An example of curved array[START_REF] Ermilov | Laser optoacoustic imaging system for detection of breast cancer[END_REF]

Fig. 2

 2 Fig.2.5 A planar FP polymer film sensing interferometer (E. Z.Zhang et al., 2009)



  s  for biological tissue ), th  is the time scale for the heat dissipation of absorbed optical energy.The stress confinement can be approximated by is the time for the stress to transit the heated region, where v s is the speed of sound. For example, to achieve a spatial resolution d c =150 µm,

Fig. 2

 2 Fig. 2.7 Illustration of delay and sum algorithm. Transducer array (elements are noted by i, i+1,…)that has a closed geometry (evenly distributed across a circle). ∆t i is the time delay of signal acquired by i th element for point r in the medium. ∆r i is the distance between the point in question r and the position of i th elements of transducer. The equation is the relation between ∆r i and ∆t i . with v s sound speed in the medium.

Fig. 2

 2 Fig. 2.8 A linear array is positioned at z=0. x-axis represents the linear array. z-axis represents depth.The recorded signals can be expressed as p(x, 0, t).

Fig. 2

 2 Fig. 2.9 Illustration of notations to be used in TR algorithm. Transducer array has a closed geometry (evenly distributed across a circle). d i represents the position of i th element of transducer array. A represents PA wave source. Purple circle represents surface of medium, noted by S.

  is the detected acoustic pressure at the circular surface with radius= d r , () P  is the Fourier transform of ( , ) d p r t . IFFT represents inverse Fourier transform. Relative notation has been explained in Fig. 2.9.

  system should consist of closed-shape geometry of transducer array with point-like point spread function. Fig 2.10 shows an example to illustrate the difference between DAS and TR algorithm that is considered as representative reconstruction algorithms described above. Fig 2.10(a) shows the initial pressure distribution (p(r,0)

  Fig 2.10(b) and Fig 2.10(c) show the recovered results of TR and DAS algorithm, respectively. The two methods can provide the positions of wave sources (absorbers in biological tissues). Fig 2.10(d) shows horizontal cross slice through the line, x=0 mm, of initial pressure, recovered results of TR and DAS algorithm. Fig 2.10(

Fig. 2 .

 2 Fig. 2.10 Illustration the difference between DAS and TR method. (a) Initial pressure distribution and the geometry of transducer array (black points). the recovered results of TR (b) and DAS (c). (d) Horizontal cross slice through the line, x=0 mm, of initial pressure, recovered results of TR and DAS.

  model of absorbed energy, r refers to an arbitrary point in medium, ˆa  and ' ˆs  are recovered absorption and scattering coefficient map, respectively.

  analytical model of PA waves for multiple cylinders ............................. 3.3.2 Pure numerical model of PA waves ................................................................

Fig. 3

 3 Fig. 3.1 Horizontal cross section of the geometry. Cylindrical inhomogeneities of absorption coefficient embedded in an otherwise homogeneous infinite turbid medium. i is index of i th cylinder and i' represents an arbitrary remaining cylinders. Cylindrical coordinate systems have their origins at the center of each inhomogeneity. The source is noted by S. P indicates an arbitrary point within the medium, and its coordinates values with respect to each cylinder coordinate system are indicated. The z-axis is along the axis of each inhomogeneity and comes out of the cross section.

  Bessel functions of the first and the second kind, respectively.

,

  i   and i   are, respectively, the smaller and the larger values of i s  and i  , where i s  is the radial distance of the light source and i  is an arbitrary position (indicated by p) in the medium with respect to the i th inhomogeneity center. n and p arise from the separation of variables.

Fig. 3 . 2

 32 Fig. 3.2 Physical interpretation of the different orders of scattering waves by the i th inhomogeneity. Left: the incident wave produces the first-order scattered wave by the i th inhomogeneity; Right: the sum of first-order scattered waves from the remaining cylinders produces the second-order scattered wave by the i th inhomogeneity

  optical fluence distribution difference between the result calculated by the analytical model, x, and the reference result calculated by NIRFAST, y. The MRE is normalized to the local reference result y. Therefore, this parameter evaluates local errors. It is the more relevant way of evaluating the ability of the model for use in PA imaging. Global parameters would neglect deep parts of the inhomogeneities that receive lower optical fluence, and are then the more challenging part of PA imaging. Two inhomogeneities with different depth This set of simulation examples shows the interaction of two inhomogeneities with different depth. Two columns in Fig. 3.3 show two phantoms and the corresponding results. The position of inhomogeneities (denoted by c1 and c2) and light source are shown in the top rows, Fig. 3.3(a) and Fig.

Fig. 3 4

 34 Fig. 3.3 . Comparison between the analytical model and NIRFAST for two inhomogeneities with different depth. Top row ((a) and (d)): Horizontal cross section of the first numerical phantom of two cylindrical inhomogeneities with different depth (depths of the c1 center in (a) and (e) are 2.1 cm and 2. 3cm, depths of the c2 center in (a) and (e) are 1.5cm) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) and (e)): Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the positions of the inhomogeneities c1 and c2. Third row ((c) and (g)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ].

Fig. 3 . 5

 35 Fig. 3.5 Comparison between the analytical model and NIRFAST for two inhomogeneities with different radius. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different radius (radius of c1 in (a) and (e) is 0.5cm and 0.7cm, respectively, radius of c2 in (a) and (e) is 0.3cm) in a homogeneous background (µ a =0.2cm -1 and

Fig. 3 . 7

 37 Fig. 3.7 Comparison between the analytical model and NIRFAST for two inhomogeneities with different absorption. Top row ((a) and (e)): Horizontal cross section of the numerical phantom of two cylindrical inhomogeneities with different absorption coefficient (µ a of c1 in (a) and (e) is 0.4cm -1 and 0.8 cm -1 , respectively, µ a of c2 in (a) and (e) is 0.8 cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). Second row ((b) (f)): Horizontal cross section of the optical fluence distribution through y=-0.5cm; the dashed lines indicate the positions of the inhomogeneities c1. Third row ((c) (g)): Horizontal cross section of the optical fluence distribution through y=0.5cm; the dashed lines indicate the position of the inhomogeneity c2. Fourth row((d) and (h)): MRE corresponding to the area [-1.2cm, 1.2cm ] ×[-1.2cm, 1.2cm ].

Fig. 3 . 9

 39 Fig. 3.9 Comparison between the analytical model and NIRFAST. (a)Horizontal cross section of the first numerical phantom with five cylindrical inhomogeneities (µ a =0.8cm -1 and µ' s =10cm -1 ) in a homogeneous background (µ a =0.2cm -1 and µ' s =10cm -1 ). (b) Horizontal cross section of the optical fluence distribution through y=-0.7cm; the dashed lines indicate the positions of the inhomogeneities c1 and c3. (c) Horizontal cross section of the optical fluence distribution through y=0cm; the dashed lines indicate the position of the inhomogeneity c5. (d) MRE corresponding to the area [-1.2cm, 1.2cm] ×[-1.2cm, 1.2cm].

Fig. 3 .

 3 Fig. 3.11 One update of x (0) by use of V-cycle. The grid levels are noted by (0) to (3). Red broken line represents the path which multigrid algorithm traverses. The dotted lines represents grid with different spacing.

  the oscillation components in solution (wave).Therefore, one can adopt a set of multigrid as follows,1) The finest grid should meet the condition,

Fig. 3 .

 3 Fig. 3.13 Illustrations of grid in finite difference. red circles represent the positions of space points.(x i ,y j ) represents coordinate values of an arbitrary point in space.

  ) Neumann or flux boundary condition, which states a constant flux at the boundary. (3) The flux leaving the external boundary is equal to the fluence rate at the boundary weighted by a factor that accounts for the internal reflection of light back into the tissue.

  PA signals produced by multiple cylinders. Two simulation examples will be presented to show the PA signals considering the absorbed energy distribution inside objects. Figure 3.15 and Fig. 3.16 show the geometry and simulation results.

  Figure 3.15(a) shows the three absorbers at different depths.

  Figure 3.15(b)

Fig. 3 .

 3 Fig. 3.15 Simulation of PA signals with pure analytical model. (A) Horizontal slice of schematic diagram of the simulation domain. Black squares represent ultrasound (US) transducers. Red circle represents point light source. Red discs represent inhomogeneities embedded in an homogeneous turbid medium. (B) Simulated PA signals recorded by the transducer.

Fig. 3 .

 3 Fig. 3.16 Simulations of PA signal with pure analytical model. (a) Horizontal slice of schematic diagram of the simulation domain. Black squares represent ultrasound (US) transducers, the position of upper, middle and lower transducer are (x=2,y=-0.5), (x=2,y=0.2) and (x=2,y=0.9) , respectively. The point light source position is (x=2, y=0). c1 c2 and c3 are inhomogeneities embedded in an homogeneous turbid medium. (B) Simulated PA signals recorded by three transducers by three US transducers.

Fig. 3 .

 3 Fig. 3.17 An example simulating PA signals with pure numerical method. Upper part represents configuration in which each pixel has a volume of 0.2mm×0.2mm×0.2mm. The absorber is spherical volume in figure and the transducer is placed on the border of medium (linear shape in figure). Lower part represents the signal recorded by one of elements in linear array.

Fig. 4 .

 4 Fig. 4.1 Geometry and absorbed energy distribution produced with NIRFAST. (A) geometry of the numerical phantom, c1, c2 and c3 represent three inhomogeneities with different parameters (see Table 4.1). Red circle represents the position of the point light source. (B) absorbed energy distribution produced by FEM.

  the position and boundary of the absorbers by use of a 'Sobel' filter. The closed contours are then obtained by a boundary-filling algorithm.

Figure 4 .

 4 2 shows segmented results based on the absorbed energy of Fig 4.1(B). Fig. 4.2(A) shows the detected border of the absorbers, Fig. 4.2(B) shows the result obtained by the boundary-filling algorithm, Fig. 4.2(C) shows the labeled regions, in which different colors represent different regions. In the segmented results, it can be seen that the absorbers are clearly detected by our algorithm.

Fig. 4 .

 4 Fig. 4.2 Segmentation of the absorbed energy distribution shown in Fig 4.1(B). (A) detected borders of the absorbers. (B) results obtained after boundary-filling method. (C) labeled regions.

.

  represents the iteration number. The error for n th iteration is defined by, The iteration process can be stopped by required precision. In following simulation, Eq. (4.5) is used as the criterion to stop iteration.

Fig. 4 .

 4 Fig. 4.3 Quantitative reconstruction of absorption coefficient from noise-free absorbed energy shown in Fig. 4.1(B). (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm (absorber c2). (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm (absorber c3). (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm (absorber c1).

Fig. 4 . 4

 44 Fig. 4.4 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise, signal-to-noise ratio (SNR) is ~50 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm (absorber c2). (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm (absorber c3). (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm (absorber c1).

  Fig. 4.4 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise, signal-to-noise ratio (SNR) is ~50 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm (absorber c2). (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm (absorber c3). (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm (absorber c1).

Figure 4 Fig. 4 .

 44 Fig. 4.5 Geometry of the numerical phantom and absorbed energy distribution produced with NIRFAST for a light source with a wide beam. (A) geometry used to produce the simulated absorbed energy, c1, c2 and c3 represent three inhomogeneities with different parameters (see in Table 4.1). Crosses are used to approximate wide beam. (B) absorbed energy distribution produced by NIRFAST.

Fig. 4 .

 4 Fig. 4.6 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise with SNR of ~50 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm. (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm. (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm.

  Fig. 4.6 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise with SNR of ~50 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=-0.7cm. (C) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0 cm. (D) horizontal profiles through the true (black) and recovered (red) absorption distributions through the line y=0.7cm.

Fig. 4 .

 4 Fig. 4.7 Quantitative reconstruction of absorption coefficient from the absorbed energy that includes noise with SNR of ~25 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal

  Fig. 4.7 Quantitative reconstruction of absorption coefficient from the absorbed energy that includes noise with SNR of ~25 dB. (A) true absorption coefficient distribution (see Table 4.1). (B) horizontal

Fig. 4 .

 4 Fig. 4.8 Simulation tests of quantitative reconstruction of absorption coefficient by use of fixed grid and multigrid algorithm. (A) recovered absorption coefficient. True (black line) and recovered absorption coefficient profile (red and green lines) plotted along the line x= -0.5 cm (B) and along the

Figure 4 .

 4 10(A) is the reconstructed absorption map.

Figure 4 .

 4 10(B)-(D) are horizontal profiles that pass through the centers of absorbers.

  Figure 4.10(E) and (F) indicates the advantages of multigrid-based algorithm.

Fig. 4 .Fig. 4 .

 44 Fig. 4.9 Comparison of the convergence speed of multigtrid-based algorithms with different number of grid levels

Figure 5 .

 5 Figure 5.1 shows the experimental system. A pulsed light beam from a Nd:YAG laser (wavelength: 532 nm, pulse duration: 6 ns) is incident on the phantom by an optical subsystem. PA signals generated inside the phantom are acquired by a transducer (1 MHz central frequency) situated next to the phantom wall. The phantom is immersed in a water tank. A rotary stage rotates the receiver relative to the center of the phantom. The incident laser beam diameter is 5.0 cm at the phantom surface and the incident optical fluence is about 10 mJ•cm -2 .

Fig. 5

 5 Fig. 5.1 Experimental system for PA signals acquisition

Fig. 5

 5 Fig. 5.4 PA image formed by delay and sum algorithm

Fig. 5 . 5

 55 Fig. 5.5 Image formed by time reversal algorithm

  5.2) here, n j and n H are Bessel functions of the first and second kinds. Based on the same boundary conditions as those mentioned in chapter 3, one can resolve the unknown coefficient,

Fig. 5

 5 Fig. 5.6 Segmentation of the absorbed energy distribution shown in Fig 5.5. (A) detected borders of the absorbers. (B) results obtained after boundary-filling method. (C) labeled regions.

  4. Quantitative reconstruction of absorption coefficient distribution is shown in Fig. 5.7. Figure 5.7(A) -(C) show the recovered absorption map and vertical profiles that pass through the centers of targets.

  Figure 5.7(E) and (F) show the curves that represent convergence speed of multigrid and fixed grid based algorithm.

Fig. 5

 5 Fig. 5.7 Quantitative reconstruction of absorption from experimental data by use of fixed grid and multigrid algorithm. (A) Recovered absorption coefficient. True (black line) and recovered absorption

  To prove Eq. (3.4), both sides of the diffusion equation need to be expressed in cylindrical coordinates, represents the radial function. Substituting the three equations into Eq. (2) yields an equation for the radial function  

Fig. 1 A

 1 Fig.1 A typical shape of the curve of convergence factor λ(k) versus mode frequency k

  𝛷(𝑧) = 𝛷 0 𝑒 -𝜇 𝑎 𝑧 [in the ballistic regime] 𝛷(𝑧) ≈ 𝛷 0 𝑒 -𝜇 𝑒𝑓𝑓 𝑧 [in the diffusive regime] , Pedley, & Beard. In vivo high-resolution 3D photoacoustic imaging of superficial vascular anatomy. Physics in Medicine and Biology, 54(4), 1035-46 (2009). Zhang, Piao, & Bunting. Photon diffusion in a homogeneous medium bounded externally or internally by an infinitely long circular cylindrical applicator. III. Synthetic study of continuouswave photon fluence rate along unique spiral paths. Journal of the Optical Society of America A, 29(4), 2012.

Table 4

 4 .1). Crosses are used to approximate wide beam. (B) absorbed energy distribution produced by FEM. ................................................................................................................... 74 Fig. 4.6 Quantitative reconstruction of absorption coefficient from absorbed energy that includes noise with SNR of ~50 dB. (A) true absorption coefficient distribution (see Table 4

	Shengfu LI
	Photoacoustic Imaging
	INSA Lyon

  Lastly, the analytical model of optical fluence is to be incorporated into PA wave model to obtain a pure analytical model of PA signal. In chapter 4, we present our contribution to QPAT reconstruction algorithms among three directions: (1) optical forward models;

	optical fluence.
	Shengfu LI
	Photoacoustic Imaging
	INSA Lyon	5

in biological tissues

  Absorber in biological tissues ............................................................................ 1.1.2 Multi-Physics Coupling: optics and acoustics .................................................

	Chapter 1
	1 Physical models for photoacoustic effect
	Contents	PART I
	1.1 Photoacoustic effect .1.2 Optical propagation model
	part II Background
	Shengfu LI
	Photoacoustic Imaging Shengfu LI
	INSA Lyon Photoacoustic Imaging	6
	INSA Lyon

......................................................................................... 1.1.1 .......................................... 1.2.1 Statistical method: Monte Carlo ...................................................................... 1.2.2 Differential equation (1): Radiative transfer equation (RTE) .......................... 1.2.3 Differential equation (2): Diffusion equation .................................................. 1.2.4 Methods resolving the diffusion equation ....................................................... 1.2.4.1 Analytical model...................................................................................... 1.2.

Shengfu LI Photoacoustic Imaging INSA Lyon 13

Shengfu LI Photoacoustic Imaging INSA Lyon 14

Table 2 .

 2

	1 MPE for skin exposure to a laser beam at 400-1400 nm
	Wavelength λ(nm

  SA , for a given point in medium is obtained by the sum of the delayed signals of all the elements of transducer,

	e signal amplitude, p 0 N ( ) SA i p r   	( w p t t i i i  D	)
	i p is the signal that has been received by the i th transducer element, ∆t i is the time delay
	applied to that signal as presented in Fig. 2.7, e N is the total number of element. i w is the
	weighting factor.
	Shengfu LI
	Photoacoustic Imaging
	INSA Lyon		32

  .1.1 Geometry ......................................................................................................... 3.1.2 Incident model ................................................................................................. 3.1.3 Scattering model .............................................................................................. 3.1.4 Boundary conditions ........................................................................................ 3.1.5 Validations .......................................................................................................

	Chapter 3
	3 PA forward model
	Contents	PART II
	3.1	Analytical optical models in PAI .
		part II Contributions
	Shengfu LI Shengfu LI
	Photoacoustic Imaging Photoacoustic Imaging
	INSA Lyon INSA Lyon	37

..................................................................

33.2 Multigrid solver of the diffusion equation .................................................... 3.2.1 Multigrid idea briefly sketched ........................................................................ 3.2.2 Multigrid solution of diffusion equation ......................................................... 3.2.3 Implementations of Multigrid solver ............................................................... 3.2.4 Validation ........................................................................................................

Table 3 .1 Parameters of the second numerical phantom

 3

	Inhomogeneity	Position (x, y)	Radius (cm)	Absorption (cm -1 )
	c1	(0.5, 0)	0.5	0.8
	c2	(-0.5, 0.7)	0.4	0.6
	c3	(0, -0.7)	0.3	0.4

  .13. Red circles represent the positions of space points and (x

	Shengfu LI
	Photoacoustic Imaging
	INSA Lyon	59

i ,y j ) represents coordinate values of an arbitrary point in space. For convenience, (x i ,y j ) is replaced with (i,j). Then the optical fluence at the point (x i ,y j ) is noted , ij  .

Table 4 .1 Parameter of the numerical phantom

 4

	Inhomogeneity	Center position (x,y)	Radius(cm)	Absorption(cm -1 )
	C1	(0.7, 0.7)	0.5	0.3
	C2	(0.7, -0.7)	0.4	0.2
	C3	(-0.7,0)	0.3	0.15

Table 4 .2 Reconstruction results obtained by segmented results with the proposed algorithm

 4

	Shengfu LI
	Photoacoustic Imaging
	INSA Lyon

Table 4 .3 Reconstruction results obtained by assuming that the position of absorbers is known exactly

 4

	Inhomogeneity	Absorption(cm -1 )
	C1	0.3001
	C2	0.2
	C3	0.1501

  Calibration ....................................................................................................... 5.4.2 Analytic-based reconstruction ......................................................................... 5.4.3 Multigrid-based algorithm ...............................................................................

	Chapter 5 5 Experiment validation Contents 5.1 Experiment system ......................................................................................... 5.2 Phantom and Recorded PA signals ............................................................... 5.3 Absorbed energy distribution ........................................................................ 5.4 Quantitative PA reconstruction .................................................................... 5.4.1 Chapter 5 Experiment validation
	Shengfu LI Shengfu LI
	Photoacoustic Imaging Photoacoustic Imaging
	INSA Lyon INSA Lyon	82

5.5 Conclusion .......................................................................................................

Table 5

 5

	.1. The

Table 5 .2 Reconstructed and true absorption

 5

	Inhomogeneity	Reconstructed (cm -1 )	True (cm -1 )
	C3(pure absorption)	0.17	0.2
	C2 (absorption & scattering)	0.21	0.2

  MG(H) represents the modeled data obtained by multigrid solver, J represents Jacobian matrix. Here, it can be rewritten as, is usually defined as the ratio of the variances of the measurement data and optical properties. I is the identity matrix.

													(	 I X J MG H H ) ( ) TT J J   D  	m	
			  		10 log	  MG H 1 ( ) 1 a			10 ln log  MG H 2 ( ) 1 a		  10 log MG H ( ) 1 NI a		;			10 log  MG H 1 ( ) 1 s				10 log	  MG H 2 ( ) 1 s			  10 log MG H ( ) 1 NI s		  
	J		   		10 log	  MG H 1 ( ) 2 a			10 log	  MG H 2 ( ) 2 a			  10 log MG H ( ) 2 NI a		;		10 log	  MG H 1 ( ) 2 s				10 log	  MG H 2 ( ) 2 s			10 log	( MG H NI   s	2	)		   
																															
			   		10 log	  MG H 1 ( ) N a			10 log	  MG H 2 ( ) N a			10 log	  MG H ( ) N NI a		;		10 log	  MG H 1 ( ) N s			log	10	  MG H 2 ( ) N s			log	10	  MG H ( ) N NI s		   
	1   ,  2 1 2 ,..., , , ,..., NI NI     D  D D X D D D D	T
						a	a					a	s		s		s
		Shengfu LI
		Photoacoustic Imaging
		INSA Lyon																									98

 . NI represents the number of region, J represents Jacobian matrix. The Jacobian defines the relationship between changes in absorbed energy H  , resulting from small changes in optical properties   , as X   . λ is the regularization parameter, which

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2015ISAL0021/these.pdf © [S. Li], [2015], INSA de Lyon, tous droits réservés

L'imagerie photoacoustique (IPA) des tissus biologiques permet de combiner les avantages des imageries optique et ultrasonore. Le principal contraste endogène pour l'IPA provient des vaisseaux sanguins en raison de la forte absorption de l'hémoglobine par rapport aux tissus environnants.

De plus, les vaisseaux sanguins sont à peu près cylindriques et la concentration d'hémoglobine peut être supposée uniforme à l'intérieur des veines. Comme première contribution, nous avons développé dans cette thèse un modèle analytique de fluence optique pour plusieurs inhomogénéités cylindriques parallèles incorporées dans un milieu turbide. Les modèles n'existent que pour les cas simples. Pour traiter des situations plus complexes, comme les tissus biologiques, les méthodes numériques sont nécessaires. La deuxième contribution de cette thèse consiste à développer un solveur multigrilles de l'équation de diffusion optique et donc de proposer une méthode numérique efficace pour résoudre la fluence optique.

Enfin, notre troisième contribution concerne la reconstruction de la tomographie quantitative photoacoustique (TQPA). Basée sur les modèles efficaces présentées dans les première et seconde contributions, nous avons proposé une méthode de reconstruction basée sur le modèle direct analytique pour les cas simples et une méthode d'inversion basée sur multigrille pour les cas plus réalistes. Les avantages de la méthode d'inversion basée sur multigrille sont présentés à la fois en terme de temps de calcul et de vitesse de convergence. Une validation expérimentale est présentée dans le dernier chapitre de cette thèse, prouvant la validité et l'analyse des performances des méthodes développées. MOTS-CLÉS : Imagerie photoacoustique, imagerie biomédicale, fluence optique, modèle analytique, modèle numérique, problème inverse, méthodes multigrilles, reconstruction quantitative, tomographie.