Introduction

The present work was motivated by the development of alternative imaging techniques for breast cancer early diagnosis. Among various imaging techniques in the field of breast examination, photoacoustic imaging is gaining increasing importance, because it potentially couples the merits of optical imaging and ultrasound imaging that is high optical functional contrasts brought by optical probing and high spatial resolution by ultrasound detection. Since 1990s, much progress was made in the application of photoacoustic in biomedicine. Several photoacoustic tomography systems have been developed for breast cancer diagnosis showing breast tumours can be detected deep inside the tissues. However, they do only provide global visual information. A data acquisition protocol coupled with a robust inversion algorithm, able to provide truly quantitative and interpretable parameters for the physician, was still to be developed.

Personal achievement

Within this context, the present work seeks to model properly the photoacoustic multiwave phenomenon and incorporate it in an efficient reconstruction algorithm to solve the inverse problem.

The inverse problem consists in the recovery of 3D interior maps of physical properties that have to be able to monitor quantitatively breast cancer activity.

It is also part of the AVENTURES (Analyse Variationnelle EN Tomographies photoacoustique, thermoacoustique et UltrasonoRES, AVENTURES-ANR-12-BLAN-BS01-0001-01) project financially supported by ANR (Agence Nationale de la Recherche), the French National Research Agency. This project is both theoretical and experimental, and uses the joint models (acoustics, optics and mathematics), for numerical analysis and experiments. Mathematically, the models in Ultrasonic Tomography (UT), the Thermo-Acoustic Tomography (TAT) and Photo-Acoustic Tomography (PAT) are the same and these techniques are often studied similarly. However, the models can be refined to specify the source term. Indeed, the emitted signals are different and depend on the wavelength of the excitation signal. The purpose of this project is to establish models for photoacoustic tomography (PAT) and ultrasound tomography (UT), which take into account the micro and mesoscopic optical effects, develop various methods for satisfactory image reconstruction (inverse problems) techniques and compare them with the traditional methods, and generate effective codes both for the direct and the inverse problems, then conduct experiments and validate models in UT, TAT and PAT. Finally optimize the experimental bench. This project is cooperated by four universities or research institutes in France: MAPMO Laboratory from university d'Orléans, ISAE in Toulouse, Laboratory of Mechanics and Acoustics (LMA) and Institut Fresnel (IF) in Marseille. The work of our laboratory in this project is to model the photoacoustic phenomenon, and develop the experimental bench and conduct the PAT experiments.

In the present work, models were built for both the forward and the inverse problems for PAT.

A multiphysics forward solver has been firstly implemented coupling light and sound propagations. The propagation time units being extremely different resolution methods have been adopted for solving the governing partial differential equations (PDEs), resulting in a hybrid forward solver software: the classical Finite Element Method (FEM) was chosen for solving the optical propagation equation, while a semi-analytical method based on Fourier transforms calculations (kspace method)was preferred for solving the acoustic propagation equation. Several techniques have been studied to solve inverse problem. In a direct, non-quantitative approach, classical techniques, inherited from ultrasound probing have been adopted: from the acoustic pressure measurements, the problem can be seen as an acoustic sources localization problem, and the Time Reversal Method (TRM) can be used for reconstructing initial pressure distribution maps. We call this method of probing 'passive' because it is analogous to the problem of detecting earthquakes in seismology. The forward solver allows generating synthetic data that can be eventually blurred by introducing noise.

A preliminary study has been conducted on the optimization of the detectors nature and positions as a function of the size of the objects to be reconstructed (tumours).

To achieve a real quantification, an 'active' probing is introduced and tested: the idea is to interrogate optically, locally, at different positions, the sample. By doing so, the problem becomes more analogous to echography, localization of cracks in concrete, and has some similarities with Diffuse Optical Tomography. This multiple optical sources probing conducted to the proposition of a new experimental protocol for the PA signals acquisition with the experimental bench developed at IF. In the last step, experimental data were acquired and compared to simulations.

Presentation of the manuscript

Chapter 1 studied the principle of breast cancer, the conventional techniques for breast cancer diagnosis, the application of photoacoustic imaging in biomedicine, and the existing photoacoustic imaging systems for breast cancer.

Chapter 2 presents the governing equations in the PA phenomenon. First, the main optical contrasts of the biological tissue are defined. The Radiative Transport Equation (RTE), describing the light propagation in tissues is described. The thickness of the samples is large compared to the transport mean free path, and the Diffusion Approximation (DA) to the RTE holds: the optical propagation model is reduced to a diffusion equation. The photoacoustic phenomenon is then described and the photoacoustic equation established. Different methods of resolution can be used for solving the forward and inverse models. They are presented and discussed.

Chapter 3 describes the multi-physics model implementations. The resolutions of the forward models (optical and acoustic phenomena) and adopted inverse models are presented and tested.

Chapter 4 presents the two optimization studies in PAT conducted within this thesis:

-the optimization of the acoustic detectors distribution was firstly studied, we considered the size of one single detector, the number of the detectors, the coverage of the detectors (full or sparse angular coverage), and the limited frequency bandwidth of the detectors both in 2D and 3D. -the optimization of the illumination/detection protocol in Quantitative PAT (QPAT). It is shown that a quantification is not possible with a single measurement; the robustness of the approach was tested under different potential experimental situations with single wavelength illumination, point scanning single element transducer and limited view angle examination. The study reveals that the method, although promising, is highly sensitive to noise.

Chapter 5 presents the experimental setup and trials of comparison between measurements and simulations are presented. The proposed protocol was tested, however the experimental setup needs additional improvements to meet the reconstruction method requirements.

Chapter 1. Introduction to photoacoustic imaging

. Data for breast cancer

Breast cancer is the leading cancer in females, with nearly 1.7 million new cases diagnosed in 2012 in the world (second most common cancer overall). This represents about 12% of all new cancer cases and 25% of all cancers in women (Statistics from World Cancer Research Fund International). Table1.1 lists the countries with the top 20 highest incidence of breast cancer in 2012.

France ranked the top 3 of all the countries with the highest incidence of breast cancer, and

European countries usually have higher rate of breast cancer. One in a thousand people have the chance of breast cancer on the average. There are 232,670 new cases of breast cancer in female and 2,360 in male in the United States in 2014, and there are 40,000 cases of death in female and 430 in male.

Importance of breast cancer early diagnosis

Thus, breast cancer has high rate of incidences. Breast cancer is categorized into different stages: Stage I, II (A or B), III (A, B, or C), or IV (National Cancer Institute (NCI), America). The classification of the stages is based on the size of the tumour and whether the cancer has spread.

Stages I, IIA, IIB, and IIIA are considered as "early-stage" breast cancer and refer to as cancers that may have spread inside the breast but not to other parts of the body. People have more chance of recovery if the breast tumours are removed before they spread, so early diagnosis of breast cancer is of vital importance. Early stage of the breast cancer means smaller size of the tumour (in Stage I, the tumour size < 2 cm, while in Stage II, the tumour size < 5 cm), it's also a challenge to the techniques of breast cancer diagnosis. Physicians and scientists are focusing on the development of better technologies for breast cancer early detection.

Besides the improvement of the techniques for breast cancer early diagnosis, people who have higher risk of breast cancer, like women who are over 40 years old, people who have family history of breast cancer, and people who have personal history of breast cancer should undertake regular breast examination even if no symptom of breast cancer.

The french programs organized since 1987 by CNAM (Caisse Nationale d'Assurance Maladie) comprises palpation, mammography and biopsy. A statistical analysis published in 2010 (source:

INVS, CépiDC Inserm, 2010) conducted over 10 years shows: i) the number of new cases increases; ii) mortality decreases over all age categories (especially for young and middle aged women) except for women >85 years old. These statistics highlight two points:

-Efforts must be followed on the development of imaging techniques for early stage and small size cancers diagnosis; -Improve available therapies efficacy at the detection step.

Physics of breast cancer

The breast is the tissue overlying the chest muscles. The basic units of the breast are the lobulars, which produce the fatty breast milk. They give the breast its offspring-feeding functions as a mammary gland, the milk travels through a network of tiny tubes called ducts. The ducts connect and come together into larger ducts, which eventually exit the skin in the nipple. The dark area of skin surrounding the nipple is called the areola.

Breast cancer is a type of cancer originating from breast tissue, most commonly from the ducts or the lobules. Cancers originating from ducts are known as ductal carcinomas, while those originating from lobules are known as lobular carcinomas.

In addition, there are two types of breast tumours that are not cancer but increase the risk of breast cancer (Cancer Prevention and Treatment Fund, America): http://en.wikipedia.org/wiki/Mammary_gland Ductal carcinoma in situ (DCIS) -DCIS is a condition in which abnormal cells are found in the lining of breast ducts. These cells have not spread outside the duct to the surrounding breast tissue. But some cases of DCIS become breast cancer over time. So DCIS is sometimes called Stage 0 breast cancer. Since it's not possible to know which cases of DCIS will become breast cancer, it's important to get treatment for DCIS. Women with DCIS often are treated with breast-sparing surgery and radiation therapy. Radiation therapy lowers the chance that DCIS will come back or develop into breast cancer. If a large area of DCIS is found or it is found in more than one location, some women will choose to have a mastectomy. Underarm lymph nodes usually are not removed in the treatment of DCIS. The drug tamoxifen, which stops the growth of breast tumours that depend on estrogen, is also sometimes used in the treatment of DCIS. Tamoxifen may decrease the risk of a breast cancer developing in the same breast after treatment or in the opposite breast.

Lobular carcinoma in situ (LCIS) -LCIS is a condition in which abnormal cells are found in breast lobules. In contrast with DCIS, there is no evidence that the abnormal cells will become cancerous. However, having LCIS means that a woman has an increased risk of developing breast cancer in either breast. Despite this increased risk, most women with LCIS will never get breast cancer. Most women with LCIS are followed closely with regular checkups and mammograms. Some women choose to take tamoxifen to decrease their risk of developing breast cancer. Rarely, women with LCIS choose to have both breasts removed as a preventive measure, but most doctors think this approach is inappropriate. Almost all malignant tumours show the formation of new blood vessels from pre-existing vessels, this phenomenon is called angiogenesis. Approximately 95% of all human cancers are carcinomas that originate in situ lesion which is not usually angiogenic at the outset. In situ is not yet a cancerous stage and the tumour can stay in this phase from months to years. Switching to a malignant phenotype is followed by neovascularisation around the tumour. Usually, only when neovascularization has been developed, the tumour shows differences from normal surrounding tissues, new blood vessels generate and then grow and invade surrounding tissue using this extra blood supply.

As blood is one of the main light absorber in the tissue, in the visible wavelength range, formation of extra blood vessels means higher light absorption in tumours than surrounding normal tissues, besides, the malignant breast tumour also has higher dielectric conductivity [START_REF] Lazebnik | A large-scale study of the ultrawideband microwave dielectric properties of normal breast tissue obtained from reduction surgeries[END_REF] due to the increase of the water content, water is also an absorber in tissues in the infrared.

Conventional techniques for breast cancer diagnosis

Several techniques have developed focusing on breast tumour, seeking for early stage detection. All of them have their advantages and drawbacks and bring concurrent information in terms of morphology or physiology: X-ray mammography uses the X-ray to scan the breast, it produces images with very good resolution. However, exposure to X-ray radiation is not risk-free; MRI (Magnetic resonance imaging) provides strong anatomical contrast between soft tissues, but it is relatively expensive and time consuming; MBI (Molecular breast imaging) uses the nuclear medicine scanner to scan the breast, it is sensitive to breast tumour, but the resolution is poor compared with other techniques, like X-ray , it also suffers from longer detection time; ultrasound has high penetration ability, but no contrast agents (based on nano-or micro-bubbles) are available yet, it is more suitable to detect tissue boundaries; Optical imaging benefits high optical contrast between normal tissue and breast tumour, and it is relatively inexpensive and risk-free, but light propagating in tissue is highly scattered, resulting in poor resolution. Additional details on these techniques are provided hereafter.

X-ray Mammography

X-ray mammography is the most commonly used method in early detection of breast cancers. It uses low dose X-ray and high contrast, high-resolution film to do the examination. In X-ray mammography, the breast is placed at an X-ray machine and firmly compressed with a plate. Two Xrays images are taken of each breast at different angles.

The x-ray source is turned on and x-rays are radiated through the compressed breast and onto a film cassette positioned under the breast, thus exposing the film with an image of the internal structures of the breast. The film then treated in a dark room. Nowadays digital X-rays detectors are generally used instead [START_REF] Bansal | Digital radiography. A comparison with modern conventional imaging[END_REF]. http://fitsweb.uchc.edu/student/selectives/Luzietti/ Breast_mammography.htm Breast compression is necessary in order to:

-Reduce the breast thickness so that all of the tissue can be visualized.

-Spread out the tissue so that small abnormalities are less likely to be obscured by overlying breast tissue. -Allow the use of a lower X-ray dose since a thinner amount of breast tissue is being imaged.

-Hold the breast still in order to minimize blurring of the image caused by motion.

-Reduce x-ray scatter to increase sharpness of picture.

As the X-rays pass through the breast, they are attenuated by different types of tissues. Fat is very dense and absorbs a great amount of the x-rays. The tissue around the breast ducts and fat is less dense and absorbs far less x-ray energy. It is these differences in absorption and the corresponding varying exposure level of the film that create the images which can clearly show different tissues, such as fat, breast ducts, and so on. The breast is made of fat, fibrous tissue and glands. Breast masses (these include benign and cancerous lesions) appear as white regions on mammogram film. Fat appears as black regions on a mammogram film. Everything else (glands, tumours and other significant abnormalities) appear as levels of white on a mammogram.

X-ray mammography can detect small abnormal tissue (0.1 mm) growths in the breast. These early tumours do not harm patients if they are removed at this stage.

Breast MRI (Magnetic resonance imaging)

MRI is a non-invasive medical test, it uses a powerful magnetic field, radio frequency pulses and a computer to produce detailed pictures of organs, soft tissues, bone and virtually all other internal body structures. The magnetic field is produced by passing an electric current through wire coils in most MRI units. Other coils, located in the machine and in some cases, placed around the part of the body being imaged, send and receive radio waves, producing signals that are detected by the coils. A computer then processes the signals and generates a series of images, each of which shows a thin slice of the body. The images can then be studied from different angles by the interpreting radiologist.

To do the MRI exam, the patient should lie on a special table inside the MRI system where a radio signal is turned on and off. It consists in probing the nuclear magnetic resonance of water protons present in the body, hence it probes essentially the water concentration. The absorbed energy is echoed or reflected back out of the body. These echoes are continuously measured by the MRI scanner. A digital computer reconstructs these echoes into images of the breast. To highlight the abnormalities, contrast-enhancing material called Gadolinium DTPA is injected into the patient's bloodstream before or during the exam to improve the quality of the images. MRI has many advantages: MRI allows images to be taken in any plane and from any orientation, and it's sensitive to small abnormalities. MRI can also image the breast implants; MRI may also be useful in determining whether breast cancer has spread into the chest wall or not. MRI also has its drawbacks: First, MRI cannot distinguish between cancerous and noncancerous abnormalities efficiently. Furthermore, MRI is an expensive detecting method, it also takes much longer time: about 30-60 minutes compared to 10 to 15 minutes for mammography. MRI is not a replacement for X-ray mammography or ultrasound imaging but rather a supplemental tool.

Molecular breast imaging (MBI)

Molecular breast imaging (MBI) is a new nuclear medicine technique that utilizes small γ-cameras in a mammographic configuration to provide high-resolution functional images of the breast.

In MBI, a radioactive tracer (Tc-99m Sestamibi) is injected into the body to lights up any areas of tumour inside the breast.

Breast tumour cells tend to take up much more radioactive substance than normal cells do. Breast is lightly compressed between two opposing cadmium zinc telluride (CZT) detectors; a nuclear medicine scanner then scans the breast, two 10-min images are adopted of each breast, thus the total imaging time is 40 min.

MBI is sensitive to small breast tumours, and the patients suffer from less pain because of the lighter compression force. The disadvantage is that the imaging time is a little bit longer.

Ultrasound imaging

Ultrasound imaging, also called ultrasound scanning or sonography, involves the use of a small transducer and ultrasound gel to expose the body to highfrequency sound waves. When a sound wave strikes an object, it bounces back, or echoes. By measuring these echo waves, it is possible to determine how far away the object is and its size, shape and consistency. In an ultrasound examination, a transducer both sends the sound waves and receives the echoing waves.

In ultrasound imaging, the patient lies on an examination table, a water-based gel is applied to the area of the body being studied. The gel will help the transducer contact with the body and http://www.auntminnie.com/index.aspx?s ec=spt&sub=mbi&pag=dis&itemID=91850 http://www.wishmd.com/services/ultrasound.htm l eliminate air between the transducer and the skin. The transducer is placed on the body and moved back and forth over the area of interest until the desired images are captured.

Ultrasound imaging is non-invasive and does not use any ionizing radiation, but is more suitable to detect tissue boundaries and provides high quality morphological information. Though, contrast agents for ultrasounds based on nano-or micro-bubbles are being developed and functionalized [START_REF] Zheng | Nano/microparticles and ultrasound contrast agents[END_REF]. The use of contrast agents in ultrasound diagnosis may bring additional quantitative information related to tumours activity.

Optical mammography

Optical mammography uses near-infrared light (NIR) to detect breast tumour. The hemoglobin concentration and its oxygen saturation are two key parameters that may allow not only the detection of cancer, but also the discrimination between malignant and benign breast lesions using optical methods. The hemoglobin saturation values are much lower in tumour than that in normal tissue.

Optical mammography takes advantage of high optical contrast between blood vessels and bloodrich areas (such as the tumour) in the breast, as the tumours are always associated with increased vascularisation.

Optical imaging uses light propagation models to reconstruct the optical properties of the illuminated tissue. By varying the wavelength of the optical source the spectroscopic dependence of optical properties can be obtained. Optical imaging is relatively inexpensive and risk-free. However, light propagating in tissue is highly scattered, resulting in poor resolution.

Photoacoustic imaging combines the merits of optical imaging and ultrasound imaging, optical contrast on the one hand, and low scattering of ultrasound in breast tissue on the other and could potentially provide quantitative images with higher resolution than optical imaging.

Photoacoustic imaging techniques in biomedicine

The photoacoustic (PA) effect, also known as thermoacoustic effect or photoacoustic effect, describes the conversion of electromagnetic radiative energy into acoustic energy in a medium, e.g. http://medicalphysicsweb.org/cws/article/research/27191 the generation of sound from the absorption of light [START_REF] Bialkowski | Photothermal Spectroscopy Methods for Chemical Analysis[END_REF]. The effect was firstly discovered by [START_REF] Graham | On the production and reproduction of sound by light[END_REF] in 1880 in his search for wireless communication while working on an invention he called the photophone.

The photophone comprised a transmitter consisting of a voice-modulated mirror which reflected a beam of sunlight towards a receiver with a selenium cell connected to a loudspeaker. The invention allowed wireless communication over a distance of 400m but Bell realized that the range would be difficult to extend due to the lack of appropriate light sources. During his experiments however Bell noted that an audible sound was emitted from the selenium cell itself.

However, research on the photoacoustic effect made little progress primarily due to the lack of appropriate light sources. It was not until 1970s, after the invention of lasers, that photoacoustic regained interest, when Kruezer (L. B. [START_REF] Kreuzer | ultralow gas concentration infrared absorption spectroscopy[END_REF]) detected gas constituents using photoacoustic effect by a laser -induced PA effect.

Although the applications of the PA effect in biomedicine began in the 1970s, not much progress has made in photoacoustic until 1990s. Kruger [START_REF] Kruger | Photoacoustic ultrasound[END_REF] introduced the idea of using time-resolved photoacoustic signals for imaging the optical absorption in optically diffusive media. Karabutov (A.A. Karabutov et al., 1996) developed the theory of photoacoustic signal generation and ultrasonic wave propagation in liquid media. Oraevsky (A. A. Oraevsky, S. L. Jacques, R. O. Esenaliev, 1996) proposed using time resolved photoacoustic signals in analogue to echographic imaging as a medical diagnostic imaging tool in a US patent application in 1996. Oraevsky et al [START_REF] Oraevsky | Measurement of tissue opticalal properties by time-resolved detection of laser-induced transient stress[END_REF] demonstrated that it was possible to use the temporal profile of photoacoustic signals to obtain estimates of the average optical properties of different biological tissues that were of the correct order of magnitude. Wang [START_REF] Wang | Microwave-induced acoustic imaging of biological tissues[END_REF] presented thermoacoustic images using microwave excitation of a sample. In their work, pulsed microwave radiation was used to illuminate the samples. Absorbed microwave energy caused thermoelastic expansion that radiated acoustic waves. The imaging contrast is based on the difference in the dielectric constants among biological tissues. Because of the large contrast in microwave absorption among different tissue types, microwave-induced acoustic tomography could potentially provide a new modality for detecting early-stage cancers.

Generation of photoacoustic waves

When a pulsed laser illuminates the optically absorbing objects in the turbid medium, the objects will absorb some amount of energy, the absorbed energy converts into heat and the temperature of the objects increase, if the thermal confinement is met, the optically absorbing objects will emit the absorbed optical energy as a transient photoacoustic pressure wave.

Information related to the location, size, shape, and optical properties of the absorbing objects is present in the measured pressure waves. The pressure signals acquired over a number of transducers can be processed by a reconstruction algorithm and initial pressure distribution maps of the optically absorbing objects can be produced. The photoacoustic imaging combines the merits of optical and acoustic imaging, high optical contrast brought by optical probing and high spatial resolution by acoustic detection.

Contemporary Photoacoustic Imaging Modalities and Applications

Enormous applications of photoacoustic effect have been studied in the 21st century. Oraevsky and Karabutov (Vo-Dinh, 2003), [START_REF] Xu | Photoacoustic imaging in biomedicine[END_REF] have introduced some basic principles, imaging modalities and reconstruction methods. There are mainly two photoacoustic imaging modalities, related to different probing depth:

-photoacoustic microscopy, that seeks to achieve deeper probing than conventional microscopy, with moreover a quantitative information and higher resolution; -photoacoustic tomography, that aims to provide the same quantitative information as Diffuse Optical Tomography but with a higher resolution because ultrasounds are unscattered.

As PAT relies on the measured photoacoustic pressure and reconstruction methods, both can introduce some image artifacts. Despite of this defect, large efforts are invested in developing PAT for various applications in biomedicine.

Breast cancer imaging is one of the main applications of PAT: the two PAT breast imaging systems, described in the next section, are ongoing through clinical trials. Kruger [START_REF] Kruger | Breast cancer in vivo: contrast enhancement with thermoacoustic CT at 434 MHz-feasibility study[END_REF] also developed a PAT system for breast cancer diagnosis, In their study, they produced three-dimensional maps of tissue absorption at 434 MHz with spatial resolution of 1.5-4.0 mm, they also performed a limited clinical trial.

PAT can be applied to molecular imaging with the injection of optical absorbent biomarkers or targeted molecular probes to enhance the absorption contrast. Sentinel lymph node (SLN) was detected by using methyline blue as the contrast agent as shown in Fig. 1-1. [START_REF] Li | Photoacoustic imaging of lacZ gene expression in vivo[END_REF] imaged the expression of lacZ gene by using a chromogenic assay. [START_REF] De La Zerda | Carbon nanotubes as photoacoustic molecular imaging agents in living mice[END_REF] used bioconjugated carbon nanotubes to target mice tumour. PAT has the potential to visualize dynamic and functional properties of the nervous system.

The functional representations of whisker stimulations in the cerebral cortex of rats was visualized in [START_REF] Wang | Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain[END_REF].In their work, they also imaged hyperoxia and hypoxia-induced cerebral hemodynamic changes.

Besides the above mentioned applications, PAT can be also used as temperature monitoring [START_REF] Larina | Real-time photoacoustic monitoring of temperature in tissues[END_REF], biological cell detection [START_REF] Weight | Photoacoustic detection of metastatic melanoma cells in the human circulatory system[END_REF].

Hereafter, a review on the different PA imaging techniques is proposed.

PAM (Photoacoustic microscopy)

Photoacoustic microscopy uses a focused ultrasonic transducer, detects photoacoustic waves from the focal zone. The laser beam is focused into a rather small volume, and no scattering occurs in PAM, so photoacoustic signal can be used directly for the formation of the photoacoustic image, the axial resolution is derived from the time resolved detection, the lateral resolution can be determined by acoustic resolution (AR-PAM) or optical-resolution (OR-PAM) (Fig. 1-2). In OR-PAM, the laser beam is tightly focused to a diffraction-limited spot, the lateral resolution is determined by the optical focal spot size [START_REF] Hu | Second-generation opticalal-resolution photoacoustic microscopy with improved sensitivity and speed[END_REF] as shown in Fig. 1-2 (a); While in AR-PAM, the laser beam is loosely focused to fulfilling the entire acoustic detection volume shown in Fig. 1-2 (b), the lateral resolution does not depend on the tissue's optical scattering properties, but the ultrasonic focusing.

The resolution of the PAM is limited by the duration of the light pulse and the absorption of high-frequency ultrasound in tissue, but the high-frequency ultrasound attenuates severely in the tissue. Besides, the illumination may suffer from strong acoustic reverberations due to strong absorption close to the surface, this photoacoustic signals generated by these strong absorbers can reverberate inside the detection system, and overshadow late arriving signals. [START_REF] Maslov | In vivo dark-field reflection-mode photoacoustic microscopy[END_REF]introduced dark-field confocal AR-PAM that can help relieve this problem. The acoustic lens provides high lateral resolution of 45-120 μm in this system. An ultrasonic detection system provides high axial resolution of about 15 μm. The system is capable of imaging optical-absorption contrast as deep as 3 mm in biological tissue. The image of the vascular distribution in rat skin obtained by this system is shown in Fig. 123. ( [START_REF] Guittet | High-frequency estimation of the ultrasonic attenuation coefficient slope obtained in human skin: simulation and in vivo results[END_REF] studied high-frequency estimation of the ultrasonic attenuation coefficient slope obtained in human skin performed at 40 MHz, the attenuation coefficient slope ranging from 0.8 to 3.6 dB/cm MHz. Instead of using high frequency transducers, OR-PAM uses optical focusing to provide the lateral resolution. [START_REF] Maslov | Opticalal-resolution photoacoustic microscopy for in vivo imaging of single capillaries[END_REF] developed the first OR-PAM system, which provides a lateral resolution of 5 μm and an imaging depth of 0.7 mm. Fig. 1234shows an image of microvasculature in a nude mouse ear obtained by this system. However, the high optical scattering in the skin limits the imaging depth of OR-PAM.

Photoacoustic spectroscopy

By using multiple wavelengths, Photoacoustic spectroscopy can reveal more functional information about the tissue. A photoacoustic spectrum of a sample can be recorded by measuring the sound at different wavelengths. This spectrum can be used to identify the absorbing components of the sample. Photoacoustic spectroscopy can detect the spectrally dependent absorption parameters in tissue. Haemoglobin is the primary absorber in the blood. There are two kinds of haemoglobin: oxygenated haemoglobin (HbO 2 ) and deoxygenated haemoglobin (Hb), the absorption coefficient is related to HbO 2 and Hb:

)] ]( )[ ( ) ]( )[ ( [ 303 . 2 ) , (

HbO

are molar concentrations of two kinds of hemoglobin (mol L -1 ). And the haemoglobin oxygen saturation (SO2), is defined as:

] [ ] [ ] [ 2 2 2 HbO Hb HbO SO + = (1. 2)
If the absorption coefficient is measured via a PA measurement and the Hb and 2 HbO molar extinction coefficients are known, the molar concentrations of two kinds of haemoglobin can be obtained only when multiple wavelength illumination (at least 2) is adopted. Knowing the molar concentrations of two kinds of haemoglobin, the SO2 can be known also. Fig. 1-5 shows images of SO 2 in single subcutaneous blood vessels in a rat.

PAT (PhotoAcoustic Tomography)

PAT relies on a computed image reconstruction: the photoacoustic signals generated by the sample are collected by the acoustic transducers positioned all around the sample. By using the collected signals and reconstruction algorithms, the image of the physical properties of the sample can be reconstructed.

The acoustic transducers and the reconstruction algorithms play an important role in PAT.

Various types of acoustic detectors can be used: point or small surface area detectors, line or large plane detectors acoustic detector arrays. The single-element finite size flat piezoelectric transducer, is widely used in PAT [START_REF] Wang | Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain[END_REF], besides transducers made of piezoelectric materials, optical methods can be used to provide high sensitivity acoustic detection. In [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], an imaging system based on a planar Fabry-Perot film sensing interferometer (FPI) was developed, this kind of transducer can provide a practicable alternative to piezoelectric-based photoacoustic systems for high-resolution structural and functional imaging of the skin microvasculature and other superficial structures.

In addition to transducers that are much smaller than the scanning region, large planar and line detectors have also been used. [START_REF] Haltmeier | Thermoacoustic computed tomography with large planar receivers[END_REF]) used a piezoelectric plate detector and a line optical fiber-based Fabry-Perot interferometer was used as transducer in [START_REF] Paltauf | Photoacoustic tomography using a Mach-Zehnder interferometer as an acoustic line detector[END_REF].

Acoustic transducer arrays are also used or specifically developed, the array in a commercial ultrasound imaging system was used in [START_REF] Witte | Enhanced photoacoustic neuroimaging with gold nanorods and PEBBLEs[END_REF], a high frequency piezoelectric ultrasound transducer array was developed for PAT [START_REF] Zemp | Photoacoustic imaging of the microvasculature with a high-frequency ultrasound array transducer[END_REF], in [START_REF] Ermilov | Laser photoacoustic imaging system for detection of breast cancer[END_REF],a hemicylindrical cup shaped acoustic detector array was adopted.

Existing systems in photoacoustic mammography

In photoacoustic mammography, the breast is exposed to short pulses of non-ionizing laser.

Some part of the laser energy can be absorbed and converted into heat, leading to thermal expansion and thus ultrasound emission. The generated ultrasonic waves are then detected by ultrasonic transducers to form images. As blood has higher orders of absorption than surrounding tissues, and tumour development always associates with angiogenesis, it provides an intrinsic marker in photoacoustic imaging for distinguishing tumours from normal tissues.

There exists two photoacoustic tomography systems, one is LOIS-64 [START_REF] Ermilov | Laser photoacoustic imaging system for detection of breast cancer[END_REF], designed, fabricated and tested especially for breast cancer detection; the second one is the Twente Photacoustic Mammoscope (S. [START_REF] Manohar | The Twente Photoacoustic Mammoscope: system overview and performance[END_REF], (Susanne E. [START_REF] Vaartjes | First clinical trials of the Twente photoacoustic mammoscope (PAM)[END_REF], [START_REF] Jose | Passive element enriched photoacoustic computed tomography (PER PACT) for simultaneous imaging of acoustic propagation properties and light absorption[END_REF]developed at the University of Twente. Both of these systems have undergone through several clinical trials for breast cancer diagnosis.

The laser optoacoustic imaging system (LOIS)

It was designed and fabricated by Fairway Medical Technologies, Inc. and The University of Texas Medical Branch, as shown in Fig. 1-7, this system was designed to image a single breast slice with an arc-shaped array of 64 rectangular acoustic transducers. In this system, a Q-switched Alexandrite laser was used as an NIR light source. The laser emitted 75 ns pulse (750 mJ per pulse) at wavelength of 757 nm with the repetition rate of 10 Hz, the maximum fluence of the laser is 10 mJ cm -2 , the output laser beam was coupled into a custom-made optical fiber bundle and expanded with the lens system, producing a Gaussian profile on the surface of the breast. The laser beam incident on the surface of the breast was 70 mm in diameter and had a maximum fluence of 10 mJcm -2 .

The hemicylindrical cup shaped acoustic detector has the radius of 70 mm and width of 90 mm, it consisted of 64 rectangular polyvinylidene fluride (PVDF) transducers (20x3x0.11 mm) and was used for reconstruction of 2-D images of the breast tissue. The electrical signals generated by each transducer were amplified by a two-stage amplifier and sent to an analog-to-digital converter (ADC). The data acquisition system communicated with a PC, the software of the system controlled the data acquisition, system configuration (number of transducers, signal processing), and the realtime display of the signals. The reconstruction algorithm they used was the radical back-projection algorithm. With this system, they have done some preliminary clinical studies on breast cancer patients. This system has demonstrated high-contrast modality for imaging of breast cancer (Fig. 1-8). This system was designed only for visualization of breast cancer and does not achieve any quantitative image and can visualize only a single breast slice, to see a large volume of the breast, several acquisitions have to be taken.

The Twente Photocoustic Mammoscope

Another system was developed by the University of Twente as shown in Fig. 123456789. In this system, the light source is a Q-switched Nd: YAG laser with a pulse duration of 5 ns and a repetition rate of 10 Hz at wavelength of 1064 nm [START_REF] Vaartjes | First clinical trials of the Twente photoacoustic mammoscope (PAM)[END_REF] , the light beam is coupled into a light delivery system. The breast is mildly compressed between a glass window and a flat detector matrix. The PVDF detector sheet is 110 μm thick and approximately 90 mm in diameter, thin gold electrodes have been deposited on the front and rear faces of the PVDF sheet, these 590 electrodes are 2 x 2 mm square and arranged in a roughly circular grid with a centre-centre spacing of 3.175 mm. The system is controlled by a digital input-output card in the PC. 3D images are reconstructed by a delay-and-sum beamforming algorithm.

With this system, they did their first clinical trials on the suspect breast of a 54 year

Conclusion

This A review of possible methods of resolution is finally presented and discussed.

Physical properties of biological tissues

Optical contrasts in biological tissues

Absorption

The absorption properties of a biological medium are described by the absorption coefficient a µ , defined as the probability of photon absorption in the medium per unit path length, it's related to absorption cross section a σ and the number density of the absorbers in the medium a N [START_REF] Wang | Biomedical Opticals: Principles and Imaging[END_REF] :

a a a N σ µ = . Its reciprocal, ( ) 1 - = a a l µ
, is called absorption mean free path length, it corresponds to the mean distance that the photon travels before it is absorbed. In a non scattering homogeneous environment, when a light beam collimated through an environment with a depth of L , the relationship between the transmitted light intensity I and the incident light intensity 0 I is:

) exp( 0 L I I a µ - = (2. 1)
Water, hemoglobin and melanin are the main absorbers in the tissue. Fig. 1-1 displays the absorption coefficient of the main absorbers in tissue.

Scattering

The scattering properties of a biological medium are described by the scattering coefficient s µ , it is defined as the probability of photon scattering in the medium per unit path length, it's related to scattering cross section s σ and the number density of the scatters in the medium s N ,

s s s N σ µ = . Its reciprocal, ( ) 1 - = s s l µ
, is called scattering mean free path length, it corresponds to the mean distance that the photon travels before it scatters. In a non-absorbing homogeneous environment, when a light beam collimated through an environment with a depth of L , the relationship between the transmitted light intensity I and the incident light intensity 0 I is:

) exp( 0 L I I s µ - = (2. 2)
The extinction coefficient t µ , also called total interaction coefficient is defined as:

s a t µ µ µ + = (2. 3) Its reciprocal t t l µ 1 =
is the mean free path between interaction events. s to be scattered in the observation direction s (Da Silva, 2011).

Usually, the phase function ) , ( ' s s p is normalized at 1 by integration over the complete space:

1 ) , ( ' 4 ' = ∫ s s p d π ω (2. 4)
Here ' ω d is the unit solid angle of observation. The mean value of the cosine of the scattering angle θ is defined as the anisotropy factor:

' ' 4 ' ) , ( cos s s s s ⋅ = = ∫ p d g π ω θ (2. 5) 0 = g
stands for a perfectly isotropic scattering and 1 = g stands for a total forward scattering (no scattering). In biological tissues, the scattering is preferentially forward, the anisotropy factor varies from 0.8 to 0.98. For biological tissues, usually, the Henyey-Greenstein phase function is a good approximation [START_REF] Henyey | Diffuse radiation in the galaxy[END_REF]:

2 / 3 2 ' 2 ' ) 2 1 ( 1 4 1 ) , ( g g g p + ⋅ - - = s s s s π (2. 6)
Considering the anisotropy, the reduced scattering coefficient is defined as:

s s g µ µ ) 1 ( ' - = (2. 7)
Consequently, a reduced transport coefficient can be defined as:

' ' s a t µ µ µ + = (2. 8) Its reciprocal ' ' / 1 t t l µ =
is referred to as the transport mean free path. [START_REF] Cheong | A review of the opticalal properties of biological tissues[END_REF] Two other parameters are frequently used: the diffusion coefficient D and the effective attenuation coefficient , defined as follows:

) ( 3 1 ' s a D µ µ + = (2. 9) ) ( 3 ' s a a a eff D µ µ µ µ µ + = =
(2. 10)

Photometric quantities

The spectral radiance v L is the most general physical quantity to describe light propagation in tissue, it's defined as the energy flow per unit normal area per unit solid angle per unit time per unit temporal frequency bandwidth. The radiance L is defined as the spectral radiance integrated over a narrow frequency range [ ]

v v v ∆ + , : v t L t L v ∆ = ) , , ( ) , , ( s r s r [ 1 2 - -sr Wm
], r is the position, s is the direction unit vector, t is the time.

The fluence rate (intensity) φ is defined as the energy flow per unit area per unit time regardless of the flow direction; it is expressed as the radiance integrated over the π 4 solid angle:

∫ = π ω φ 4 ) , , ( ) , ( d t L t s r r ( 2 - Wm ) (2. 11)
Fluence F is defined as the time-integrated fluence rate:

∫ +∞ ∞ - = dt t F ) , ( ) ( r r φ ( 2 - Jm ) (2. 12)
Current density is defined as the net energy flow per unit area per unit time, it's defined as:

∫ = π ω 4 ) , , ( ) , ( d t L t J s r s r ( 2 - Wm ) (2. 13) eff µ

The Radiative Transport Equation

In the radiative transport equation, the optical properties (such as the absorption coefficient, the scattering coefficient, the anisotropy factor) are assumed to be time invariant but space variant, and only elastic scattering is considered in this model. The propagating energy per unit volume per unit solid angle can be obtained by taking the ratio between the radiance and the speed of light. As a result, we can write the energy variation in the considered volume element, within the infinitesimal solid angle per unit time, as:

ω ω dVd t c t L dsdAd t c t L dP ∂ ∂ = ∂ ∂ = / ) , , ( / ) , , ( s r s r (2. 14)

Here

is the volume of the infinitesimal cylinder.

The energy variation dP is actually given by the algebraic sum of terms that account for the increase and decrease of energy inside the cylinder for different terms: divergence, extinction, scattering, and source.

src sca ext div dP dP dP dP dP + + - - = (2. 15)

Contribution 1: divergence

If the photon beam is not collimated, divergence is nonzero. Energy diverging in the considered volume element, within the solid angle per unit time is given as:

dAds dV = dV d s t L dP div ω ∂ ∂ = ) , , ( s r (2. 16)
The contribution is positive for actual divergence and negative for actual convergence.

In divergence form, Eq. (2.16) becomes

dV d t L dV d t L dP div ω ω ] ) , , ( .[ ) , , ( s s r s r s ∇ = ∇ ⋅ = (2. 17
) , ( ) , ( ' ' s s s s p p = (2. 20) Because s s s N σ µ = ω ω µ π d d p t L dV dP s sca ] ) , ( ) , , ( )[ ( ' 4 ' ∫ = s s s r (2. 21) Contribution 4: Source ] ) , , ( )[ ( ω µ dAd t L ds dP t ext s r =
Energy produced by a source in the volume element within the solid angle element per unit time is given by

ω dVd t Q dP src ) , , ( s r = (2. 22)

Final form of radiative transfer equation

Substituting Eq. (2.17), Eq. (2.18), Eq. (2.21), and Eq. (2.22) into Eq. (2.15), we get

) , , ( ) , ( ) , , ( ) , , ( ) , , ( ) , , ( 1 '
' 4 t Q d p t L t L t L s t t L c s t s r s s s r s r s r s r + + - ∇ ⋅ - = ∂ ∂ ∫ ω µ µ π (2. 23)
Eq. (2.23) is the radiative transfer equation.

For the time-independent responses,

0 / ) , , ( = ∂ ∂ t c t L s r
(2. 24)

The Diffusion Equation

The radiative transfer equation is difficult to solve because it has six independent variables:

t z y x , , , , , φ θ
. However, it can be simplified provided some assumptions are fulfilled.

The radiance can always be expanded to a series of spherical harmonics m n Y , of the first order:

) ( ) , ( ) , , ( , 0 , s r s r m n n n n m m n Y t L t L ∑ ∑ +∞ = - = ≈ (2. 25)
Within the so-called 1 P approximation, the series is truncated at order 1 and reduces to:

s r r s r ⋅ + ≈ t) , ( 4 3 t) , ( 4
1 t) , , ( J L π φ π (2. 26)
Noting the fluence rate can be expressed as

) ( t)Y , ( L 4 t) , ( 0,0 0,0 s r r π φ =
and the net energy flow as

) ( ) , ( 3 4 t) , ( , 1 1 1 , 1 s r s r J m m m Y t L ∑ - = = ⋅ π .
Within the Diffusion Approximation, the radiance is supposed to be nearly isotropic after sufficient scattering. This can be achieved when the biological tissue is a highly scattering medium (

s a µ µ ′ <<
) and the probed medium is large with its interfaces are far from any sources. If the medium is the breast probed with a light source whose wavelength is in the red or near infrared, these assumptions are met. The source can be assumed to be isotropic. Substituting Eq. (2.26) into Eq. (2.23) and integrating over the π 4 solid angle, we obtain the following scalar differential equation:

) , ( ) , ( ) , ( ) , ( 1 t Q t J t t t c a r r r r = ⋅ ∇ + + ∂ ∂ φ µ φ (2. 27)
Substituting Eq. (2.26) into Eq. (2.23), multiplying both sides by s , and integrating over the π 4 solid angle, we get the vector differential equation:

0 t) , ( 3 1 t) , ( ) ( t t) , ( 1 ' s a = ∇ + + + ∂ ∂ r r r φ µ µ J J c (2. 28) Assume that the fractional change in ) , ( t J r is small, ) , ( ) ( ) , ( ' t J t c t J s a r r µ µ + << ∂ ∂ , one obtains the Fick's law: ) , ( ) , ( t D t J r r φ ∇ - =
(2. 29) Substituting Eq. ( 2.29) into Eq. ( 2.27), we get the Diffusion Equation )

, ( )] , ( [ ) , ( ) , ( 1 t Q t D t t t c a r r r r = ∇ ⋅ ∇ - + ∂ ∂ φ φ µ φ (2. 30)

Photoacoustic phenomenon

When the tissue absorbs energy, the absorbed energy converts into heat and the temperature of the tissue increases, this causes thermal expansion, then the acoustic pressure is generated. There are two important timescales in light heating:

The thermal relaxation time

th c th d α τ 2 = , th α is the thermal diffusivity ( 1 2 - s m
) and c d is the characteristic dimension of the heated region;

The stress relaxation time can be expressed as

s c s v d = τ , s v is the sound speed.
If the laser pulse is much shorter than th τ , the heating process is said to be in thermal confinement, which means that the heat conduction is negligible during the light heating. If the laser pulse is much shorter than s τ , the heating process is said to be in stress confinement, which means that the stress propagation is negligible during the light heating.

In this discussion, we assume that the laser pulse meets both thermal confinement and stress confinement; the tissue is assumed to be acoustically homogeneous and isotropic, which means that its mechanical properties, like density and speed of sound are assumed to be constant in the medium.

To get the expression of the initial pressure, the fractional volume expansion is introduced, it's expressed as:

T p V dV β κ + - = (2. 31)
Here, κ is the isothermal compressibility, β is the thermal coefficient of volume expansion, p is the change of photoacoustic pressure and T is the change of temperature. The isothermal compressibility κ can be expressed as

v s p C v C 2 ρ κ = (2. 32)
Here, ρ is the mass density, p C and v C are specific heat capacities at constant pressure and volume, respectively. Under the thermal and stress confinements, the fractional volume expansion are negligible and the initial pressure can be expressed as:

κ βT p = 0 (2. 33)
From thermodynamics, the relation between temperature variation and internal energy

variation E ∆ is T mC E v = ∆ , m is mass. E ∆ depends also on the absorbed intensity H , HV E = ∆ , leading to v v C H mC HV T ρ = =
. Substituting this expression to Eq. (2.33), we get:

H H C v H C p p s v Γ = = = β κρ β 2 0 (2. 34)
Here, Γ is the Grüneisen parameter (dimensionless).

When the tissue is heated by a laser pulse, an initial pressure distribution will be generated due to the thermal expansion, this initial pressure starts immediately propagating into the medium.

In order to establish the photoacoustic equation, we begin with the thermal expansion equation:

) , ( ) , ( ) , ( t T t p t r r r β κ ε + - = ⋅ ∇ (2. 35)
Here, ε is the medium displacement. We also need the linear inviscid force equation, expressed as:

) , ( ) , ( 2 2 t p t t r r ε -∇ = ∂ ∂ ρ (2. 36)
The left-hand of the equation is the mass density times the displacement, while the righthand of the equation represents the applied force per unit volume. The divergence of Eq. (2. 35) is:

) , ( ] ) , ( [ 2 2 2 t p t t r r ε -∇ = ⋅ ∇ ∂ ∂ ρ (2. 37)
We substitute Eq. (2.35) into Eq. (2.37), get

( ) ( ) 2 2 2 2 2 , , t t T t p t ∂ ∂ - =         ∂ ∂ - ∇ r r βρ κρ (2. 38)
As the speed of sound can be rewritten as

κρ / 1 = s v , Eq. (2.38) becomes ( ) ( ) 2 2 2 2 2 2 2 , , 1 t t T v t p t v s s ∂ ∂ - =         ∂ ∂ - ∇ r r κ β (2. 39)
This is the general form of photoacoustic equation. The left-hand side of the equation describes wave propagation, while the right-hand side represents the source term. As it is assumed that we are under thermal confinement, we can use the thermal equation:

) , ( ) , ( t H t t T C v r r = ∂ ∂ ρ (2. 40)
Here, ) , ( t H r is the heating function, defined as the thermal energy converted from light per unit volume and per unit time. For optical illumination, the heating function is the product of absorption coefficient and light fluence:

) , ( ) ( ) , ( t t H a r r r Φ = µ
(2. 41)

If we substitute Eq. (2.40) into Eq. (2.39), we get

( ) ( ) 2 2 2 2 2 2 , , 1 t t H C t p t v p s ∂ ∂ - =         ∂ ∂ - ∇ r r β (2. 42)
Here, source depends on the first derivative of ( ) t H , r .This means that the photoacoustic wave is generated only if the heating function varies in time. So to generate photoacoustic pressure, only pulsed or modulated light is used as the light source.

Resolution methods

Resolution of the Diffusion Equation

Analytical method

The analytical solution to the Diffusion Equation is related to the geometry of the medium. [START_REF] Arridge | The theoretical basis for the determination of opticalal pathlengths in tissue: temporal and frequency analysis[END_REF] established the Green function solutions to the Diffusion Equation for various simple geometry media, in the time and frequency domains: infinite medium, semi-infinite halfspace, infinite slab, 2D circle, finite cylinder, infinite cylinder, and sphere. As an example, the expression of the time domain Green's function in an infinite medium is:

[ ] { } { } )] ( 4 / ) ( [ exp ) ( ) 4 ( 2 / ) , , ( G ' 2 ' 2 / 5 ' 2 / 3 ' inf t t cD d t t c t t cD d t t d a - + - - - = µ π (2. 43)
D is the distance between source and detector. Analytical solutions to the RTE can be found but are restricted by the geometry of the medium, and only solution to certain simple geometries can be obtained.

Numerical methods

Monte Carlo method

The Monte Carlo method is a widely used numerical method to solve the radiative transport problem, it's a stochastic method, it records photons histories as they are scattered and absorbed.

The flow chart of Monte Carlo is shown in Fig. 23. The Monte Carlo method begins by launching a photon into the tissue, the photon's direction is chosen randomly from all possible directions. Then the photon propagated with small, fixed stepsize, the stepsize must be small relative to the mean free pathlength of a photon in the tissue:

s a t s µ µ µ + = << ∆ 1 1 (2. 44)
The possibility of internal reflection occurs when the photon is propagated across a boundary into a region, the possibility is determined by the Fresnel reflection coefficient:

] ) ( tan ) ( tan ) ( sin ) ( sin [ 2 1 ) R( 2 2 2 2 i t i t i θ θ θ + - + + - = (2. 45)
Here, i θ is the angle of incidence on the boundary and t θ is the angle of transmission. After each propagation step, the photon is split into two parts: one part is absorbed and the other is scattered, the part that is absorbed:

s a s s a a absorbed f µ µ µ µ µ µ + - = + = 1 (2. 46)
The photon will never reach zero, and propagating a photon with a minuscule weight yields little information. A technique called roulette is used to terminate a photon once its weight drops below a specific minimum. The Monte Carlo method that describes each photon propagation is accurate but computationally inefficient as soon as the scattering level is high and/or the medium is large.

Finite Element Method (FEM)

The FEM divides the geometry into small subdomains, called meshes, then find approximate solutions to the equations to be solved.

Hereafter the method is illustrated through the resolution of the Diffusion Equation. We seek a continuous and piecewise linear approximation h φ ofφ . The geometry Ω into subdomainsτ , Nodes ) ,..., 1 ( P j N j = are attached to the element vertices, h φ at each point r within an element i τ is

given by a linear interpolation of nodal values j

φ , ) ( ) ( ) , ( | r r j N j j h t t i j ψ φ φ τ ∑ ∈ =
where j ψ are linear nodal shape functions with support over all elements which have the node j N as a vertex. The weak formulation of the diffusion equation using the Galerkin approach is given by [START_REF] Arridge | A finite element approach for modeling photon transport in tissue[END_REF])

Ω = Ω + ∇ ⋅ ∇ - ∂ ∂ ∫ ∫ Ω Ω d Q d t r D t c j h a ) ( ) ( ) , ( )] ( ) ( 1 )[ ( j r r r r r ψ φ µ ψ (2. 47)
Integration by parts and substitution of Eq. (2.78) leads to

) ( ) , ( ) ( 1 ) ( ) ( )] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( 1 ) ( ) ( [ Ω ∂ Γ - Ω = Ω + ∇ ⋅ ∇ + ∂ ∂ ∫ ∫ ∫ Ω ∂ Ω Ω d t c d Q d t t D t t c j j i i j a i i i i i j ξ ξ ψ ψ φ ψ ψ µ φ ψ ψ φ ψ ψ r r r r r r r r r r (2. 48)
This can be written in matrix form:

β φ µ φ - = + + ∂ ∂ Q C K B )] ( ) ( [ a D t (2. 49)
where

Ω ∇ ⋅ ∇ = ∫ Ω d D K i j ij ) ( ) ( ) ( r r r ψ ψ Ω = ∫ Ω d t Q t Q j j ) , ( ) ( ) ( r r ψ Ω = ∫ Ω d c C i j a ij ) ( ) ( ) ( r r r ψ ψ µ (2. 50) ) ( ) , ( ) ( Ω ∂ Γ = ∫ Ω ∂ d t j j r r ψ β Ω = ∫ Ω d M i j ij ) ( ) ( r r ψ ψ
Besides FEM, there are Finite Difference Method (FDM), approximating the solutions to differential equations using finite difference equations to approximate derivatives, and Finite Volume Method (FVM), "Finite volume" refers to the small volume surrounding each node point on a mesh.

In the FVM, volume integrals in a partial differential equation that contain a divergence term are converted to surface integrals, using the divergence theorem. These terms are then evaluated as fluxes at the surfaces of each finite volume. FEM, FDM, and FVM are easy to implement, but the accuracy strongly depends on the mesh. The FEM is chosen in this manuscript to solve the DE.

Acoustic resolution

Analytical

When the medium is acoustically homogeneous, the general way to solve the photoacoustic equation is by using the Green function. The Green function is defined as the response to a spatial and temporal impulse source:

) ( ) ( ) , ; , ( ) 1 ( ' ' ' ' 2 2 2 2 t t t t G t v s - - - = ∂ ∂ - ∇ δ δ r r r r (2. 51)
Here, ' r is the position of the source and ' t is the time of the light pulse. In infinite medium, the Green's function is given by

' ' ' ' ' 4 ) ( ) , ; , ( r r r r r r - - - - = π δ s v t t t t G (2. 52)
Which represents an impulse diverging spherical wave, and

) , ; , ( ) , ; , ( ' ' ' ' t t G t t G - - = r r r r
.

Numerical

The acoustic wave equation is also a partial differential equation, it can be solved by finite element method or finite difference method. For both finite element method and finite difference method, they become slow when solving large scale, high-frequency acoustic waves, because it requires many mesh points for one wavelength, and small time steps to minimize the dispersion. For the gradient of the field, the FD should fit a high order polynomial to get an accurate solution.

The k-space method (Cox et al., 2007a) fits a Fourier series to all the data and the gradient can be calculated simply as

{ } { } (x) 1 - p F ik F
x , F stands for the Fourier transform, x k is the spatial wave vector in k-space, and the basis functions in Fourier transform are sinusoidal, so only two points per wavelength are needed to describe a wave.

photoacoustic resolution

The photoacoustic is the combination of optics and acoustic. The distribution of the photoacoustic initial pressure is related to the light fluence in the light propagation.

Applying the Green function approach, the photoacoustic equation can be solved as

) , ( 1 4 )
, ( 1 4 ) , ( ' ' ' ' ' ' ' ' s p s p v t H dr t C t v t H d C t p r r r r r r r r r r r r - - - ∂ ∂ = ∂ - - ∂ - = ∫ ∫ π β π β (2. 54)
For a short laser pulse that meets the acoustic stress confinement, the heating function can be treated as a delta function:

) ( ) ( ) , ( ' ' ' ' t A t H δ r r ≈
, so Eq. (2.54) can be expressed as [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF])

)] ( ) ( 1 [ 4 )] ( ) ( 1 [ 4 ) , ( ' ' 0 ' 2 ' ' ' s s s s s p v t p d t v t v v t A d t v t C t p r r r r r r r r r - - ∂ ∂ = - - ∂ ∂ = ∫ ∫ δ π β δ π β (2. 55)
From Eq. (2.55), we can see that, the photoacoustic pressure at position r and time t comes from sources over a spherical shell centered at r with a radius of t v s . . The sound speed is 1 1500 - ms . The waveform is called 'N-shape' wave.

Fig. 2-5 influence of sound speed to photoacoustic pressure from a uniform spherical source at a distance of 20 mm with a radius of 2 mm, with the sound speed of 1500 ms -1 and 1200 ms -1 , respectively.

Fig. 2-5 tested the influence of sound speed to the photoacoustic pressure, the sound speed just influences the arrival time of the pressure, the higher the speed, the earlier the arrival time; the sound speed does not influence the amplitude of the pressure. Fig. 2-7 tested the influence of the detection distance to the photoacoustic pressure; we can see that, the amplitude of the pressure is inversely proportional to the distance from the source.

In a word, the amplitude of the pressure is approximately proportional to the radius of the source and inversely proportional to the distance from the source. The width of the waveform is proportional to the radius of the source. The zero point of the waveform indicates the distance from the source. Many reconstruction methods for solving the inverse problem exist, hereafter are listed the most used.

Back projection method

Among various techniques for photoacoustic image reconstruction, back projection method is the most commonly used one. In this method, the projections of the target are taken from every angle around the target, to reconstruct the image, the projections are equally re-distributed back to its original path, and the amount of the projection is the highest of all the projections. If we backproject from all angles over °360 , we will get an image similar to the original one. After the backprojection, we will get a blurred version of the original image. To sharpen the image, we apply a negative wing before backprojection to filter the projections, after this process, we will get a clear image. This image reconstruction algorithm is called the Filtered Back projection Algorithm. For parallel beam geometry, the total attenuation can be represented by a line integral ( ) ) , (

- + = = ∫ ∫ ∫ ∞ ∞ - ∞ ∞ - θ θ δ θ θ (2. 56) D1 D2 s 1 t v s * 2 t v s *
+ - ∞ ∞ - ∞ ∞ - - ∞ ∞ - ∫ ∫ ∫ = = (2. 57)
Several Filters can be used in this step such as: Ramp (standard), Cosine, Hamming, and SheppLogan. [START_REF] Faria | Performance study of filtered back-projection algorithms implemented on GPUs[END_REF]. If the ramp filter defined as ( ) ω ω ρ = is used, the projection becomes,

( ) ( ) ( ) ω ρ ω ω θ θ S M =
. The inverse Fourier transform of ( )

ω θ M : ω ω πω θ θ d e M t Q t i 2 ) ( ) ( ∫ ∞ ∞ - = (2. 58)
The backprojection is represented as:

( ) θ θ θ π θ d y x Q y x f ) sin( ) cos( ) , ( 0 + = ∫ (2. 59)
A temporal back projection reconstruction algorithm in various geometries is discussed in [START_REF] Xu | Time-domain reconstruction algorithms and numerical simulations for thermoacoustic tomography in various geometries[END_REF], [START_REF] Burgholzer | Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface[END_REF] use specially shaped detectors that are larger than the imaged object to acquire acoustic signals. This approach is based on the fact that such a detector receives a signal that is not an approximate projection over a spherical surface but rather an exact projection over an area that is determined by the shape of the detector itself.

Delay-and-sum beam-forming algorithm

The idea of delay-and-sum beam-forming algorithm is easy to understand: If a propagating signal is present in an array's aperture, the sensor outputs, delayed by appropriate amounts and added together, reinforce the signal with respect to noise or waves propagating in different directions (Hoelen and de Mul, 2000). The delay-and-sum beam-forming algorithm points the array's sensitivity in the direction of view can be steered toward the acoustic source by adjustment of the sensor arrays. The signal from each voxel of the imaged volume is expressed as (A. [START_REF] Kharine | Towards photoacoustic mammography[END_REF]:

( ) M t S w t S M i f i i f i f ∑ + = τ ) ( (2. 60) Where ) (t S f is the signal from th f focus; ( ) t S i is the signal from th i detector; f i τ is the time delay of
the fth focus to the th i detector; f i w is the amplitude weightfactor of th f focus to th i detector; M is the number of detectors. In 3D, where the sources are in the near field of the detector array, the transients cannot be treated as plane waves as often is so convenient in beam-forming algorithm applications [START_REF] Dudgeon | Fundamentals of digital array processing[END_REF], we should modify the delay-and-sum beam-forming algorithm:

( )

∑ ∑ + = M i f i M i f i i f i f w t S w t S τ ) ( (2. 61)
Where the summed signal is normalized to make the output independent of the actual set of transducers. Klemm at al. [START_REF] Klemm | Improved Delay-and-Sum Beamforming Algorithm for Breast Cancer Detection[END_REF] proposed a modified delay-and-sum beamforming algorithm for breast cancer detection, the algorithm uses an additional weight factor calculated at each focal point to improve image quality.

Iterative reconstruction method

Unlike the analytical reconstruction algorithms where each projection sample is weighted, filtered, and backprojected to build an image, iterative reconstruction arrives at the final solution in an iterative manner: the initial reconstructed images are refined and modified iteratively until certain criteria are met. The criteria are written in the form of a cost function to be minimized, which measures the fit of the image to the data according to a model of the imaging system. [START_REF] Hsieh | Recent Advances in CT Image Reconstruction[END_REF]. However, the large size of the optimization problem as well as the complexity of the model makes it difficult to solve. Each iteration involves the forward and back-projection of an image. The forward-projection simulates the light interactions with the object and produces a set of synthesized projections. The synthesized projections are compared against the real projection measurements, and the differences between the two are attributed to the error. The error is then used to update the intermediate image volume to reduce the discrepancy between the image and the acquired data.

The new intermediate image is then forward-projected again in the next iteration, provided the choice of an adequate cost function and a globally convergent optimization algorithm.

Time reversal method

The principle of time reversal reconstruction method lies in that, for a sufficiently large time T , the photoacoustic pressure ( ) t r p , can be set be zero for any T t > . So one can impose zero condition on ( ) t r p , for T t = and boundary conditions to solve the wave equation back in time. Xu (Y. Xu and L. V. [START_REF] Xu | time reversal and its application to tomography with diffracting sources[END_REF] proposed an exact reconstruction for three-dimensional broadband diffraction tomography. As time reversal is exact only in constant sound speed in odd dimension, (Y. [START_REF] Hristova | reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media[END_REF] proved that it works well even under the most unfavorable circumstances (variable, and even trapping sound speed). (Treeby et al., 2010) tested time reversal reconstruction method with success in absorbing acoustic media.

Although back projection is the mostly commonly used method in photoacoustic image reconstruction, it works only when the sound speed is constant; the iterative method always needs enormous computation; while the time reversal method works for arbitrary geometry of the closed surface, the sound speed can be variable, it's seen as the least strict reconstruction method, so in our work, we used time reversal method to do the reconstruction of the initial pressure distribution.

Quantitative PAT (QPAT)

The first step of the inverse problem in PAT is to reconstruct the initial pressure map, this process is well known and much work has been done in this area. The initial pressure map,

)) ( ), ( , ( ) ( ) ( ) ( 0 r D r r r r r p a a µ φ µ Γ =
,is the product of the Grüneisen coefficient, the optical absorption coefficient and the optical fluence which depends itself on the local optical properties of the medium, in terms of absorption and diffusion. Hence, in order to go further in the quantification through an accurate estimation of the absorption coefficient, uncoupling the effects of absorption and diffusion is necessary. This is the main goal of quantitative photoacoustic tomography (QPAT: complete uncouple the maps of the absorption, scattering and Grüneisen parameters).

Karabutov's work [START_REF] Karabutov | Time-resolved laser photoacoustic tomography of inhomogeneous media[END_REF] might be the first that refers to the QPAT, in their work, the distribution of light absorption was obtained from the temporal course of acoustic pressure; In Oraevsky's work [START_REF] Oraevsky | Measurement of tissue opticalal properties by time-resolved detection of laser-induced transient stress[END_REF], the absorption and effective scattering coefficients of beef liver, dog prostate, and human aortic atheroma at three wavelengths were deduced from laser-induced stress profiles with additional measurements of total diffuse reflectance.

the distribution of absorbed energy and the absorption coefficient in the medium were measured from the photoacoustic waves with a flat top laser beam profile in [START_REF] Paltauf | Pulsed photoacoustic characterization of layered media[END_REF]; [START_REF] Jaeger | Diffraction-free acoustic detection for photoacoustic depth profiling of tissue using an opticalally transparent polyvinylidene fluoride pressure transducer operated in backward and forward mode[END_REF] presented an photoacoustic detection method suitable for depth profiling of optical absorption of layered or continuously varying tissue structures. However, the above mentioned works are limited to 1D case. In 2D, Cox et al. [START_REF] Cox | Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method[END_REF]) used a fix-point iterative inversion scheme for recovery of absorption coefficient when the scattering distribution is known; [START_REF] Yuan | Quantitative photoacoustic tomography: Recovery of opticalal absorption coefficient maps of heterogeneous media[END_REF] recovered the absorption coefficient from a finite element reconstruction method; [START_REF] Ripoll | Quantitative point source photoacoustic inversion formulas for scattering and absorbing media[END_REF] recovered the absorption map with a point source using an iterative diffusion regime inverse method, in their work, then also accounted for instrumental factors such as the source strength, the shape of the optical pulse, and the impulse response and finite size of the transducers; [START_REF] Banerjee | Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of opticalal absorption coefficient from the reconstructed absorbed energy map[END_REF] proposed a noniterative method to recover the absorption coefficient directly from boundary pressure measurements; [START_REF] Habib Ammari | Mathematical modelling in photo-acoustic imaging of small absorbers[END_REF]reconstructed small absorbing regions and absorbing energy density inside a bounded domain from boundary data; (Bal et al., 2010) developed and used an inverse transport theory with internal measurements to extract information on the optical coefficients from knowledge of the deposited thermal energy map.

By introducing prior knowledge on the spatial frequency of the two quantities, [START_REF] Rosenthal | Quantitative Photoacoustic Signal Extraction Using Sparse Signal Representation[END_REF]extracted both the absorption coefficient (high frequency) and the fluence (low frequency) from the photoacoustic image. As this algorithm is not based on the explicit solution the theoretical light transport equation, it does not require explicit knowledge of the illumination geometry.

The above noted literature has considered only a single optical illumination, (Zemp, 2010a) recovered the absorption perturbation when the optical properties of the turbid background are known using a noniterative reconstruction method with multiple optical illuminations; [START_REF] Ranasinghesagara | Photoacoustic technique for assessing opticalal scattering properties of turbid media[END_REF] describes a simple method for estimating optical scattering properties of turbid media using multiple illumination.

With multiple illumination, both the absorption and diffusion coefficients can be reconstructed simultaneously. [START_REF] Tarvainen | Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography[END_REF] reconstructed both the absorption and diffusion coefficients, by considering both the radiative transport and diffusion equations as light transport models; [START_REF] Bal | On multi-spectral quantitative photoacoustic tomography in diffusive regime[END_REF] recovered the absorption, diffusion, and Grüneisen coefficients mathematically when the propagation of radiation is modeled by a second-order elliptic equation; [START_REF] Shao | Estimating opticalal absorption, scattering, and Grüneisen distributions with multiple-illumination photoacoustic tomography[END_REF] estimated the absorption, scattering and Grüneisen distributions from reconstructed initial pressure map with multiple optical sources. In their updated work (Shao et al., 2012a), they recovered the absorption and diffusion coefficients directly from measured acoustic pressure.

The absorption and diffusion coefficients can be reconstructed simultaneously by introducing prior knowledge (spectra of known species, diffusion spectral profile) with multiple wavelength illuminations; [START_REF] Bal | On multi-spectral quantitative photoacoustic tomography in diffusive regime[END_REF] have proven that, when multiple wavelength data are available, the absorption, diffusion and Grüneisen coefficients can be reconstructed simultaneously under additional minor prior assumptions.

We classified the approaches in QPAT into four categories:

-the absorption coefficient supposing the diffusion coefficient is known.

-the absorption coefficient and the fluence.

-the absorption and diffusion coefficients simultaneously by introducing active probing with multiple optical illuminations.

-the absorption and diffusion coefficients simultaneously by introducing prior knowledge (spectra of known species, diffusion spectral profile) with multiple wavelength illuminations.

Conclusion

In this chapter, we firstly introduced and defined the physical properties of the biological tissue, and establish the equation governing the photoacoustic phenomenon. The forward multiwave photoacoustic model is a combination of optical and acoustic propagation equations. For the optical propagation, we chose the finite element method; for the acoustic propagation, instead of the FEM that can be very slow when solving large time steps, the k-space method is preferred. The In the forward model the optical and acoustic propagation are weakly coupled through the source term of the acoustic wave equation. For solving the optical propagation, the Finite Element Method (FEM) is chosen, this leads to the calculation of the initial pressure distribution map; and the k-space resolution method is used for the acoustic propagation.

In the inverse problem, reconstruction techniques conventionally used in ultrasound imaging, like the time-reversal method, can be used to reconstruct the initial pressure map from real or synthetic data. However to have access to the reconstruction of optical and acoustic parameters and to achieve a real parameter quantification (QPAT stands for Quantitative PAT), more sophisticated approaches have to be adopted, including modeling of the optical path in the tissues. We discuss two approaches introduced very lately by [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF][START_REF] Zemp | Quantitative photoacoustic tomography with multiple opticalal sources[END_REF]: one consists in recovering the uncoupled absorption and diffusion parameters from the reconstructed initial pressure map; the other allows us to obtain directly the optical parameters maps directly from the measured photoacoustic pressure.

Multiphysics models implementation

Introduction to the Finite Element Method

We solve the optical equation (

) ( )) ( ) ( .( ) ( ) ( r r r r r Q D a = ∇ ∇ - φ φ µ
) with finite element method (FEM).The starting point of the FEM is the mesh. In 1D, the mesh elements are discrete points, called mesh vertices; in 2D, the geometry is discretized into triangular or quadrilateral mesh elements, if the boundary is curved, these elements represent only an approximation of the original geometry; in 3D, the geometry is discretized into tetrahedral, hexahedral, prism, or pyramid mesh elements whose faces, edges, and corners are called mesh faces, mesh edges, and mesh vertices, respectively, the boundaries in the geometry are discretized into triangular or quadrilateral boundary elements. Once the mesh is created, an approximation to the variables is sought. Usually, the variables are approximated with a certain number of parameters, called degrees of freedom (DOF).

Inserting this approximation into the weak form of the equation generates a system of equations for the degrees of freedom. Once the DOF are fixed, the corresponding basis functions are chosen, which are certain piecewise linear functions or piecewise polynomial functions. Now consider an example of 1D problem P1 to better illustrate the process of the finite element method.

P1:

   = = = , 0 ) 1 ( ) 0 ( ) 1 , 0 ( ) ( ) ( ' ' u u in x f x u (3. 1)
Here f is given, u is unknown.Firstly, we convert P1 into its weak form:

dx x v x u dx x v x f ) ( ) ( ) ( ) ( 1 0 1 0 ' ' ∫ ∫ = (3. 2)
Here, ) (x v satisfies the boundary conditions, that is,

0 = v at 0 = x and 1 = x . Hence: ∫ ∫ ∫ ∫ - = - = - = = 1 0 ' ' 1 0 ' ' 1 0 ' 1 0 1 0 ' ' ) , ( ) ( ) ( ) ( ) ( | ) ( ) ( ) ( ) ( ) ( ) ( v u dx x v x u dx x v x u x v x u dx x v x u dx x v x f φ (3. 3)
Now is the discretization of the problem. We should replace the infinite dimensional linear problem:

Find 1 0 H u ∈ , 1 0 H is a suitable space, such that ∫ = - ∈ ∀ fv v u H v ) , ( , 1 0 φ With a finite dimensional version: Find V u ∈ such that ∫ = - ∈ ∀ fv v u V v ) , ( , φ
Where V is a finite dimensional subspace of 1 0 H .

We take the interval

) 1 , 0 (
, choose n values of x with 1 ... 0

1 1 0 = < < < < = + n n x x x x
and we define V by:

} 0 v(1) v(0) and 0,...n, k for linear | , : R 1] [0, : {v ] x , [x 1 k k = = = → = + v continuous is v V
Now, we should choose a basis function, for each control point k x , we should choose the piecewise linear function k v in V whose value is 1 at k

x and zero at every j x , k j ≠ like: , and if we let

         ∈ - - ∈ - - = + + + - - - otherwise x x x if x x x x x x x if x x x x x v k k k k k k k k k k k 0 ] , [ ] , [ ) ( 1 1 1 1 1 1 (3. 4) If we write ) ( ) ( 1 x v u x u k n k k ∑ = = and ) ( ) ( 1 x v f x f k n k k ∑ = = , taking ) ( ) ( x v x v j = , for n j , ... , 1 = , then Eq. (3- 3), becomes dx v v f v v u j n k k k j k n k k ∑ ∫ ∑ = = = - 1
) ( ij L L = and
) ( ij M M =
be matrices whose entries are ) , (

j i ij v v L φ = and dx v v M j i ij ∫ =
Then we may express Eq. (3.5) as

f u M L = - (3. 6) It's not necessary to assume ) ( ) ( 1 x v f x f k n k k ∑ = = . For a general function ) ( x f , Eq. (1-75) with ) ( ) ( x v x v j =
, for n j , ... In a word, we choose a mesh for the domain to be solved, then is the discretization of the weak form of the equation to be solved, the discretization is performed by the chosen basis function, then we can form a matrix form of the problem, by solving the matrix, we can get the variable we want.

In FEM, instead of solving the equation, the weak form of the equation is usually solved.

Weak form means, instead of solving a differential equation of the underlying problem, an integral function is solved. The integral function implicitly contains the differential equations, however it's much easier solving an integral function than solving a differential function. Also, the differential equation of system poses conditions that must be satisfied by the solution (strong form), whereas, the integral equation states that those conditions need to be satisfied in an average sense (weak form).

Resolution of the optical forward model

The photoacoustic phenomenon can be divided into optical propagation and acoustic propagation. As the light speed is much higher than the sound speed, the optical propagation can be seen as instantaneously compared with acoustic propagation, so the optical propagation can be solved in continuous wave (CW) domain.

) ( )) ( ) ( .( ) ( ) ( r r r r r Q D a = ∇ ∇ - φ φ µ (3. 8)
We solved Eq. (3.8) based on FEM, by using the COMSOL multiphysics software (2008),according to the following steps:

Step (1) Choose a model, in our case, it's a coefficient form of partial differential equation (PDE), the basis function is Lagrange basis function;

Step (2) Draw the geometry, Fig. 3-1 gives an example of geometry in 2D; Fig. 3-1 the breast has a radius of 25 mm, we used water to simulate the surrounding environment, this environment is a square with a side length of 60 mm, the tumour located (10 mm, 10mm) away from the centre of the breast, the optical source located 2 mm from the boundary of the breast.

Step (3) Set the piecewise constant parameters of the equation in different sub domains, Table 3.1 listed the optical parameters to do the simulation; The absorption and reduced scattering coefficients of the surrounding medium corresponds to the average optical properties of the breast at certain wavelength. Water is supposed to be non absorbing and non scattering. The tumour is supposed to absorb twice more than the surrounding medium (healthy tissue), the diffusion coefficient

) 3( 1 D ' s a µ µ + =
, the optical source is modeled as a point source.

Note: Implementation of point source

In our work, the optical source is treated as a point source, and can be modeled as a δ function, so we can rewrite the optical equation as: We followed the recommendation of Comsol for implementation of a point source (http://www.comsol.com/model/implementing-a-point-source-using-poisson-s-equation-73). The easiest way to describe a point source is using the weak form of the FEM. To get weak form of Eq. (3-9), we should multiple these equation with a test function test u and integrate over the computation domain.

   Ω ∂ = ∇ ⋅ + Ω - = ∇ ∇ - 0 ) ( 2 ) ( ) ( )) ( ) ( .( )
     = ⋅ ∀ ⋅ - = ⋅ ∇ ∇ - ∫ ∫ ∫ Ω ∂ Ω Ω 0 ) ( ) ( ) ( ) ( ) ( ))) ( ) ( .( ) ( ) ( ( 0 r r r r r r r r r r test test test s test a u u u u S u D δ φ φ µ (3. 10) As ∫ Ω = ⋅ - ) ( ) ( ) ( s test test s u u r r r r δ
, the following set of equations is now considered Step (5) Choose a mesh, Fig. 3-2 showed the mesh we used in our simulation;

       Ω ∂ = ∇ ⋅ + = ⋅ ∀ = ⋅ ∇ ∇ - ∫ ∫ Ω ∂ Ω 0 ) ( 2 ) ( 0 ) ( ) ( ) ( ) ( ))) ( ) ( .( ) ( ) ( (

The K-space method

The acoustic wave equation is also a partial differential equation, it can be solved by the finite element method or the finite difference method. The flexibility in the choice of basis functions in FEM allows domains and heterogeneities of any shape to be well approximated (Cox et al., 2007a).

In finite difference method (FD), the derivative of the partial differential equations are approximated by differences, finite difference method is less flexible than finite element method, because it usually requires a regular mesh. Both finite element method and finite difference method become cumbersome and slow when solving large scale, high-frequency acoustic waves.

The pseudospectral (PS) method can solve the problem of being cumbersome and the kspace propagators can solve the problem of being slow.

To calculate the spatial gradient of the pressure, the FD should fit a high order polynomial to get an accurate solution. The PS method fits a Fourier series to all the data and the gradient can be calculated simply as

{ } { } (x) 1 - p F ik F x
, here, F is the Fourier transform and 1 -F is the inverse Fourier transform. And the basis functions in Fourier transform are sinusoidal, so only two points per wavelength are needed to describe a wave. Both PS time-domain (PSTD) and k-space method calculate the spatial gradients using FFT, the PSTD method approximates the temporal gradient with a finite difference, in k-space method, the temporal gradient is solved by a k-space operator [START_REF] Treeby | k-Wave A MATLAB toolbox for the time domain simulation of acoustic wave fields User Manual[END_REF], this operator provide exact solution.

Consider the wave equation in a homogeneous and lossless media:

0 ) , ( 1 ) , ( 2 2 2 2 = ∂ ∂ - ∇ t p t v t p s r r (3. 12)
A pseudospectral solution can be derived by taking spatial Fourier transform and discretising the time derivative using a second-order accurate central difference:

) , ( ) ( ) , ( ) , ( 2 ) , ( ˆ2 0 2 t p k v t t t p t p t t p s k k k k - = ∆ ∆ - + - ∆ + (3. 13) Here 2 2 2 2 z y x k k k k + + = ⋅ = k k ,
where k is the wavenumber.

However, the finite difference approximation of the temporal derivative introduces errors and can only be controlled by limiting the size of the time-step. The k-space methods try to relax this limitation to allow large time-steps to be used. Using an exact solution to the wave equation in the homogeneous and lossless media for an initial pressure distribution (Tabei et Here, u is the acoustic particle velocity, 0 ρ is the ambient density, ρ is the acoustic density.

In k-space method, the computation of the spatial derivatives via the FFT causes waves leaving one side of the domain to reappear at the opposite side (Treeby et al., 2010). To avoid the wave wrapping, a perfectly matched layer (PML) is implemented. A PML is an artificial absorbing layer for wave equations; it strongly absorbs outgoing waves from the interior of a computational region without reflecting them back into the interior. The use of a PML requires the propagating density or pressure to be artificially divided into Cartesian components: The heating function can be modeled as a δ function:

) ( ) ( ) , ( t H t H δ r r = (3. 20)
The initial pressure can be expressed as:

) ( ) ( ) ( ) ( 2 0 r r r H H C v p p s Γ = = β (3. 21)
The wave equation Eq. (3.19) can be solved by Green's function:

' ' ' ' ' 0 ' ' ) ( ) ( ) , ; , ( ) , ( dt d t H t t G C t p V p r r r r r δ β ∫ ∫ ∞ = (3. 22)
The Green's function equation corresponding to Eq. (3.19) is

) ( ) ( 1 ' ' 2 2 2 2 t t t G v G s - - - = ∂ ∂ - ∇ δ δ r r (3. 23)
The Green's function is written as:

' ' ' ' ' 4 )] ( [ ) , , , ( r r r r r r - - - - = π δ t t v t t G s (3. 24)
To get an expression of the Green's function in k-space, we take a 4D Fourier transform of Eq. (3.23):

' ' ) , ( ) , ( ˆ2 2 2 t i i s e e G v G k ω ω ω ω r k k k ⋅ - - = + - (3. 25)
The Green's function in k-space becomes:

k r r r r k d d v k e e t t G s t t i i ω ω π ω ∫ ∫ - = - - - ⋅ 2 2 ) ( ) ( 4 ' ' ) / ( ) 2 ( 1 ) , ; , ( ' ' (3. 26)
Using the property of the δ function

) ( ) ( ) ( 0 ' 0 ' t f dt t f t t - = - ∫ δ , Eq. (3.22) becomes ' ' ' ' ) , ; , ( ) ( ) , ( r r r r d t t t G x H C t p V p ∂ ∂ = ∫ β (3. 27)
We can also express Eq. (3.27) with the initial pressure,

' ' ' 0 2 ) , ; , ( ) ( 1 ) , ( r r r r r d t t t G p v t p V s ∂ ∂ = ∫ (3. 28)
Eq. (3.26) can be solved by using Cauchy's residue theorem (J. D. [START_REF] Jackson | classical electrodynamics[END_REF] :

k r r r r k d e k kt v v G i s s ) ( 3 ' ' ) sin( ) (2 ) ; t , ( - ⋅ ∫ = π (3. 29)
The derivative of G becomes:

Resolution of the acoustic propagation

Let us remind the photoacoustic equation: 3.2 listed the chosen values of the acoustic parameters used for the simulations: Once the fluence map ) (r φ known, the acoustic equation can be solved in FEM or by k-space method.

) (t) ) ( ( ) ( ) ( ) ( ) , ( ) ( 1 ) , (

Resolution of the acoustic propagation in k-space

As the k-space method is developed to solve the acoustic equation, we first have to introduce the initial pressure distribution from FEM to k-space. As the mesh in FEM is irregular, while the mesh in K-space is regular grids, we should interpolate the data in FEM to a regular mesh, in our case, it's a 60*60 grid. Acoustic parameters

1 4 4 - - = K e β 1 _ 1 3 4 K Jkg e C p - = 1 1485 - ⋅ = s m v s
The initial pressure is expressed as:

) ( ) ( ) ( ) ( 0 r r r r φ µ a p Γ =
, shown in Fig. 3456. Firstly we have to define the simulation geometry, then discretize this geometry into regular mesh. The number of the mesh points and the distance between mesh points is defined by us.

Secondly we have to define the acoustic properties of the medium: the sound speed, the mass density, the power law absorption exponent and the power law absorption coefficient. If the absorption parameters are not set, the simulation is assumed to be lossless.

Then we have to define the acoustic source term in the medium. One can choose between three types of sources. The first one is the initial pressure distribution, this kind of source is used to simulate pulsed photoacoustic tomography, which is our case. But it can also account for time varying pressure sources, corresponding to a mass source, appearing as a source term in the mass conservation equation. The third kind of sources is a time varying particle velocity source, corresponding to a force source term, in the momentum conservation equation.

The last parameter to define is the acoustic sensors that are used to record the acoustic signals. The positions of the sensors have to be defined.

Once these parameters defined, the simulation is ready for run. Fig. 34567shows a comparison between the photoacoustic pressure solved by FEM and k-space respectively measured at different detectors positions D1, D2, D3 and D4. The acoustic pressures solved by k-space method are more stable, have less oscillations; the accuracy of FEM is mesh-dependent, when solving large time during of the PDE, the FEM becomes very slow, from the comments of last paragraph, we know that k-space method is much faster than FEM method, so finally we used k-space method to solve the acoustic propagation.

In k-space, the computation grid in 60*60, the dimension of each grid is 1mm*1mm. The acoustic detectors are focused at the centre of the breast with a focused distance of 25 mm, both the number and the length of the detectors are variable. In the K-space model, the breast can furthermore be treated as an acoustically absorbing medium, the absorption coefficient follows a frequency power law of the form (Treeby and Cox, 2010):

y ω α α 0 =
. Here, 0 α is the absorption coefficient ω is the angular frequency, and y is the power law exponent. In our case,

cm MHz dB / / 75 . 0 0 = α and 5 . 1 = y
, which is corresponding to the absorption parameters in the breast.

Detectors geometry

In the k-space simulations we considered, the acoustic detectors are located at the periphery of the breast as shown in Fig. 345678, both the length and the number of the detectors can be vary , the effect of the distribution of the detectors are further discussed in the next chapter. -Each detector is composed of elementary points (a pixel), the number of the points is defined by the user, the choice of the number of points is discussed in Chapter 4.

-The pressure generated by these points is collected at each elementary point and averaged over the total length of the detector. Hence, for each detector, there is only one single value output, that is the averaged data as shown in Fig. 3456789.

Fig. 3-9 each detector is composed of several points, the output of the points are averaged, so for each detector, there is only one output.

Then we will show examples of acoustic pressure collected by two kinds of detectors.

Sound speed inhomogeneity

In the above work described, we assumed that the acoustic speed to be constant, 1485 m s -1 everywhere. In k-space, both the inhomogeneity of sound speed and density are considered. We can rewrite the equations as: Here, S is the source term, for both the sound speed s v and density ρ , instead of using the homogeneous one s v 0 and 0 ρ , inhomogeneity can be considered in s v and ρ . Now we showed an example of the photoacoustic pressure considering the sound speed inhomogeneity. Here the acoustic speed is supposed to be inhomogeneous, 1485 m s -1 in water, 1540 m s -1 in the breast, 1620 m s -1 in the tumour. [START_REF] Goss | Comprehensive compilation of empirical ultrasonic properties of mammalian tissues[END_REF] tumour 1570-1620 [START_REF] Kossoff | Average velocity of ultrasound in the human female breast[END_REF] Fig. 3-14 shows the acoustic pressure computed by considering inhomogeneous acoustic speed. There are two N-shape signals when the acoustic speed is not constant, the former pressure is generated by the tumour, the latter pressure is transferred from source to tumour, then from the sound is not constant, the Huygens' Principle does not hold anymore. However, there is still some decay in the bounded domain. To avoid inverse crime, that is, when the same meshes are used in the forward and inverse models, one should choose a mesh in the inverse model slightly different from the mesh in the forward model.

αu u - ∇ - = ∂ ∂ ρ p t
In the forward model, the computation area is 60 mm* 60 mm. We used the same computation area in the inverse model, but the computation grid is 59*59. In the inverse model, the number of detectors is the same as in the forward model.

Quantitative photoacoustic tomography (QPAT) (active)

The initial pressure map:

)) ( ), ( , ( ) ( ) ( ) ( 0 r D r r r r r p a a µ φ µ Γ =
is the product of the Grüneisen coefficient, the optical absorption coefficient and the optical fluence which depends itself on the local optical properties of the medium, in terms of absorption and diffusion. Hence, in order to go further in the quantification through an accurate estimation of the absorption coefficient, uncoupling the effects of absorption and diffusion is necessary. This is the main goal of quantitative photoacoustic tomography (QPAT: complete uncouple the maps of the absorption, scattering and Grüneisen parameters).

We tried two approaches to uncouple the absorption and diffusion coefficient, one approach is to recover these two parameters from the reconstructed initial pressure map; the other is to recover the two parameters from synthetic measured photoacoustic pressure. It's not possible to recover two parameters simultaneously from one single measurement. To recover both, we have to use multiple optical sources or multiple wavelength illuminations. In our work, we used multiple optical sources, because in our experiment setup, only one wavelength optical source is available.

Recovering the optical parameters from the reconstructed initial pressure map

First of all, we express the light fluence considering the perturbation of absorption and diffusion: , we have that [START_REF] Wang | Biomedical Opticals: Principles and Imaging[END_REF]:

' ' ' 0 0 ' ' ' ' 0 0 ' ) , ( ) , ( ) ( ) , ( ) , ( ) ( ) , ( r r r r r r r r r r r r r r d G D D d G D s j s j a s j SC φ δ φ δµ φ ∇ ⋅ ∇ + - = ∫ ∫ (3. 34) Here 0 D is the diffusion coefficient of the background ) , ( ' 0 r r l G
is the Green's function representing propagation from ' r r → j . We can discretize the integral as , we take a ratio of reconstructed initial pressure map at optical source location i and n : We suppose that the light fluence distribution is slowly varying, the absorption and diffusion Eq. (3.42) can be written in matrix form:

) ( ) ( ) , ( ' } { ' } { l l s l ij l a l a l ij s j SC D W W r r r r δ δµ φ ∑ ∑ + = (3. 35) Here, 0 ' 0 ' 0 } { / ) , ( ) , ( D V G W s l l j a l ij ∆ = r r r r φ 0 ' 0 ' 0 } { / ) , ( ) , ( D V G W s l l j s l ij ∆ ∇ ⋅ ∇ = r r
)) ( ) ( ) , ( (
)) ( ) ( ) , (
perturbation is small, if H is linear, ≅ ) , ( ˆr r j p )) ( ) ( ( ) , (
b u Q = δ (3. 42)
Matrix Q is usually not invertible, and it's not injective, there might be more than one solution that satisfies equation (3.43).If the pseudoinverse of Q exists (inverse of the Hessian matrix

Q Q H T = ),
the unknown perturbation can be estimated by Moore-Penrose pseudoinversion

b Q Q Q u T T 1 ) ( - = δ .
Alternatively, we may introduce a singular value decomposition to decompose the Hessian matrix

T T V U Q Q H ∑ = =
, the inversion is given as

i i T i i T v b u b U V u σ δ ∑ = ∑ = + , +
∑ is the pseudoinverse of ∑ . When the pseudoinverse of Q doesn't exist, Tikhonov regularitation or other regularization have to be used to avoid instabilities due to poor matrix conditioning.

The inversion is definitely affected by the accuracy of the reconstruction of the initial pressure distribution, so we sought for a second approach: recover the optical parameters from the measured photoacoustic pressure directly.

Recovering the optical parameters from measured photoacoustic pressure

In the second approach, we rewrite the distribution of the initial pressure as:

) , ( ) ( ) , ( 0 u r r u r φ µ a p Γ = , T a ] [ D µ u = .
The photoacoustic pressure ) , , ( t p u r then propagates at the speed of sound in the medium s v as:

)) ( ( ) , ( ) , , ( 0 2 2 2 2 t t p t p v t s δ ∂ ∂ =         ∇ - ∂ ∂ u r u r (3. 43)
which solution in free space is well known and can be expressed as [START_REF] Li | Photoacoustic tomography and sensing in biomedicine[END_REF]:

)] ( ) , ( 1 [ 4 1 ) , , ( ' ' 0 ' 2 s s s v t p dr t v t v t p r r u r u r - - ∂ ∂ = ∫ δ π (3. 44)
To avoid the calculation of time derivative, the velocity potential v p is also frequently used: , becomes [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]:

) ( 4 ) , ( 1
∑ = ∆ - Γ = L l l s l s dl dl d s V E v R t R t p 1 model ) ( 4 ) , , ( τ τ δ π u r (3. 47)
And the Jacobian matrix:

∑ = = = ∂ ∂ = L l j s l dl d s j s d t p 1 0 , , model 0 )) ( ( ) , , ( )) ( ( 0 r u u u r r u J u u β α τ τ τ (3. 48) ) ( 4 , s dl dl l dl v R t R V - Γ∆ = τ τ δ π α
is a time evolution transfer matrix propagating the pressure created at initial time at position l r to point detector located at d r ;

)) ( ( 0 , j s l r u β represents the sensitivity of the energy instantaneously deposited due to the optical perturbation u δ located at positions j r , j r belonging to the reconstruction mesh: (3. 49) lj δ is the Kronecker symbol. Invoking again the perturbation approach, knowing the Green's function

              ∂ ∂ ∂ ∂ + =         ∂ ∂ ∂ ∂ = = = = 0 0 0 ) ( ) , ( ) ( ) ( ) , ( ) ( ) , ( ) ( )
0 G of the diffusion equation ( ) ( )) ( ) ( .( ) ( ) ( r r r r r Q D a = ∇ ∇ - φ φ µ
), that can be computed for the specified geometry and optical properties of the reference medium: ).

) ( ) ( ) , ( ) ( ) ( ) , ( ) , ( ) , ( 0 0 0 0 j j l s j a j a l s l s l s D D r r u r r r u r u r u u u r u u u u δ φ δµ µ φ φ δ φ = = ∂ ∂ + ∂ ∂ + ≈ + = (3. 50)        ∆ - ∇ ⋅ - ∇ = ∂ ∂ ∆ - - - = ∂ ∂ = = 0 0 0 0 0 0 0 0 / )) ( , ( )) ( , ( ) ( ) , ( / )) ( , ( )) ( , ( ) ( )

Optical energy deposition matrix:

)) ( ( 0

, j s l r u β
The computation of this matrix is obtained through the calculation of several matrices involving the resolution of sets of diffusion equations:

Resolution of the inverse problem

Different reconstruction methods were tested: Singular value decomposition (SVD) method for which the Hessian matrix was decomposed into:

T V S U H × × =
. S is a diagonal matrix with same dimension as H containing nonnegative singular values i σ ordered in decreasing order, U and V are unitary matrices. The inversion of )

2 2 λ σ σ + i i
, λ is the Tikhonov regularization parameter.

The LSQR algorithm [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF], a classical iterative algorithm for sparse linear equations and sparse least squares, that searches for the minimization of the L2-norm of The algebraic reconstruction technique (ART) [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF] that searches iteratively for the solution u δ such that:

( ) [ ] H H u H b λ u u 2 1 , - + 〉 〈 - + = n n n n δ δ δ , n
λ is a hyperparameter that can be tuned in order to regularize or to speed up the reconstruction. In all what follows, n λ was kept constant and equal to 1. Non negative values constraint is also introduced in the algorithm. These reconstruction methods are tested for different noise levels. For small noise levels, SVD provides sufficiently good results but is known to be unstable and very sensitive to noise, while the ART behaves better, especially when constrained.

In the second approach, as we used the photoacoustic pressure directly to recover the optical parameters, this avoid the inaccuracy brought by the reconstruction of the initial pressure map, but the measured photoacoustic pressure is time-varying, we have to sample the pressure at different time steps, this is a burden to the computation. The reconstructions will be shown in the next chapter.

Conclusion

In this chapter, the implementation of the photoacoustic tomography model was discussed.

This model refers to both the forward and the inverse models. The forward model refers to as the simulation of the generation and propagation of photoacoustic pressures, while the inverse model aims to reconstruct the initial pressure map from the collected photoacoustic pressure in the forward model and recover of the absorption and diffusion coefficients.

The optical model is solved by the FEM while the acoustic model is solved by k-space method. The coupling of the FEM and the k-space method is proven efficient to solve the forward problem.

For the inverse model, the time reversal method in k-space is used. Then we introduced the implementation of the QPAT. We tried two approaches; one is to reconstruct the absorption and diffusion coefficients from reconstructed initial pressure map, the second is to reconstruct the two parameters from synthetic measured photoacoustic pressure, as the efficiency of the first approach is affected by the accuracy of the reconstruction of the initial pressure map, we finally choose the second approach.

Introduction

In this chapter, the theoretical tools presented in the previous chapter are employed and two simulation studies are presented:

The first one is related to the optimization of the acoustic detectors distribution in photoacoustic tomography (PAT). Through 2D simulations, we study the influence of the sizes and angular distribution of the detectors on the quality of the images of the initial pressure distribution.

According to the size of the object to be reconstructed (a tumour), the distribution of detectors have to be adapted. The acoustic detectors are used to collect the acoustic signals, these acoustic signals are then used to reconstruct the initial pressure map, to locate the abnormality inside the tissue. To optimize the detectors distribution, we varied the size of the detector, the number of the detectors, the angular coverage of the detectors (full coverage or sparse coverage), and included the limited frequency bandwidth of the detectors. We also optimized the sources distribution in this part.

The second one is related to the optimization of the optical sources and the acoustic detectors distribution on quantitative photoacoustic tomography (QPAT) and their influence on the reconstruction of the absorption and diffusion coefficients maps. Quantification in photoacoustic imaging has proven to be challenging. The multiple optical illumination sources approach is used to recover the optical parameters.

Effect of the detectors distribution on PAT

Effect of the detectors distribution to the reconstruction of point objects

The acoustic detectors we used in the experiment are immersion transducers produced by Olympus Company as showed in Fig. (4.9). An immersion transducer is a single element longitudinal wave transducer with a 1/4 wavelength layer acoustically matched to water. Focusing lens can be added to increase transducer sensitivity and performance in a particular area of a part. Immersion transducers can be focused either spherically or cylindrically upon request within the allowable focal range for a given frequency and element size. In our case, the detectors we used are focused detectors with element diameter of 6 mm and 10 mm respectively with central frequency of 3.5

MHz. To study the detectors distribution, in 2D, the detectors are simulated as arcs, in 3D, the detectors are simulated as discs.

To evaluate the effect of the detectors distributions, we firstly created an initial pressure map composed of 5 point sources, with equal amplitude 1, then used the acoustic pressure collected by these detectors to reconstruct the initial pressure map. We evaluated the distributions of the detectors through the quality of the reconstructed amplitude of the points. Here, all the simulations were done based on k-space method for the forward model and Time Reversal Method for the inverse problem. Fig. 4-1 shows the distribution of detectors and objects. To avoid the inverse crime, we used two different meshes in the forward and inverse models resolution. For both models, the computation size is the same, 47.2 mm. For the forward model, there are 472*472 computation grids, so each grid has a dimension of 0.1 mm * 0.1 mm; while in the inverse model, only 462*462 grids. The detectors are supposed to be focused; they focused at the centre of the geometry. They are modeled by considering them as an arc, the radius and the center of the arc is the same as the simulation geometry. In the forward model, each detector was composed of several points elements, the points are evenly distributed along the arc. The data collected by these points is averaged, and assigned to each point to do the reconstruction.

We borrowed the idea of the detectors distribution from [START_REF] Allen | Optimising the detection parameters for deep-tissue photoacoustic imaging[END_REF]. Two kinds of detectors are considered: intense and sparse detectors distributions. For the sparse distribution, the number of detectors was constant and the length of each detector was varied. 45 detectors are used here. For the intense detectors distribution, the total length of the detectors covers angularly the object and the number of detectors varies inversely with the length of the detectors. Fig. 4-2

shows the distribution of these two kinds of detectors. The focusing of the detectors influences the amplitude of the reconstructed values in the reconstruction maps. The purpose of this study is to optimize the detectors distribution in order to obtain the best reconstructed maps.

Optimization of the number of point elements

The detectors in the forward model are composed of points. First, we searched for the necessary number of points to model a detector of a given length. For two very different lengths of detectors, 0.8 mm and 2.4 mm, the number of elementary element points constituting the detector is increased from 2 to 32, the initial pressure points map are reconstructed. From Fig. 4-3 we can see that, for all kind of detectors, even with large length, 8 points per detector are sufficient, taking more than 8 do not modify visibly the reconstructions. In the following simulations, for all the detectors even with large length, we used 8 points in the forward model.

Effect of the length of the detectors

The length of the detectors was varied from 0.1 mm to 2.4 mm. Each of the 5 points should be reconstructed with the same amplitude 1. Fig. 4-4 shows the effect of the length of the detectors to the reconstructed initial pressure map. For the sparse detector, the reconstructed initial pressure increases as the length of the detector increases and it goes with a higher discrepancy between the According to these results, for sparse detectors distribution, the larger the detector the less discrepancy between the values reconstructed. In the intense distribution case, the reconstructed values are always higher than in the previous case, the smaller the detector the less discrepancy between the reconstructed values. The optimized length was found to be 0.2 mm. These figures show the detector needs to be as small as possible, even for sparse detectors geometry.

Effect of the detectors distribution as a function of the size of the abnormality in 2D.

In this part, we tested the effect of detectors distribution considering our breast geometry shown in Fig. 4-6, with a tumour whose size is varying. We also considered two kinds of detectors:

sparse and intense detector distribution. For the sparse distribution, there are 22 detectors. The radius of the tumour varied: 0.5mm, 1.5mm, 2.5mm, 5mm, and 7.5mm. The quality of the reconstructed images is evaluated through the value of the quadratic error:

∑ - = j i ij reij ij p p p d , 2 0 2 0 ) ( , 0
p is the initial pressure, re p is the reconstructed initial pressure. The image obtained (Fig. 45678) by intense distribution has better quality (d closer to 0) than by sparse detector; for sparse detector, as the detector length increases, the image quality increases; while for intense detector, as the detector length increases, the image quality decreases.

The second remark is the larger the tumour to be reconstructed, the better the results for both configuration. In this evaluation, there are two variables: the length of the detector and the radius of the tumour, whatever the length of the detector, as the radius of the tumour increases, the quadratic error follows the same trend, so we can evaluate the image with two variables.

Besides quadratic error, there are other criteria to evaluate the reconstructed image, such as contrast, resolution. We may extend our study at this approach. Fig. 4-9, the cube with the edge-length of 50mm was used to simulate the breast. The focused detectors located around the breast were selected based on the detectors we used in the experiment. The focus distance of the detectors are 25mm, the radius of the detectors varied at 2mm, 3mm, 4mm, 5mm, respectively. We put the tumour at the centre of the breast, which is the best situation. For the sparse detectors, the number of the detectors is constant (4 detectors) and only the radius of each detector is varied; for the intense detectors, the number and the radius of the single detector are varied, but the product of the number and the radius of the detector is constant. The height of the tumour is 50mm, the radius of the tumour varied at 1mm, 2mm, 3mm, 5mm, and 7mm, respectively.

In this work, sparse and intense detectors distributions were considered, the radius of the detectors and the radius of the tumour were varied for each situation and the quadratic error at the middle plane of x-y direction, x-z direction and y-z direction between true and reconstructed images were calculated and reported in Figure 4-10. mm, the radius of the tumour is 7mm; in(c), there are 26 detectors, the radius of each detector is 3 mm, the radius of the tumour is 1 mm, in(d), there are 16 detectors, the radius of each detector is 5 mm, the radius of the tumour is 7 mm.

In both 2D and 3D, intense detectors always behave better than sparse detectors. For intense detectors, in 2D, when the length of the detector increases, the reconstructed image quality get worse; while in 3D, when the length of the detector increases, the reconstructed image quality get better, that is because in 3D, when the length of the detector increases, the total coverage surface of the detector increase.

Effect of the detectors frequency bandwidth in 2D

To -40 f min f max detectable frequency decreases; as the source radius increases, the maximum detectable frequency decreases.

Fig. 4-13 maximum frequencies at different distance with source radius of 0.1 mm, 0.2 mm, 0.3mm, 0.5 mm, and 1 mm, respectively.

Effect of the optical sources distribution in 2D

In this part, we tested the influence of the number of the optical point sources. Fig. 4-14 is the simulation geometry, the detectors located just at the periphery of the breast. In the following simulations, the number of detectors is 135, and the length of each detector is 0.8 mm.

Optimization of the illumination/detection protocol in QPAT

In this study, the multiple illumination approach described in Chapter 3 is applied in order to test the robustness of the approach and to define the more suitable illumination/detection protocol for uncoupling absorption and diffusion coefficients uncoupled maps.

Sets of synthetic data were generated by solving the forward problem and introducing a fraction of random noise with normal distribution. Hereafter are presented the results of simulations performed on the synthetic phantom with different reconstruction methods. The purpose of the study is to highlight the effects of different sources and detectors distributions and geometries on the reconstruction of the perturbation of the optical parameters

T a ] [ D µ u δ δ δ =
, in terms of accuracy and robustness to noise. To that end, the following aspects were considered: -number of point sources distribution: the point source illumination has been specifically studied in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF], here the results are shown as reference;

-point detectors distributions: the purpose is to understand which transducers distribution would be the best, for a given number of measurements;

-using extended sources instead of the point sources: point sources have been studied but, in practice, extended sources are easier to handle experimentally;

-limited angle examination: this corresponds to actual experimental situations, here the transmission geometry, where sources and detectors belong to different half-spaces, is studied.

Our simulations were conducted in 2D. The chosen test object (Fig. 4-17) was basically the same as in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]): a 6 cm × 6 cm square, in the centre area of which a 2 cm × 2 cm square was selected as the reconstruction area (Fig. [4][5][6][7][8][9][10][11][12][13][14][15][16][17]. The background absorption and diffusion coefficients are respectively , positions: (0.6cm;1cm) and (1.6cm;1cm) , dimensions: 0.3cm×1.1cm and 0.3cm×0.3cm) and one diffuser ( D δ , position : (1cm;1cm), dimensions: 0.5cm×0.5cm).

Different reconstruction methods

Four point sources and sixty detectors (15 detectors on each side of the object) are considered for this study. A first remark that can be done is the poor conditioning of the Hessian:

after performing the SVD, and examining the spectrum of the singular values i σ of the Hessian matrix H(in Eq. (3.47)): the condition number obtained for this configuration is the introduction of non-negative constraint, ART seems to have better potential than the two others in the recovery of the optical parameters at higher noise levels. It is interesting to notice that, though ART reconstructs the perturbations with a higher contrast, LSQR seems to provide reconstructions that look stable as the noise level increases. In all what follows, ART was chosen as the reconstruction algorithm.

Hence, the aim of the following studies is to improve the conditioning of H that the

Sources distributions

In this study, 60 point detectors were distributed around the object and the sources number was varied (2, 4, 8 and 16 sources considered). Examining the condition number of H as a function of the number of sources (Fig. [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], the reconstructions should be indeed more stable as the number of sources increases. Although the inversion is mathematically demonstrated, two sources are not enough to obtain robust reconstructions. Fig. 4-21 shows the reconstructions obtained: with 2 sources the diffusion coefficient perturbation is almost not detected, but interestingly, the absorption perturbations are well recovered. An improvement of a factor 10 in the condition number is obtained by adding only two sources. However, the improvement in increasing the number of sources might be reduced after 8 sources and the gain in the reconstructions might be small even if the number of measurements (and hence the duration of the examination) is increased. This is all-the-more visible in the cross-plots shown in Fig. 4

Example of simulations for experimental situations with limited angle examination

Examining the object under all angles may be experimentally complicated or simply not feasible. Here are reported reconstructions obtained under constrained examination geometry: the medium is probed exclusively under transmission geometry, sources and detectors can rotate and probe the object on its four sides. A first remark is the reconstructions obtained with point or line sources illuminations (Fig. 4-30 (a) and (b)) are improved compared to those for which the detectors are collecting signal from all the sides of the object: when the object is probed in transmission, although the condition numbers are several orders of magnitude higher, still one can get reconstructions of satisfying quality with illumination scheme seems to play a major role in the quantification, while the detection in the localization with high resolved time-of-flight measurement. An important result is that probing the sample under transmittance geometry only provides better reconstructions than with full angular coverage detection.

-Different simple experimental situations were tested and good quality reconstructions were obtained with a number of measurements extremely reduced.

Experimental validations may support these findings that may influence future experimental setups designs. The present study has highlighted the importance of the data acquisition protocol. A smart combination between different types of measurements will undoubtedly improve the quality of the reconstructions. An example of such a process would be to combine one measurement taken under a full angle illumination, with a field of illumination as wide as possible, and a set of data acquired under different illumination angles. The first measurement would serve in a first reconstruction loop that would provide initialization and constraint values for a second reconstruction loop, performed with the second series of data.

Conclusion

In this chapter, we presented two studies on the optimization of the sources and detectors distribution in PAT and QPAT.

In PAT, we studied two kinds of detectors distribution: sparse (number of the detectors constant length of each detector varied) and intense (total length of the detectors maximum and number of detectors varied inversely with the length of the detectors) detectors distribution. The intense distribution has better performance than sparse distribution because intense detectors usually have more angular coverage surface than sparse detectors. In 2D, as the single element size of the detector increases, the sparse distribution behaves better, while the intense distribution behaves worse. In 3D, the results are slightly different: as the single element size of the detector increases, both the sparse and the intense distribution behave better. In terms of frequency bandwidth, the maximum detectable frequency is inversely proportional to the detection distance, and is inversely proportional to the source radius. With more evenly distributed optical sources, we can get a more homogeneous reconstructed initial pressure map. This study has been presented at the 1 st EOS Topical Meeting on Frontiers in Optical Imaging (FOI 2013).

We used synthetic photoacoustic pressure to test the MI-PAT (Multiple Illumination PAT).

The results of the study show multiple and wide field illumination, with large angular coverage of the acoustic detectors, can better recover the perturbations of the absorption and diffusion coefficient. We also proved that, with a very restricted sources and detectors distribution as in our setup, we can still do the QPAT and obtain uncoupled reconstructions of the absorption and diffusion coefficients. This study has been reported in a peer review paper (submitted 17/06/2014, to Biomedical Optics Express).

Chapter 5. Modeling the experimental situations This chapter is related to photoacoustic tomography simulations, in soft tissues, such as the breast. The experimental setup developed at Institute Fresnel is firstly presented and simulations are compared to experimental data. In the last part simulations in QPAT are presented. They are related to realistic situations that could be reproduced with this experimental setup. In our system, the laser is Nd: Yag (Quantel Eazy) laser with pulse of 5 ns at repeat frequency of 10 Hz at wavelength of 1064 nm, the pulse energy is 330 mJ, after the laser pulse is generated by the laser device, it was transferred to the breast phantom via optically system as Shown in Fig. 5-3.

IF experimental set up

The acoustic transducer is a focused Panasonic transducer with focus radius of 2.5 cm and the central frequency is 3.5 MHz, the breast phantom is a cylinder with radius of 2.5 cm. The phantom composed of a mixture of Agar, Tio 2 powders water and black ink. The position of the source is fixed, the phantom is rotated as shown in Fig. 5-5, the laser was reduced to a laser beam, the diameter of the laser beam illuminated on the surface of the phantom is about 0.5 cm. During the tomography measurement, the position of the source and detector is fixed, and only the phantom is rotated every 10 degree. Table 5-1 is the parameters we used in the simulation. As the abnormality in the phantom is almost opaque, in the simulation the absorption coefficient of the tumour is much higher than the normal breast. We modulate the parameters to match with the experiment results. From Fig. 56789, we can see that the noise of the experimental data is mainly the high frequency oscillations, so we fit an 8 order Butterworth low pass filter with the cutoff frequency of 10 MHz. Whatever the rotation angle is, the distance between the optical source and the detector is constant, we can also see it from Fig. 5-13, the N-shape pressure located at about 6 cm is generated from the source, the distance indicates the distance between the source and the detector. In Fig. 5-13 (a), the tumour is close to the source, the pressure generated by the tumour is much higher than by the source, the big N-shape pressure is from the tumour. In most of the cases, as the phantom is far from the source, the pressure generated by the tumour is hardly visible. Because of our experimental setup, the rotation center is not exactly the center of the phantom, there are some mismatches between the experimental data and the simulation.

QPAT simulation related to IF experimental set up

In this part, we mainly discussed if it's possible to do the quantitative reconstruction with the existing experimental setup. It's well known that with a single illumination, it's not possible to reconstruct the absorption and diffusion coefficient simultaneously. The problem becomes how to rotate the phantom, to get the multiple measurements. Rotating the phantom is the same effect as rotating the source-detector pair. Fig. 5-14 gives the 2D illustration of the problem. The breast has a radius of 26 mm, the 4 arc sources located at a circle of radius 25 mm, we introduced a 1 mm extrapolated boundary from the arc sources, the Dirichlet boundary condition is applied of the extrapolated boundary(M. [START_REF] Schweiger | The finite element method for the propagation of light in scattering media: Boundary and source conditions[END_REF]. we selected a computation area of 25mm*25mm (25*25 grids) to do the QPAT, there are two perturbation inside the computation area, the absorption perturbation located at (4mm, 4mm), has a radius of 2 mm, the diffusion perturbation located at (-4mm, 4mm), has a radius of 2 mm, too, there are 32 point detectors located at a square with a side length of 61 mm, there are 8 point detectors at each side, outside the breast is water, the geometry is a square with a side length of 100 mm. Table 5.2 lists the optical parameters used in the simulation, here, the absorption perturbation is 0.01 cm -1 , and the scattering perturbation is 1 cm -1 , so the diffusion perturbation is 0.003 cm.

Conclusion

In this chapter, we introduced our experimental setup, and compared the photoacoustic pressure generated by the experiments and simulations. However, the measurement of the optically properties could not be conducted for 1064 nm. This comparison highlighted the difficulties when facing experimental measurements (misalignment, mismatch with optically and acoustic properties).We also simulated the QPAT under experimental setup, it's proven that it's possible to do the QPAT, what's interesting is to use our experimental data to do the QPAT directly. With very limited number of measurements, with sources and detectors positioned as proposed here, information on the position and quantification of the abnormalities can be obtained. Due to the lack of time, the study is incomplete. Additional measurements and simulations have to be done. This PhD work was devoted to the implementation and study of a new numerical method for solving the forward and inverse problems in Quantitative PhotoAcoustic Tomography, especially for tissues physical properties and geometries similar to those of the breast. The method relies first in a specific illumination/detection protocol: it consists in interrogating the medium with multiple illumination sources. By doing so, in a way that has some similarities with Diffuse Optical Tomography, it is demonstrated that it is possible to achieve the uncoupling of the optical properties reconstructions.

The physical bases of this multi-physics phenomenon are introduced first. The biological tissue is illuminated with a pulsed source with wavelength belonging to the so-called "therapeutic window" (600-1000 nm), where the tissues are less absorbing but highly diffusing. A diffusing wave is generated and propagated through the medium, at the speed of light, with a diffusion equation. The energy deposited is dissipated mostly into heat that produces locally an adiabatic expansion. An ultrasonic pressure wave is generated. This initial pressure wave then propagates according to a wave equation at the speed of sound in the medium. Speeds of sound and light are so different that the heating step can be considered as instantaneous from the point of view of the acoustic detection. The collected signal vehicles both the optical and the acoustic properties of the medium.

More specifically, the initial pressure distribution map hence contains the information related to the optical properties of the medium. The main difficulty relies in uncoupling the absorption and diffusion attenuations.

Within this work, the k-space method introduced by (Cox and Beard, 2005) [START_REF] Cox | k-space propagation models for acoustically heterogeneous media: Application to biomedical photoacoustics[END_REF] for solving the acoustic propagation with a given initial pressure distribution map has been coupled

with the resolution of the diffusion optical equation by the Finite Element Method (FEM). This forward solver is then used to generate synthetic data that are eventually blurred with controlled amount of noise. The simulations were conducted by considering experimental situations when possible. To this purpose, focused detectors are modeled, disks (2D) or cylindrical (3D) phantoms are efficient data acquisition protocols (multiple and wide field illumination), and we believe that they will serve for improving the design of the setup.

Chapter 1 .

 1 Introduction to photoacoustic imaging.............................................................. photoacoustic mammography ........................................... 26 1.4.1. The laser photoacoustic imaging system (LOIS) .........................................

Fig. 1

 1 Fig.1-1(Song et al., 2008) Non-invasive in vivo photoacoustic images of the SLN in a rat. (A) Photograph with hair removed before photoacoustic imaging. (B) Photograph with skin removed after photoacoustic imaging. (C) Control photoacoustic image without methylene blue injection. Bright parts represent optical absorption, blood vessels. x and y denote B-scan and 3D imaging. (D) After-injection photoacoustic image.

Fig. 1

 1 Fig.1-2[START_REF] Xia | Photoacoustic Tomography: Principles and Advances (Invited Review)[END_REF] schema of OR-PAM and AR-PAM (a) Second-generation optical-resolution photoacoustic microscopy system (2G-OR-PAM)[START_REF] Hu | Second-generation opticalal-resolution photoacoustic microscopy with improved sensitivity and speed[END_REF], employs a wavelength-tunable laser system; (b) Dark-field acousticresolution photoacoustic microscopy (AR-PAM)[START_REF] Maslov | In vivo dark-field reflection-mode photoacoustic microscopy[END_REF], A Q-switched pulsed Nd:YAG laser operating at 532 nm, delivers 300 mJ of energy per pulse to a 0.60-mm-diameter optical fiber. The laser pulse width is 6.5 ns, and the pulse repetition rate is 10 Hz.

Fig. 1

 1 Fig.1-3(Maslov et al., 2005) Vascular distribution in rat skin. (a) Four in situ consecutive photoacoustic B scans: 1, 2, vessels perpendicular to the image plane; 3, in-plane vessel; (b)the in situ photoacoustic projection image taken from the epidermal side; (c) a photograph taken from the dermal side with transmission illumination; (d) In vivo noninvasive photoacoustic projection image taken from the epidermal side.

Fig. 1 - 4 (

 14 Fig.1-4(Maslov et al., 2008) microvasculature in a nude mouse ear obtained by OR-PAM

Fig. 1

 1 Fig.1-5 (Zhang et al., 2007)PAM imaging of variations in SO2 in single subcutaneous blood vessels in a 200 g Sprague Dawley rat in vivo. (a) Structural image reflecting the total hemoglobin concentration acquired at the 584 nm optical wavelength under hyperoxia. (b) Static SO2 image within the marked region in panel A under normoxia, where arteries and veins are pseudocolored red and blue, respectively, based on the imaged SO2 values. (c) Image of the SO2 changes from normoxia to hypoxia (hypoxia value-normoxia value). (d) Image of SO2 changes from normoxia to hyperoxia (hyperoxia value-normoxia value). (e) Typical imaged values of SO2 in venous and arterial bloods under all three physiological states, where different trends of variation are observed.

  Fig.1-6 example of different kinds of transducers

Fig. 1

 1 Fig.1-7 (Ermilov et al., 2009) Schematic diagram of LOIS-64

Fig. 1

 1 Fig.1-8 (Ermilov et al., 2009)(a) Digital mammography image (mediolateral projection) of the breast, showing a round, apparently benign tumour. (b) Ultrasonic imaging with the Doppler showing the tumour boundaries and the focal blood flow increase. (c) Photoacoustic image (mediolateral projection) showing the enhanced contrast around the area of the focalblood flow increase detected by the Doppler ultrasound.

Fig. 1

 1 Fig.1-9 (S. Manohar et al., 2005) Schematic of the Photoacoustic Mammoscope

  Caucasian woman. The X-ray mammogram and breast sonogram found an area of suspect as shown in Fig.1-10 (a) the area within the red dot line, and then the photoacoustic mammography was performed on this area, where they found the vascular signature of the cancer due to higher absorption in the breast shown in Fig.1-11 (black dot line).

  Fig.1-10 (Vaartjes et al., 2007)54 year old woman with tumour in right breast: (a) x-ray mammogram (b) sonographic scan image

Fig. 2

 2 Fig.2-1 the absorption coefficient of the main absorbers in tissue (http://omlc.ogi.edu)

  Scattering properties are classified into three categories defined by the size of the scattering object relative to the wavelength (Vo-Dinh, 2003): (1) Rayleigh Scattering, where the size of the particles are small compared to the wavelength of the light; (2) Mie Scattering, where the size of the particles are comparable to the wavelength of the light; (3) Geometrical Optics, where the size of the particles are large compared to the wavelength of the light.For the Rayleigh scattering, the particles have the same possibility to scatter forward or backward, for Mie scattering and Geometrical Optical, the particles mainly scatter forward. phase function, is defined to describe the scattering probability of a photon in direction'

  Consider a cylindrical volume element as showed in Fig.2-2, here, sd is the length of the cylinder along the photon propagation direction s , the area dA is perpendicular to the direction s , the solid angle ω d is considered in order to evaluate the changes in photon energy.

Fig. 2 -

 2 Fig.2-2 schematic of a cylindrical volume element

)N

  Contribution 2: ExtinctionEnergy loss per unit time in the volume within the solid angle due to absorption and scattering is given by: is the number density of the scatters and s σ is the scattering cross section of a scatter, dV N s is the number of scatters in the volume element; of light with propagation direction ' s being scattered into ω d around direction s .The phase function depends only on the angle between the scattered and incident directions:

Fig. 2

 2 Fig.2-3 (Prahl et al., 1989) flow chart of Monte Carlo method solving radiative transport problem

  Fig.2-4(a) the photoacoustic pressure detected by a point detector at time t from PA source over a spherical shell centered at the detector position and a radius of vst, (b) an example of photoacoustic pressure from a uniform spherical source at a distance of 20 mm with a radius of 2 mm, the sound speed is 1500m s -1 .

Fig. 2

 2 Fig.2-6 influence of radius to photoacoustic pressure from a uniform spherical source at a distance of 20 mm with radius of 2 mm and 1mm, respectively, the sound speed is 1500 m s -1 .

Fig. 2 -

 2 Fig.2-7 influence of distance to photoacoustic pressure from a uniform spherical source at a distance of 20 mm and 40 mm respectively with radius of 2 mm, the sound speed is 1500 m s -1 .

  The inverse problem in photacoustic tomography firstly refers to as the reconstruction of the initial pressure map or heating function from the collected photoacoustic pressure. The arrival time, duration and amplitude of the pressure indicates the distance from the source to the detector, the size and the strength of the source, respectively. These parameters permit reconstruction of the source geometry from signals measured in different positions. Fig.2-8gives an example to locate the source position.

Fig. 2 -

 2 Fig.2-8 there are two detectors collected the pressure, the pressure arrival time defines the radius of a sphere where PA source is located. The crossing point of the two spheres has the highest probability of the PA source presence.

  photoacoustic resolution is the coupling of the optical propagation and acoustic propagation. The inverse model aims in a first step to reconstruct the initial pressure map from the measured photoacoustic signals, then, in a second step, to decouple the optical and acoustic parameters from the initial pressure map. Different reconstruction methods are discussed to reconstruct the initial pressure map and different approaches are finally introduced for QPAT.

  by solving the linear Eq. (3.7), we can get the unknown u .

  , n ˆis the normal to the boundary Ω ∂ atξ , r is the position. In our case, instead of solving the RBC, we introduced an extrapolated boundary ext Ω ∂ outside the physical boundary Ω ∂ , inside the physical boundary, Dirichet boundary is applied.

  Set the boundary conditions, for the external boundary, it fits the Dirichlet boundary conditions, for the internal boundary, it fits the Neumann boundary conditions; is equivalent to a perfect absorbing medium surrounding the domain Ω . (M. Schweiger et al., 1995), while The Neumann boundary condition assumes that Ω

Fig. 3 -

 3 Fig.3-2 mesh in the simulation consists of 4296 elements

Fig. 3 -

 3 Fig.3-3 shows the light fluence computed in our case. To highlight the light fluence perturbation introduced by the presence of a tumour, we used two geometries, one geometry without tumour (homogeneous), one geometry with tumour inside (inhomogeneous), the light fluence generated by the former geometry minus that by the latter shows the perturbation fluence.

  Fig.3-3(a) is the logarithm of the light fluence in homogeneous environment (

  the k-space operator, ref c is a scalar reference sound speed.For large scale acoustic simulation, this small correction becomes important. Without the term, the finite difference approximation of the temporal derivative introduces phase errors which accumulate as the simulation runs.The acoustic wave equation can be expressed by the equation of motion, equation of continuity, and equation of state (P. M.[START_REF] Moese | theoretical acoustics[END_REF]

  such that only the components of the wave field travelling within the PML and normal to the boundary are absorbed. And the acoustic equation in Eq. (3.16) becomes: absorption, which is only nonzero inside the PML. Now we refer to the general form of photoacoustic equation, and seek to the solutions to the equation in k-space.

  Fourier transform at all positions and any time step.

  Fig.3-4 mesh projection from FEM to k-space

Fig. 3

 3 Fig.3-6(a) initial pressure in homogeneous environment ( Pa ); (b) initial pressure in inhomogeneous environment ( Pa ); (c) initial pressure perturbation introduced by the tumour ( Pa ).

  Fig.3-7(a) acoustic pressure in homogeneous environment; (b) acoustic pressure in inhomogeneous environment; (c) acoustic pressure perturbation introduced by the tumour.

Fig. 3 -

 3 Fig.3-8 location and distribution of the detectors

Fig. 3 -

 3 Fig.3-10 is the acoustic pressures collected by 18 detectors, the length of each detector is 6mm, each detector is composed of 8 points. The tumour is located as shown in Fig.3-8.

Fig. 3 -

 3 Fig.3-11 shows the photoacoustic pressures from detector of length 6 mm collected by the 5 th detector.

Fig. 3 -

 3 Fig.3-11(a) is the photoacoustic pressure in homogeneous environment; (b) is the photoacoustic pressure in inhomogeneous environment, (c) is the perturbation of the photoacoustic pressure introduced by the tumour.

Fig. 3 -

 3 Fig.3-12 is the acoustic pressures collected by 135 detectors, the length of each detector is 0.8mm, each detector is composed of 8 points.

  Fig.3-12(a) is the photoacoustic pressure in homogeneous environment; (b) is the photoacoustic pressure in inhomogeneous environment, (c) is the perturbation of the photoacoustic pressure introduced by the tumour.

Fig. 3 -

 3 Fig.3-13 shows the photoacoustic pressures from detector of length 0.8 mm collected by the 37 th detector.

  Fig.3-14 the acoustic pressure at inhomogeneous acoustic speed as the tumour size and position varies.

Fig. 3 -

 3 Fig. 3-16 flow of time reverse method

Fig. 3 -

 3 Fig.3-17 (a) initial pressure map, (b) reconstructed initial pressure map.

  Fig.3-18 (a) initial pressure map, (b) reconstructed initial pressure map.

  of the light fluence due to absorption and diffusion perturbation.For both absorption and diffusion perturbations,

  introducing Eq. (3.37), we can rewrite Eq. (3.41) as

,

  matrix. Hence, the global forward problem, for all considered sources [ ] formulated linearly with the following matrix form: H is the Hessian. The multiplication by T J allows handling an inversion process with smaller dimensions of the matrices involved (reduced to the dimension of the chosen reconstruction mesh). Here, the main problem remains in expressing the Jacobian matrix. The medium is first meshed into L voxels of volume

  surface) for the values of the optical properties is chosen to be slightly different from the modeling mesh (

.

  The diffusion equation is solved here numerically with the Finite Element Method (FEM) with Dirichlet boundary conditions at the external frontier of the domain and Neumann boundary conditions elsewhere. The optical parameters are assigned to the different subdomains, point or line sources, with unit intensity and different lengths are distributed around along the above described square. The medium is decomposed into a triangular mesh { } ∆ r (Fig. 3-19) and the equations are solved for different sources geometries in order to compute in equation (3-47). As the meshes have to be changed with the sources positions and geometry, the resulting matrices are then projected into the regular meshes { } j r and { } l largest matrices J and ∆ are computed once and are not involved in the reconstruction process.

-S

  regularization was adopted to avoid instabilities caused by poor conditioning of H: the values of the diagonal of 1 are replaced by

(

  gradient style after smart decomposition (Lanczos process) of H.

Fig. 4 -

 4 Fig. 4-1 distribution of source and detectors, in this geometry, the radius of the circular detectors is 22.5mm, Source1 is located at the centre of the geometry, source2-source5 are located at 4, 8, 12 and 16mm from the centre of the geometry. The amplitude of the each source is unit.

Fig. 4 -

 4 Fig. 4-2 (a) intense detectors Fig.4-2 (b) sparse detectors

Fig. 4 -

 4 Fig. 4-3 (a)-(c) are sparse detectors, in (a), the length of each detector is 0.8 mm, in (b), the length of each detector is 1.6 mm, in (c), the length of each detector is 2.4 mm; (d)-(e) are intense detectors, in (d), there are 135 detectors, the length of each detector is 0.8 mm, in (e), there are 67 detectors, , the length of each detector is 1.6 mm.

Fig. 4 -

 4 Fig. 4-4 effect of the length of the detectors to the reconstructed initial pressure map (a) sparse sensor (b) intense sensor

Fig. 4 -Fig. 4 -

 44 Fig. 4-5 shows some examples of reconstructed initial pressure maps from different detectors distribution.

Fig. 4 -

 4 Fig.4-6 the radius of the breast is 25mm, the tumour located at the centre of the breast, the detectors distributed over the periphery of the breast.

Fig. 4 -Fig. 4 -

 44 Fig. 4-7 shows some examples of reconstructed initial pressure maps through different detectors distributions.

Fig. 4 -

 4 Fig. 4-8 quadratic error map of (a) sparse detectors, (b) intense detectors.

Fig. 4 -Fig. 4 -

 44 Fig. 4-10 shows the value of the quadratic error (1=100%, 0=0%) with different detectors radius and different sizes of tumours. (a)-(c) are the errors caused by sparse detectors, (d)-(f) are the errors caused by intense detectors. The image obtained with intense detectors distribution has better quality than with sparse; for both sparse and intense detectors distribution, as the detector radius increases, the image quality increases.

Fig. 4 -

 4 Fig.4-12 (a) photoacoustic signal generated by a disc source with a radius of 0.1 mm detected from a distance of 10 mm, (b) spectrum of the photoacoustic signal.

Fig. 4 -

 4 Fig.4-13 displays the maximum detectable frequency at different distance with sources of different radius, from this graph, we can see that, as the detection distance increases, the maximum

Fig. 4 -

 4 Fig.4-14 the breast has a radius of 25 mm, we used water to simulate the surrounding environment, this environment is a square with a side length of 60 mm, the tumour located (10 mm, 10mm) away from the centre of the breast, the optical source located 2 mm from the boundary of the breast.

Fig. 4 -Fig. 4 -Fig. 4 -

 444 Fig.4-15 displays the light fluence distributions with 1 source, 2 sources, and 4 sources, respectively. We evaluate the influence of the sources distribution through the corresponding reconstructed initial pressure (Fig.[4][5][6][7][8][9][10][11][12][13][14][15][16]. We can see that, with more sources, we can get a more

  reconstruction methods with different noise level: this preliminary step allowed to select the proper reconstruction algorithm;

  0.3cm×1.1cm and 0.3cm×0.3cm, 0.6 cm from the centre). The sources are located 0.3 cm away from the reconstruction area, both the number and the length of the sources are variable. Point detectors are located 2 cm from this computation area, right at the periphery of the object.

Fig. 4 -

 4 Fig.4-17 Geometry of the simulated object, light sources and transducers are placed respectively 0.3 cm and 2 cm away from the object. In the reconstruction area are placed the perturbations: two absorbers ( a δµ , positions:(0.6cm;1cm)

  10 , reconstructions are expected to be unstable and highly sensitive to noise.The three different reconstruction methods, SVD (with regularization when indicated), LSQR and ART (with non-negative constraint, 1000 iterations, and to 10 -2 ). 4-18 shows the maps of the reconstructed perturbations on the absorption ( a µ δ ) and diffusion ( D δ ) coefficients through different reconstruction methods under increasing noise levels. On Fig. 4-19 are represented cross-plots extracted from these maps along a horizontal line through the middle of the reconstruction area. As expected, reconstructions are highly sensitive to noise, a high value of the regularization parameter has to be introduced in the SVD, showing only few elements of H out of noise may be useful for the reconstruction, even at low noise level. With

Fig. 4 -

 4 Fig. 4-18: Comparison of different reconstruction methods at different noise levels: results of the reconstructions of the perturbations on the absorption a µ δ and diffusion D δ coefficient maps.

Fig. 4 -

 4 Fig. 4-20 : Left: Normalized singular values of H; right: condition number of H as a function of the number of sources.

Fig. 4 -

 4 Fig. 4-21 : Influence of the number of sources: reconstructed perturbation absorption a µ δ and diffusion D δ coefficient

Fig. 4 -Fig. 4 -

 44 Fig. 4-23 : Left: Normalized singular values of H; right: condition number of H for a) even and sparse distribution around the object; b) 4-sides examination but short spacing between detectors; c) same situation but with even shorter (compressed) spacing; d) 2-sides examination horizontal detection; e) 2-sides examination vertical detection.

Fig. 4 -

 4 Fig. 4-25 : Influence of the detectors distribution: cross-plots extracted from Fig. 4-24.of reconstructed perturbation absorption a µ δ and diffusion D δ coefficient maps. (a) the detectors are evenly distributed around the object (15 detectors on each side); (b) on each side, the detectors cover a length that is half the length of the object; (c) the length covered by the detectors is half of (b); (d) 30 detectors evenly distributed on each of the two vertical sides of the object;(e) 30 detectors evenly distributed on each of the two horizontal sides of the object.

Fig. 4 -Fig. 4 -

 44 Fig. 4-26 : Left: Normalized singular values of H; right: condition number of H as a function of the number of sources.

Fig. 4 -

 4 Fig. 4-28 : Influence of the sources shape: cross-plots extracted from Fig. 4-27 of reconstructed perturbation absorption a µ δ and diffusion D δ coefficient maps with with 4 points sources; 4 line sources with a length of 0.2 mm; 4 mm; 8 mm; 26 mm and full field of view probing with one single source composed of 4 lines of 26 mm length.

Fig. 4 -Fig. 4 -

 44 Fig. 4-29 : Experimental situations: (a) 4 points sources, 15 detectors; (b) 4 line sources with length 26 mm, 15 detectors; (c) same line sources, 1 detector; (d) same line sources, 2 detectors; (e) same line sources, 3 detectors. The condition numbers of H are reported for each situation.

Chapter 5 .

 5 Modeling the experimental situations ...............................................................

Fig. 5 -

 5 Fig. 5-1 is the schema of our experimental system. It illustrated the main component of the system. Fig. 5-2 illustrates each part in detail. The position of the optical source, the acoustic detector and the phantom. In the experimental setup, there is just one source and detector pair, but the phantom can rotate. The development of the setup has been conducted in between LMA and IF, and partly funded by CNANO PACH (AAP 2010) and ANR AVENTURES (ANR-12-BS01-0001-04).

Fig. 5

 5 Fig.5-1: Schema of the system

  Fig.5-4 schema of breast phantom, acoustic detector and optical source

Fig. 5 -

 5 Fig.5-5 breast phantom

Fig. 5

 5 Fig.5-6 simulation geometry considering the experimental setup

Fig. 5 -

 5 Fig. 5-7 is the light fluence generated in the simulation.

  detector

Fig. 5 -

 5 Fig.5-7 light fluence generated in the simulation

Fig. 5 -

 5 Fig.5-8 comparison of photoacoustic pressure generated in the simulation and in the experiment

Fig. 5 -

 5 Fig.5-10 experimental data after filter

Fig. 5 -

 5 Fig.5-11 are the simulated and experimental pressures when the transducer is 1cm, 2 cm, and 3 cm away from the boundary of the phantom, respectively.

Fig. 5 -

 5 Fig.5-11 experimental and simulated pressure detected by transducers 1cm, 2 cm, and 3 cm away from the phantom

Fig. 5 -

 5 Fig.5-14 2D illustration of the rotated source-detector pair

Fig. 5 -

 5 Fig.5-15 2D illustration of the breast geometry

  Fig.5-16 showed the light fluence distribution from these four arc sources.

Fig. 5 -

 5 Fig.5-16 light fluence distribution from 4 arc sources

Fig. 5 -

 5 Fig.5-19 arc source and line detector setup

Fig. 5 -Fig. 5 -

 55 Fig.5-22 arc source and focussed detector setup

FigFig. 5 -

 5 Fig.5-28 semicircle source and focussed detector setup

  the present work is the development of non invasive methods for breast cancer diagnosis and imaging. Through a review of conventional imaging techniques, PhotoAcoustic Tomography (PAT), a new emerging imaging technique, seems to be a promising non invasive alternative to conventional imaging techniques for breast cancer diagnosis. It may combine the advantages of two other well-known non-invasive techniques, that is Diffuse Optical Tomography (DOT) on the one hand and UltraSound Tomography (UST) on the other hand, and it potentially palliates their drawbacks.

Table . 1

 . -1 the countries with the top 20 highest incidence of breast cancer in 2012

	2	Denmark	105.0
	3	France (metropolitan)	104.5
	4	The Netherlands	99.0
	5	Bahamas	98.9
	6	Iceland	96.3
	7	United Kingdom	95.0
	8	Barbados	94.7
	9	United States of America	92.9
	10	Ireland	92.3
	11	French Polynesia	92.2
	12	Germany	91.6
	13	Italy	91.3
	14	Finland	89.4
	15	Luxembourg	89.1
	16	New Caledonia	87.6
	17	Australia	86.0
	18	Malta	85.9
	19	New Zealand	85.0
	20	Switzerland	83.1
	Rank	Country	Age-standardized rate per 100,000
	1	Belgium	111.9

Source: Ferlay J, Soerjomataram I, Ervik M, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray, F. GLOBOCAN 2012 v1.0, Cancer Incidence and Mortality Worldwide: IARC CancerBase No. 11 [Internet]. Lyon, France: International Agency for Research on Cancer; 2013. Available from: http://globocan.iarc.fr, accessed on 13/12/2013.

31 2.1.1 Physical properties of biological tissues .............................................................. 32 2.1.1.1Optical contrasts in biological tissues ................................................................ 32 2.1.1.2Photometric quantities ........................................................................................ 35 2.1.1.3The Radiative Transport Equation ................................................................... 35 2.1.1.4The Diffusion Equation ....................................................................................... 38 2

  Forward and Inverse Models in PhotoAcoustic Tomography ................................. 31 2.1 Introduction .................................................................................................................... .2 Photoacoustic phenomenon ............................................................................................ 39 2.3 Resolution methods ........................................................................................................

	Chapter 2. Forward and Inverse Models in PhotoAcoustic Tomography
	Chapter 2
	chapter starts by a presentation of data facts related to breast cancer nowadays,
	illustrates the importance of breast cancer early diagnosis. Various of conventional techniques for
	breast cancer diagnosis are discussed, among these, photoacoustic imaging is gaining increasing
	importance. An overview of the state of the art in photoacoustic imaging is given in this chapter, the
	application of photoacoustic imaging is discussed in this chapter. In the last part, two existing
	photoacoustic imaging systems for breast cancer are introduced.

41 2.3.1 Resolution of the Diffusion Equation .................................................................. 41 2.3.2 Acoustic resolution ................................................................................................ 44 2.3.3 photoacoustic resolution ....................................................................................... 45 2.3.4 Inverse problems ................................................................................................... 48 2

  .4 Conclusion ......................................................................................................................

53 2.1. Introduction

	pressure distribution. Then the forward problem resolution is introduced: it consists in calculating
	the photoacoustic pressures measurements, knowing the geometrical and physical properties of the
	medium. The inverse problem is then formulated. It can be handled with multiple approaches
	according to the information one seeks to reconstruct. The simple inversion consists in the
	reconstruction of the initial pressure map. The problem is a purely acoustic problem of acoustic
	source localization. However, intuitively, one can imagine the information content of these signals is
	much more important: the signal, measured acoustically, vehicles information on both the optical
	and acoustic properties of the medium. The main problem consists in finding a way to uncouple the
	contribution of each of these parameters. It is shown that this can be partly achieved with multiple
	sources illumination. The different approaches to formulate this inverse problem are presented.
	Most biological tissues are referred to as turbid media in terms of light propagation. They are
	characterized by strong optical scattering and weak optical absorption, at least in the 400-1350 nm
	spectral region. The mean free path between photon scattering events is on the order of 0.1 mm,
	whereas the mean absorption length (mean path length before photon absorption) can extend to 10-
	100 mm. There are several models to describe photons migration in tissues. The mostly common
	used one is the Radiative Transport Equation (RTE). According to the thickness of the tissue as
	compared with the transport mean free path, this equation can be reduced to a diffusion equation
	when the Diffusion Approximation (DA) holds. Once the tissue illuminated, an initial photoacoustic
	pressure distribution is generated due to the absorption and conversion of radiative energy into heat
	and the thermoelastic expansion. The photoacoustic pressure will then propagate into the tissue
	according to a Helmholtz propagation equation.
	In this chapter are first defined the main physical contrasts producing the PA phenomenon
	and the governing wave propagation equations. The two mathematical models are described in
	details: the RTE with its diffusion approximation and the photoacoustic wave equation with the initial

Table .

 .

				3-1 optical parameters
	absorption coefficient	Reduced scattering coefficient
	μ a_water =0			μ s_water =1e -3 m	-1
	μ a_breast =4m	-1	μ s_breast =1000m	-1
	μ a_tumour =8m	-1	μ s_tumour =1000m	-1

Table . 3

 . -2 acoustic parameters

Table .

 .

		3-3 reference of the sound speed in breast
	tissue	sound speed (m s -1 )	reference
	water	1403-1555	(Martin Greenspan and Carroll E. Tschiegg, 1957)
	breast	1428-1570

Table 5

 5

					. 1 simulation parameters
	Absorption coefficient	Scattering coefficient
	μ a_water =10m	-1		μ s_water =1/3000m	-1
	μ a_breast =50m	-1		μ s_breast =100m	-1
	μ a_tumour =25000m	-1	μ s_tumour =100m	-1

Table 5

 5

	. 2 optical parameters in the simulation
	absorption coefficient	Scattering coefficient
	μ a_water =0	μ s_water =1e-5cm	-1
	μ a_breast =0.1cm

tumour to detector. The former indicates the distance from tumour to detector; the latter indicates the distance from source to tumour plus the distance from tumour to detector . This second signal can be used to localize the position of the tumour. The conventional photoacoustic tomography is to reconstruct the initial pressure map from collected acoustic pressures. According to (Y. [START_REF] Hristova | reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media[END_REF] and [START_REF] Cox | Artifact trapping during time reversal photoacoustic imaging for acoustically heterogeneous media[END_REF], time reversal reconstruction method is the "least restrictive" image reconstruction method: it works for arbitrary geometry of the closed observation surface, the source doesn't need to be inside the surface, it works with a heterogeneous sound speed distribution. So in our work, we used time reversal to do the reconstruction.

The principle of time reversal lies in that: the solution of the acoustic pressure ) , ( t p r to the initial value problem is identical to the solution for sufficiently large T . In other words, for an initial pressure with a bounded area S , there exists a moment T that the acoustic pressure ) , ( t p r vanishes inside S for any T t > , then one can Impose zero initial conditions at T t = and boundary conditions in the reversed time direction. This assumption is true when the sound speed is constant and the spatial dimension is odd, the Huygens' Principle holds. When the dimension is even, or when ). The FEM mesh typically contains 4300 elements, the modeling mesh is a regular grid of 61×61 pixels (pixel size: 0.1 cm), and the reconstruction mesh is 21×21 pixels, 1 pixel shifted in horizontal and vertical with the modeling mesh.

)) ( ( 0 , j s l r u β is then obtained by assembling the sensitivity matrices, and its dimension is This result shows the possibility of using simplified experimental PAT setups. Reducing the number of detectors points, up to 3 only here, corresponding to 12 measurements only, of course worsens the quality of the reconstructions but still allows getting better results than with complete angle detection: in comparison with Fig. 4-27 (e), obtained with 240 measurements, the number of measurements is reduced by a factor 20.

Synthesis

In this study, a multiple illumination photoacoustic tomography algorithm has been implemented according to the method suggested in [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF]. The method is quantitative in the sense that it allows reconstructing both the absorption and diffusion coefficient while the conventional approaches in PAT, with a single wide field illumination give the access to the initial pressure distribution map. A description with explicit practical implementation details was provided.

Different linear reconstruction algorithms were tested: the ART with non-negative constraints was shown to be less sensitive to noise than SVD or LSQR. Even if not involved in these reconstruction loops, the computation of the Jacobian matrix is still a burden, especially if one seeks to extend the approach to the non-linear case. Gradient-based methods may be used instead because more memory-efficient as they require only gradient information to approximate the Hessian matrix (G. [START_REF] Bal | On multi-spectral quantitative photoacoustic tomography in diffusive regime[END_REF][START_REF] Gao | Bregman methods in quantitative photoacoustic tomography[END_REF][START_REF] Saratoon | A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation[END_REF].

With this tool, a 2D simulation study on a variety of sources and detectors geometries situations was conducted. The sets of simulations were conducted through the examination of the conditioning of the Hessian matrix. The reconstructions highlighted several important results:

-Optical parameters are reconstructed with improved quality by using wide field illuminations: this study demonstrates that it is preferable to use enlarged sources instead of multiplying the number of point sources illuminations.

-The position of the detectors with respect to the perturbations to be reconstructed is crucial and deserves additional prospection. When comparing the influence of sources and detectors, the role of the illumination and detection schemes are nested and complex. The considered, with dimensions and physical properties corresponding to the breast and tumours (taken from literature).

A first study was conducted on the optimization of the distribution of the detectors. Synthetic data were generated with focused detectors distributed under dense or sparse distributions with varying number and/or length. The simulations were conducted under 2D and 3D configurations. By using the time-reversal method, the initial pressure distribution maps were obtained under the different configurations. These reconstructions showed, in 2D and 3D:

-Dense detectors distributions give higher quality reconstructions, as expected, especially if the size of the detectors is small; -Sparse detector distributions can also produce good quality reconstructions provided the angular coverage around the object is complete.

This study is incomplete but the proposed approach can be followed in order orientate the choice of the detectors and the protocol of data acquisition as a function of the size of the object to be reconstructed. These results have been presented through an oral communication at the 1 st EOS Topical Meeting on Frontiers in Optical Imaging (FOI 2013).

In a second study, through a perturbation approach (Born approximation), a linear reconstruction scheme can be established for the reconstruction of the absorption and diffusion coefficients maps [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF][START_REF] Zemp | Quantitative photoacoustic tomography with multiple opticalal sources[END_REF]. The reconstruction scheme was implemented and its robustness to noise was tested. The singular value decomposition of the Hessian shows a very bad conditioning: the reconstruction scheme is highly sensitive to noise. By playing with the illumination geometry, we sought for a better conditioning of the matrices. It is shown that:

-Two illumination sources are not sufficient to decouple absorption and diffusion -The method is really improved if wide field multiple illuminations are used instead of point sources.

This study has been reported in a paper submitted to Biomedical Optics Express.

Finally the simulations were compared to measurements performed with the PAT system developed at the Institut Fresnel. The acquisition geometry is highly constrained and simulations with limited views, under transmission mode with rotating cylindrical samples, were conducted in order to test if QPAT can be envisioned with this setup. It is shown that uncoupled absorption and diffusion coefficients maps can be obtained with only 3 illumination sources per angle. However, only few measurements were acquired with the setup that is still under development.

The theoretical tools developed within this thesis project have to be more systematically tested against experimental data but prior to that, they have to be improved in order to incorporate signal pre-and post-processing. However, they have already shown to be useful in defining new and Dans la gamme de longueurs d'onde visible, les tissus biologiques sont caractérisés par une forte diffusion et absorption optique, au moins dans la région spectrale de 400-1400nm. Le libre parcours moyen de diffusion est de l'ordre de 0,1 mm, tandis que le libre parcours moyen d'absorption peut s'étendre de 10 à 100 mm dans la gamme 600-900 nm, dites « fenêtre thérapeutique ». Il existe plusieurs modèles pour décrire la diffusion lumineuse dans les tissus. Le plus utilisé est l'Equation de Transfer Radiatif (ETR). Cette équation peut être réduite à une équation de diffusion lors que les hypothèses de l'Approximation de la Diffusion (AD) sont satisfaites. Dans la « fenêtre thérapeutique », le sein est un milieu diffusant qui satisfait ces approximations : il s'agit d'un milieu épais, dont les coefficients d'absorptions sont bien plus faibles que les coefficients de diffusion. La propagation lumineuse dans les tissus s'accompagne d'un échauffement lié à l'absorption et à la conversion de l'énergie lumineuse en chaleur. Ceci génère une dilatation thermoélastique qui produit une distribution de pression photoacoustique initiale. Celle-ci se propage alors dans le tissue selon une équation de propagation de Helmholtz. Le problème inverse peut être traité avec des approches multiples en fonction des informations que l'on cherche à reconstruire. Le plus simple consiste en une inversion de la reconstruction de la carte de pression initiale. Le problème est un problème purement acoustique de localisation de sources acoustiques.

Toutefois, intuitivement, on peut imaginer le contenu informationnel de ces signaux est beaucoup plus riche: le signal mesuré acoustiquement véhicule des informations sur les propriétés optiques et acoustiques du milieu. Le principal problème consiste à trouver un moyen de découpler la contribution de chacun de ces paramètres. Il est démontré que cela peut être partiellement réalisé avec plusieurs sources d'éclairage.

Un code de calcul multiphysique a tout d'abord été mis en oeuvre. Il permet de résoudre le problème direct, en incorporant le couplage des modèles de propagation lumineuse de la lumière et des ultrasons. Les vitesses de propagation des deux phénomènes physiques étant très différents, différentes méthodes de résolution ont été adoptées pour résoudre les équations aux dérivées partielles (EDP). Il en résulte un code de calcul hybride: la méthode des éléments finis (MEF) classique a été choisie pour résoudre l'équation de propagation optique, tandis qu'une méthode semi-analytique basée sur des transformées de Fourier (méthode K-Space) a été préférée pour résoudre l'équation de propagation acoustique. Dans une approche directe, non quantitative, des techniques classiques, héritées du palpage échographique, ont été adoptées: à partir des mesures de pression acoustique, le problème peut être considéré comme un problème de localisation de sources acoustiques, et la méthode du Retournement Temporel (MRT) peuvent être utilisée pour reconstruire des cartes de distribution de pression initiale. Nous appelons cette méthode de sondage «passive», elle est analogue au problème de la détection des tremblements de terre en sismologie, par exemple. Ce solveur direct permet de générer des données synthétiques qui peuvent être éventuellement polluées par l'introduction de bruit.

Les simulations ont été réalisées en tenant compte de situations expérimentales lorsque cela est possible. A cet effet, les détecteurs modélisés sont focalisés, les fantômes numériques sont des disques (2D) ou cylindriques (3D), avec des dimensions et des propriétés physiques correspondant à celles du sein et des tumeurs (valeurs extraites de la littérature). Dans une deuxième étude, grâce à une approche perturbative (approximation de Born), le problème direct, dont les variables sont les paramètres physiques, peut se formuler de manière linéaire, du moins en ce qui concerne les paramètres optiques, à savoir, les coefficients d'absorption et de diffusion. Des algorithmes d'inversion linéaire conventionnels peuvent être utilisés pour la reconstruction de cartes découplées de coefficients d'absorption et de diffusion [START_REF] Shao | Iterative algorithm for multiple illumination photoacoustic tomography (MIPAT) using ultrasound channel data[END_REF][START_REF] Zemp | Quantitative photoacoustic tomography with multiple opticalal sources[END_REF]. Ce schéma de reconstruction a été mis en oeuvre et sa robustesse au bruit a été testée. La décomposition en valeurs singulières de la matrice Hessienne montre un très mauvais conditionnement: l'algorithme de reconstruction est en effet très sensible au bruit. En jouant sur la géométrie d'éclairement, nous avons cherché à obtenir un meilleur conditionnement des matrices.

Dans ce travail, il est démontré que:

• Deux sources éclairages ne sont pas suffisantes pour découpler l'absorption et la diffusion;

• La méthode est vraiment améliorée si l'on utilise des sources d'illumination larges et non plus ponctuelles.

Cette étude a été rapportée dans un document soumis à Biomedical Optics Express.

Abstract

The present work was motivated by the development of alternative imaging techniques for breast cancer early diagnosis, that is photoacoustic imaging, which potentially couples the merits of optical imaging and ultrasound imaging, that is high optical functional contrasts brought by optical probing and high spatial resolution by ultrasound detection.

Our work aims at modeling the photoacoustic multiwave phenomenon and incorporate it in an efficient reconstruction algorithm to solve the inverse problem. The inverse problem consists in the recovery of interior maps of physical properties of the breast. The forward model couples optical and acoustic propagations. The Finite Element Method (FEM) was chosen for solving the optical propagation equation, while a semi-analytical method based on Fourier transforms calculations (kspace method)was preferred for solving the acoustic propagation equation. For the inverse model, time reversal method was adopted to reconstruct the initial pressure distribution, an active approach of the inverse problem was also achieved, which decoupled the optical properties from measured photoacoustic pressure, this approach is called quantitative photoacoustic tomography (QPAT), in this approach, illumination/detection protocol was studied, and the experimental set up is also take into consideration.

In the last step, photoacoustic pressure measurements obtained from experiment and simulation are studied and compared.

Key words: breast cancer, photoacoustic imaging, quantitative tomography