Ecole Sigle

Doctorale

Renée El

Melhem Bât

Pascal Blaise

R Gourdon

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

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

M Philippe
email: philippe.normand@univ-lyon1.fr

Mme Emmanuelle
email: emmanuelle.canet@univ-lyon1.fr

Étage

M Luca

Jean-Yves Buffière
email: jean-yves.buffiere@insa-lyon.fr

M Jocelyn Bonjour

Viviane Polsinelli
email: viviane.polsinelli@univ-lyon2.fr

Brigitte Dubois

J Y Toussaint

M Christian Montes

Ricca Alessandri

Patrick Un

Clémence, Tanya, Rémi, Mizpah Sarah Melvil

Louise Tom

Nina Johnatan

Maxime D Et Valentin

Signal-to-noise ratio threshold for VCA algorithm

Keywords: è possibile valutare la crescita, l', come primo passo, la discriminazione dei dati

Photoacoustic imaging is a functional technique based on the creation of acoustic waves from tissues excited by an optical source (laser pulses). The illumination of a region of interest, with a range of optical wavelengths, allows the discrimination of the imaged media. This modality is promising for various medical applications in which growth, aging and evolution of tissue vascularization have to be studied. Thereby, photoacoustic imaging provides access to blood oxygenation in biological tissues and also allows the discrimination of benign or malignant tumors and the dating of tissue death (necrosis).

The present thesis aims at developing a multispectral photoacoustic image processing chain for the calculation of blood oxygenation in biological tissues. The main steps are, rst, the data discrimination (clustering), to extract the regions of interest, and second, the quantication of the dierent media in these regions (unmixing).

Several unsupervised clustering and unmixing methods have been developed and their performance compared on experimental multispectral photoacoustic data. They were acquired on the experimental photoacoustic platform of the laboratory, during collaborations with other laboratories and also on a commercial system. For the validation of the developed methods, many phantoms containing dierent optical absorbers have been produced. During the co-supervision stay in Italy, specic imaging modes for 2D and 3D real-time photoacoustic imaging were developed on a research scanner. Finally, in vivo acquisitions using a commercial system were conducted on animal model (mouse) to validate these developments.

Aneline DOLET i Résumé L'imagerie photoacoustique est une modalité d'imagerie fonctionnelle basée sur la génération d'ondes acoustiques par des tissus soumis à une illumination optique (impulsion laser). L'utilisation de diérentes longueurs d'ondes optiques permet la discrimination des milieux imagés. Cette modalité est prometteuse pour de nombreuses applications médicales liées, par exemple, à la croissance, au vieillissement et à l'évolution de la vascularisation des tissus. En eet, l'accès à l'oxygénation du sang dans les tissus est rendu possible par l'imagerie photoacoustique. Cela permet, entre autres applications, la discrimination de tumeurs bénignes ou malignes et la datation de la mort tissulaire (nécrose). Ce travail de thèse a pour objectif principal la construction d'une chaîne de traitement des données photoacoustiques multispectrales pour le calcul de l'oxygénation du sang dans les tissus. Les principales étapes sont, d'une part, la discrimination des données (clustering), pour extraire les zones d'intérêt, et d'autre part, la quantication des diérents constituants présents dans celles-ci (unmixing). Plusieurs méthodes non supervisées de discrimination et de quantication ont été développées et leurs performances comparées sur des données photoacoustiques multispectrales expérimentales. Cellesci ont été acquises sur la plateforme photoacoustique du laboratoire, lors de collaborations avec d'autres laboratoires et également sur un système commercial. Pour la validation des méthodes développées, de nombreux fantômes contenant diérents absorbeurs optiques ont été conçus. Lors du séjour de cotutelle de thèse en Italie, des modes d'imagerie spéciques pour l'imagerie photoacoustique 2D et 3D temps-réel ont été développés sur un échographe de recherche. Enn, des acquisitions in vivo sur modèle animal (souris) au moyen d'un système commercial ont été réalisées pour valider ces développements.

Spectra of two dierent pixels: in green, pixel from the background with low PA amplitudes at all wavelengths and in purple, pixel from the ROI with signicant PA amplitudes. . . .

3.2

The spatio-spectral mean-shift principle. (a) Pixel spatial features at the rst iteration with the R S parameter, and (b) spectra of three pixels with the R λ parameter. The red pixel is the reference pixel. Although the green and blue pixels are both in the spatial area of the reference pixel, only the green pixel is used to update the features of the reference.

Indeed, the spectrum of the blue pixel is not close enough to the one of the reference sample (by at least one wavelength; the blue and red spectra are separated by a distance greater than R λ ). . Weighting factor applied to the signal acquired by the q-th transducer p q q-th transducer's RF signal s q q-th transducer position τ q (s) The PhD manuscript is organized as follows. The photoacoustic imaging principle is initially reviewed in Chapter 1. The modality, the multispectral aspect and its interests are introduced. A state of the art is focused on systems and algorithms in the context of the thesis. Chapter 2 presents the rst contribution on experimental developments. The dierent used systems, as well as the imaged samples, are described. Real-time imaging developments are also presented and the dataset used to validate the developed methods are introduced.

Travel
The proposed processing methods can be separated in two dierent parts which have dierent objectives. At the beginning, the regions of interest (ROI) need to be identied since, for each application, only some biological tissues are of interest. This rst step is called clustering step. Chapter 3 describes the clustering methods and highlights their performances. The second processing part is presented in Chapter 4 and applied only to the identied ROI. It aims at quantifying the data contained in these regions (unmixing step). The performances of several quantication methods are compared.

Finally, Chapter 5 highlights the application of the proposed methods to in vivo data for the calculation of blood oxygen concentrations. It also presents preliminary results on the evaluation of tissue death. The manuscript ends with a general conclusion, including a discussion about the limitations of the developed methods.

Aneline DOLET

Introduction

The photoacoustic eect, also called optoacoustic or thermoacoustic eect, has been discovered in 1880 by Alexander Graham Bell who has created sound by periodically illuminating an object with sun light [Bell 1880]. However, the sun light is not a satisfying optical source to allow technological developments of this eect: not stable and hard to make pulsed illuminations. So, almost no developments were made of this discovery during nearly a century. The development of stable pulsed optical excitations (pulsed lasers) in the 70 th started a new interest of the photoacoustic eect with various applications. First, it has been used on gaz [Kreuzer 1971] and, some years later, the use on solid media has been validated [START_REF] Rosencwaig | [END_REF]]. It is during the 90 th that techniques using the photoacoustic eect took their rise with an increasing interest for medical applications [Li 2009]. Since then, the number of medical photoacoustic imaging applications constantly grows up [START_REF] Esenaliev | [END_REF][START_REF] Beard | [END_REF], Castelino 2008, Arabul 2015].

Aneline DOLET 3 Chapter 1. Photoacoustic imaging

In this Chapter, the concept of photoacoustic imaging is explained. Its multispectral uses are introduced and the main used processing methods and systems are presented. Some applications of multispectral photoacoustic imaging are also highlighted. The context and the motivations of the present thesis are then described as well as the organization of this manuscript.

1.2 Combination of optical and ultrasound imaging 1.2.1 Optical imaging

When biological tissues are excited by an optical source, the light propagates inside the media. During the propagation, the light is submitted to various phenomena as scattering and absorption which means that the light is attenuated. When an optical absorber is excited by a photon, its absorption leads to the emission of a lower energy photon (uorescence phenomenon) or an increase of temperature which could permit the access to tissue functional information. These information are of particular interest for various applications.

Lots of optical imaging techniques exist and are used depending on the targeted application. The general idea is that the resolution decreases when the used optical source allows larger imaging depth.

For example, optical coherence tomography (OCT) has a resolution of some micrometers (≈ 5µm) and an imaging depth of a few millimeters while optical microscopy (OM) can only image a few hundred of micrometers with a resolution < 3µm (see Figure 1.1). announced by the ANSI (American National Standards Institute) and the AFNOR (Association Française de Normalisation). A maximum permissible exposure (MPE) for eyes and skin, depending on the used wavelength, has then to be respected [Vallet 2015]. It means that the optical energy cannot be increased to enlarge the imaging depth.

4

Aneline DOLET

Combination of optical and ultrasound imaging

Finally, even if the optical imaging investigation depth is low, its major interest is the access to functional information with satisfactory contrast. However, the spatial resolution is around or superior to 1mm for an imaging depth up to few centimeters, according to various applications (blood vascularization or tumors imaging). Thus, limited imaging depth and poor resolution do not allow its use for such applications that however require knowing tissue functional information.

Ultrasound imaging

In ultrasound imaging, the probe transmits ultrasound waves that propagate inside the tissue. Biological tissues are composed of numerous inhomogenities that are localized at random positions. When inhomogenities are submitted to an incident wave, they become ponctual sources that produce spherical waves in all directions (Figure 1.2). This corresponds to US scattering. All the inhomogenities of the imaged region become scatterers when submitted to US waves. These scatterers produce then US waves that are back-propagated to the US probe which converts them into radiofrequency signal (RF signal).

Received signals are reconstructed to image the tissue using dierent reconstruction methods.

τ op d 2 s 4D th (1.1)
with d s the specic dimension of the considered structure (in meter) and D th the thermal diusivity (in m 2 .s -1 ). Another constraint called constraint connement has to be respected to make the optical absorber volumetric expansion negligible:

τ op d s c (1.2)
with c the sound speed in biological tissues. The photoacoustic eect is optimally produced when both constraints are satised. Indeed, the optical absorbers have to come back to their initial states between each optical illumination to produce the US waves.

6

Aneline DOLET 1.3. Ultrasound wave generation with pulsed light illumination

Let us consider a biological structure of approximately 10µm (e.g., red blood cell). In biological tissue, D th is around 0.14mm 2 .s -1 [Duck 1990] and c is usually considered around 1500m.s -1 . According to equations (1.1) and (1.2), τ op has to be lower than 1.79e -4 s and 6.67e -9 s, respectively. Pulsed laser of τ op around 5ns -10ns are so considered for photoacoustic imaging.

The third step of ultrasound wave generation (Figure 1.3(c)) can be expressed using the linear equations of uid mechanics to rely the temperature variations to the pressure changes in biological tissue under the previous constraints [Morse 1987, Vallet 2015]. The initial pressure p 0 (s) created by an optical absorber at the spatial position s submitted to optical pulse illuminations is expressed by [Beard 2011]:

p 0 (s) = Γµ a (s)φ(s) (1.3)
where Γ is the Grüneisen parameter (dimensionless) which expresses the proportion of optical energy converted in pressure and depends on the temperature, µ a is the absorption coecient of the considered optical absorber (in cm -1 ) which depends on the used optical wavelength, and φ the optical uence, i.e., the optical energy per unit area (in J.cm -2 ).

Considering soft tissue in ambient temperature, the Grüneisen parameter is Γ ≈ 0.25 [Li 2009]. The initial pressure p 0 generated by a blood vessel excited at 800nm with a uence φ = 20mJ.cm -2 can be calculated knowing its optical absorption µ a ≈ 4.3cm -1 . In this case, the generated initial pressure is ≈ 21.5kP a which corresponds to a temperature increase of ≈ 0.1K.

In photoacoustic imaging, the optical excitation is designed to respect the thermal connement and constraint connement. The acoustic pressure generated is generally low compared to peak pressure of medical ultrasound imaging. Considering the ultrasound attenuation in biological tissue [Szabo 1978], the acoustic pressure acquired by the ultrasound probe is even lower. For this reason, the acquisition geometry needs to be optimized in order to reconstruct the photoacoustic image with the best signal-tonoise ratio (SNR). W q p q (τ q (s))

Photoacoustic image reconstruction

τ q (s) = s-sq c (1.4)
with D the number of detected signals, i.e., the number of active transducers of the probe, W q the weighting factor applied to the signal acquired by the q-th transducer, p q the q-th transducer's RF signal and τ q (s) the travel time from the position s in the medium to the position s q of the q-th transducer.

Post-processing classically used in ultrasound imaging, like time gain compensation or interpolation, can be applied to the reconstructed photoacoustic image I.

Aneline DOLET Chapter 1. Photoacoustic imaging

The ultrasound speed in biological tissues has to be known to reconstruct the images. In soft tissues, the ultrasound speed is close to the one of ultrasound in water. The typically used value is in the range [1480m.s -1 ; 1540m.s -1 ]. As already mentioned, in photoacoustic the temperature increase is around 0.1K, which corresponds to a few kP a pressure wave and leads to a low pressure amplitude. Pulsed lasers of around 5ns pulse duration are also commonly used. Often, a ≈ 10Hz repetition frequency optical source is considered.

Medical ultrasound systems (probes and scanners) are used for photoacoustic imaging, with similar image reconstruction methods [Wells 1977, Angelsen 2000] as the DAS algorithm. The optical absorbers of biological tissues produce acoustic waves between 1M Hz and 100M Hz depending on the size of the absorbers. The choice of the ultrasound probe used to acquire photoacoustic signal is then really important and needs to be adapted to the targeted applications.

The photoacoustic signal is linked to the optical absorption of the imaged media. Depending on the concentration of a single medium, the photoacoustic signal diers because the optical absorption varies.

For example, the dilution of a medium with a ratio of 0.5 involves an optical absorption of 0.5 compared to the one of the pure medium. The concentration of a media impacts the photoacoustic signal which allows then quantication. Finally, photoacoustic imaging is a non-ionizing and non-invasive imaging technique to access biological tissue functional information with respect of light amplitude excitation.

Interest of multispectral illumination in photoacoustic imaging

Using a single optical wavelength already allows the dierentiation between two dierent media since, e.g., oxygenated (HbO 2 ) and deoxygenated (Hb) blood have dierent optical absorptions at 650nm (see As the absorption of a single medium is dierent at each wavelength, the evolution shape of the medium absorption over wavelengths is specic to a single medium. In this way, the use of several wavelengths allows the dierentiation of more than two media and increase the medium discrimination accuracy. The spectral evolution of a medium is then considered for the dierentiation or quantication of media using multispectral photoacoustic data [START_REF] Van Veen | [END_REF], Beard 2011].

In multispectral PA imaging, the wavelength range has to be well chosen; e.g., HbO Finally, the multispectral photoacoustic evolution of a single medium does not have the exact same shape than the multispectral optical absorption evolution of this medium. The multispectral PA signal is indeed linked to the used laser energy variations over wavelengths. It means that the photoacoustic signal of a single medium can be dierent depending on the light source used, then on the photoacoustic system used to acquire the data. Developing processing methods to discriminate or quantify media using photoacoustic imaging data acquired with dierent systems is thus challenging.

Aneline DOLET

Current researches in multispectral photoacoustic imaging

Multispectral PA imaging is a powerful modality to access functional information. Current researches address the following topics: (1) PA image reconstruction [Treeby 2013[START_REF] Arabul | [END_REF], Ding 2017] taking into account scattering and absorption as a function of wavelength. This is not addressed in this thesis as a classical reconstruction algorithm is used (DAS) without considering these features. (2) PA acquisition systems and applications. (3) Data processing methods for medium discrimination and quantication.

Photoacoustic imaging systems and probes

Photoacoustic imaging can be done using classical optical sources coupled with ultrasound probes and systems. Depending on the size of the desired imaged object, dierent groups of systems exist to image structures of size from some micrometers to some centimeters. System types are characterized by their resolution and imaging depth (Figure 1.7), increasing from photoacoustic microscopy (PAM ≈ 100µm) and photoacoustic tomography (PAT ≈ mm), to photoacoustic computed tomography (PACT ≈ cm).

These are the three main system groups [START_REF] Wang | [END_REF]]. Considering the possible imaged region size:

PAT have the larger imaging size with some square centimeters, PACT image region of some square millimeters and PAM of some square micrometers.

PAT acquires the ultrasound waves all around the ROI [Cai 2011]. The acquisition principle is not extended here as an experimental PAT system is deeply presented later in Chapter 2. Briey, to detect the ultrasound waves, the US probe (or a mono-element transducer) rotates all around the imaged sample.

Most of the time, it is the imaged sample that rotates inside a xed laser/ultrasound probe geometry (Figure 1.8(a) [Liu 2016]). These particular systems cannot be used for every application. It is possible for sample like phantoms or small animals but not always for human body. Researches are then done to develop dierent probes to extend the PAT uses, like a concave arc-shape array used by [START_REF] Brecht | [END_REF]], which can be seen PACT is a larger category as it clusters the system composed of laser source classically used in optical imaging and US probe usually used in US imaging. No particular acquisition geometry is dened, the laser and US probe can be put all around the imaged sample in every planes or axes (two dierent geometries are presented in the next Chapter). The only thing is that the geometry does not change during the acquisitions by opposition to PAT systems which rotates. The US probe can be separated from some optical sources [Montilla 2013] (Figure 1.10(a)). The optical sources and the US transducers can also be grouped in a single PA probe, as it was developed in the European Union project FULLPHASE [Beckmann 2014] and for the commercial system Vevo LAZR (Visualsonics,Fujilm,Figure 1.10(b) [Needles 2013]) which is largely described later in this manuscript.

Depending on the targeted application, the desired resolution and the imaging depth, one of these techniques is used to acquire multispectral photoacoustic images. Specic probes are also developed for particular applications, as endoscopy probe for applications in interventional surgery, gastroenterology or foetal medicine [Ansari 2015], and combined with one of these systems.

Examples of biological and medical multispectral photoacoustic applications

PA imaging is of great interest for a large range of applications where functional information of biological tissue are required. From interventional surgery to drug delivery, passing through biological media detection or quantication (lipid, blood etc.), multispectral photoacoustic imaging is useful.

Interventional surgery requires discriminating between lots of biological tissues to avoid any mistake during the surgical intervention. Studies on the use of PA imaging to discriminate, e.g., between nerves and tendons, are investigated [Mari 2014, Mari 2015]. Indeed, these two biological tissues cannot be discriminated based on ultrasound imaging which is most of the time used for interventional surgery.

12 Aneline DOLET 1.5. Current researches in multispectral photoacoustic imaging The larger application domain of photoacoustic imaging is based on the access to Hb and HbO 2 concentrations which allows the calculation of the blood oxygenation rate. This is of great interest for numerous applications: placental oxygenation investigation [Arthuis 2017], oxygen saturation rate calculation (sO 2 ) [Deán-Ben 2014], dierentiation between malignant or benign tumors as part of the follow-up of carcinomas [Bauer 2010] or to evaluate tissue death [Su 2012].

The photoacoustic eect can also be used for drug delivery [Dixon 2016]. Indeed, using light to heat a structure can allow the release of some drugs rstly contained inside this structure. The main interest of photoacoustic imaging in this case is that this structure can be imaged and then tracked using average energy. It thus allows drug release, using more optical energy to heat the structure, at the exact wanted position.

Image processing methods

Each application and system requires particular image processing method to accurately assess the desired information. Indeed, the photoacoustic signal of a medium varies depending on (1) the image depth, (2) the sizes of the imaged structures, (3) the medium concentrations and (4) the acquisition system. The discrimination or quantication of media using multispectral photoacoustic imaging is then challenging and numerous strategies have been developed, which can be supervised or unsupervised, to cluster or unmix media [Cox 2012].

The algorithms mentioned in the current section are separated in two groups belonging to the targeted results. Indeed, it is important to understand the dierence made between the discrimination of data, using clustering algorithms, and the quantication of data, using unmixing algorithms. Considering a Aneline DOLET 13 Chapter 1. Photoacoustic imaging photoacoustic image acquired on a biological tissue composed of several dierent media, the discrimination aims at nding if a pixel belongs, or not, to a medium. Each pixel of the data is thus labeled and each label corresponds to a single medium. On the other hand, the quantication of data aims at assessing the concentration of the imaged media, which means that each pixel can be considered composed of (1) a single pure medium, (2) a dilution of a medium or (3) a mix of media at dierent concentrations. The unmixing algorithms aim at decomposing each pixel according to dierent criteria which have to be of interest for the targeted application. Also, the dierence between supervised and unsupervised methods has to be claried. Supervised algorithms need to compare the data to existing references. In our context, the reference data could be the reference spectra acquired on the media to cluster. The algorithms then label each pixel to the closest reference. Dierent meaning of closest can be taken into account using dierent strategies. On the contrary, unsupervised algorithms process the data without additional information.

Clustering methods

Dierent clustering algorithms have been proposed in the literature to discriminate multispectral PA data. Supervised algorithms based on least-square minimization [Kruizinga 2014] or intra-class correlation [Mallidi 2008[START_REF] Wang | [END_REF]] have been proposed for application in the detection of atherosclerotic plaques. Methods based on wavelet-packet features, with training and testing steps, have been also developed [START_REF] Zalev | [END_REF]]. For applications where only two dierent media have to be discriminated, a limited range of wavelengths can be chosen as studied by [Mari 2014, Mari 2015] for the discrimination of tendons and nerves. Unsupervised methods have been also proposed to segment photoacoustic data [Bauer 2010] and to analyse the media regarding a specic peak in their spectral evolutions [Daeichin 2016].

Unmixing methods

A supervised unmixing method called spectral-tting and two unsupervised methods called principal component analysis and independent component analysis have been compared for the detection of ICG and Cy7 inclusion in biological tissues [Glatz 2011]. An unsupervised unmixing method to calculate sO 2 was implemented in the commercial system Vevo LAZR (Visualsonics, Fujilm) [START_REF] Zhang | [END_REF], Needles 2013, Arthuis 2017]. The method solves an inverse problem considering the molar extinction coecient of Hb and HbO 2 .

All these developed methods highlight the diculties to nd a unique strategy, allowing the discrimination or quantication of media in multispectral photoacoustic imaging, which could be applied for various applications and/or acquisition systems.

Motivations of the thesis 1.6.1 Objectives

The large range of acquisition systems and applications requires developing robust processing methods which can be applied on several dataset, without acquisition system dependence. As the calculation of blood oxygenation is of great interest for dierent application, the major objective of the present thesis is the development of a multispectral PA image processing chain for its calculation in biological tissues.

The main steps are, rst, the data discrimination (clustering), to extract the ROI, and second, the quantication of the dierent media in these ROI (unmixing). 14 Aneline DOLET

Manuscript organization

The proposed method has to be tested on various acquisition systems to validate its robustness. The proposed clustering and unmixing algorithms are then compared on data acquired with three dierent systems: the PACT platform at CREATIS, a PAT system (collaboration with Prof. Zhen, MACAU) and a Vevo LAZR commercial system. Finally, a real-time photoacoustic mode should be developed to improve the experimental platform acquisition at CREATIS.

Targeted applications

The targeted application is the quantication of oxygenation in vascularization since it is of great interest for various biomedical applications. In the present thesis, two of these applications are focused on: the discrimination of malignant or benign tumors and the dating of tissue death. As the calculation of saturation in oxygen ratio (sO 2 ) has already been implemented on a commercial system, the develop method is validated comparing its result to the one given by the commercial system. The study on the dating of tissue death is done on dierent biological tissues and compare to the literature [Su 2012].

Manuscript organization

In this Chapter 1, the photoacoustic imaging concept has been presented as well as the state of the art and the thesis objectives. According to the scheme colors of Figure 1.11, the set-up and dataset are described in Chapter 2, as well as the co-supervision developments presented in section 2.4. The clustering and unmixing developments and validation are presented respectively in Chapter 3 and Chapter 4 where the proposed method performances are compared to the literature method performances. Finally, the application tests and validation are shown Chapter 5.

Introduction

The interest of multispectral photoacoustic imaging for dierent applications conducted to its considerable development in the last decade. Several experimental and commercial systems have been developed

to acquire multispectral photoacoustic data for biomedical applications with various optical excitations and ultrasound devices for reception.

In this context, the dierent processings developed during this PhD thesis have been validated on dierent experimental and commercial platforms. This Chapter presents the materials used to acquire the photoacoustic dataset. The developments done to improve the experimental platform of CREATIS, in collaboration with our partner MSDLab (Florence, Italy), are highlighted and the dataset used for the validation are nally presented.

Aneline DOLET Chapter 2. Multispectral photoacoustic imaging: systems, imaged samples and dataset

Experimental platforms

Three dierent multispectral photoacoustic experimental platforms have been used during this thesis.

The main dierence among all these set-up is the geometry between laser illumination, imaged samples and ultrasound transducer/probe. These platforms and their respective geometry are presented in this section and summarized in Appendix A.

Photoacoustic tomography (MACAU)

A multispectral photoacoustic tomography (PAT) acquisition set-up [Li 2015] was available during To acquire signals with enough photoacoustic energy, the absorbing media are set close to the surface as well as the ultrasound transducer to avoid too much optical energy attenuation. This system gives access to RF data and a classical DAS algorithm is used for the reconstruction of PA tomographic images.

Experimental platform (CREATIS)

State of the art

An experimental multispectral photoacoustic platform is available for acquisition in CREATIS In comparison with the tomographic set-up at MACAU, the platform of CREATIS allows the imaging of samples (e.g., phantoms or biological tissues) of dierent shapes and sizes. RF data are accessible and the images are reconstructed with DAS algorithm. However, as the laser energy is not perfectly stable over time, twenty images must be acquired at each wavelength. The dataset is then the average of these twenty images for each wavelength.

Improvements

The above description corresponds to the initial state of the CREATIS photoacoustic system. The experimental set-up has been improved during the present PhD. As the laser beam hitting the sample is less than 5mm diameter, the possible imaged area is really small. To optimize the optical energy, tests were conducted using dierent lenses to enlarge the region of illumination. The objective is to enlarge the laser beam while keeping enough optical energy to create detectable PA signal. The initial optical energy needs to be high enough to create photoacoustic signal up to a depth of 1cm minimum taking into account optical attenuation. Finally, to optimize the optical energy, optical bers have been installed at all the laser beam output (one for output 1 and one which can be used for output 2 and 3, Figure 2 (c) Sample inside the system with the bimodal US/PA probe ready for acquisitions.

Phantom manufacturing

Data and image processings developed in this work needed to be validated on reproducible, calibrated and stable samples. Consequently, a large variety of phantoms has been made with specic geometry and property.

Aneline DOLET Chapter 2. Multispectral photoacoustic imaging: systems, imaged samples and dataset Cylindrical phantoms made in Macau for tomographic acquisitions are not presented here, a detailed description can be found in [Li 2015]. We will focus hereafter on home-made phantoms designed for acquisitions on the experimental platforms in Lyon. These phantoms and their developments are presented in the next section.

Materials

The work of Maeva VALLET [Vallet 2015] has shown that both phantom in agar (4%) or polyvinyl alcohol (PVA) have interesting acoustic and optical properties to mimic biological tissues. During my PhD work, I preferred to develop phantoms in agar (4%) because the design and the production is much faster than PVA and optical properties are easier to change.

Agar material is a powder that needs to be melt with water to create a gel. In our case, 4% in weight of agar are put in a volume of water. First, the volume of water is heated up to 65 • C. At this temperature, the agar powder is put in the water and stirred together. Finally, this mixture is heated up to 80 • C to thicken it. The gel is cast in a mold, designed depending on the desired shape and size (Figure 2.6), and put to the fridge until stiening (some hours, depending on the phantom thickness). Colored inks have important absorption in the visible spectral domain (≈ 400nm -700nm). In order to produce calibrated photoacoustic phantoms, it is essential to accurately know their absorbing optical properties. To this end, a study of these inks was conducted to compare the absorption evolutions over wavelengths, using (1) spectrophotometry, and (2) multispectral photoacoustic imaging [Dolet 2017].

Previously, a similar study had been conducted in [Cai 2011] but with a smaller range of wavelengths and less optical absorbers.

Multispectral photoacoustic imaging

Multispectral photoacoustic signals are acquired on the experimental photoacoustic platform in CRE-ATIS. In this study, the CMUT probe is used for the US detection. Acquisitions are made from 470nm to 690nm, with 5nm steps on colored 4% agar phantoms with a dilution factor of 1:1629 (example of phantom and acquisition in Figure 2.8). To make colored 4% agar phantoms, the same procedure as the one presented before is used. The desired volume of ink is added directly after mixing the agar powder with water at ≈ 65 • C. It is necessary to carefully mix the ink inside this mixture to have homogeneous phantoms in terms of optical absorbers.

Comparison method

The data acquired with both modalities cannot be compared without previous cross-calibration. The absorbances, or spectrophotometer spectra, measured with the spectrophotometer have rst to be postprocessed to make them comparable to the multispectral photoacoustic signals, or photoacoustic spectra.

First, the photoacoustic laser energy has to be taken into account as it is not constant over wavelengths nor over time. As the absorbance spectra of black ink presents a decreasing exponential shape (Figure 2.9(a)), an optical signal is available for all of the wavelengths of the selected spectrum. Diluted black ink is then used to calculate a calibration ratio (equation (2.1)), to normalize the spectrophotometer spectra and compare them to the photoacoustic spectra:

Aneline DOLET

Ratio = A black S black (2.1)
where A black is the multispectral photoacoustic signal and S black is the absorbance measurement of black ink. The ratio is applied to the ve other absorbance spectra, corresponding to 5 dierent inks:

S P A = S × Ratio (2.2)
where S P A is the normalized absorbance measurement, which can be compared to the multispectral PA signal, and S is the initial absorbance measurement. Using this ratio, the photoacoustic laser energy is taken into account, assuming a constant energy distribution over all of the photoacoustic acquisitions.

The small dierence of ink concentrations between spectrophotometer and multispectral photoacoustic experiments (1:1601 vs 1:1629) is not taken into account. Meanwhile, a gain factor is applied to S P A in order to minimize the root-mean square error (RMSE) between both measurements. This procedure focuses on the spectrum qualitative shapes, but not on the quantication of the absorbance. RMSE is calculated as follows:

RM SE = L i=1 S P A -A 2 (2.3)
where L is the number of measured wavelengths. For each ink, the similarity between normalized spectrophotometry and multispectral photoacoustic imaging can then be analyzed.

Results

In this study, six inks were used with black ink as reference. The comparison between spectrophotometry and multispectral photoacoustic signal is done on the ve phantoms (brown, yellow, red, blue and green inks). Note that the spectrophotometry absorbance for blue and green inks exhibit logically a peak at 590nm and 600nm respectively. In the same time, in the range below 525nm for blue and 580nm for green, the measured signal by spectrophotometry is very low and the applied ratio is high.

We decided to restrict the study ranges above these wavelength values (> 525nm for blue and > 580nm

for green) avoiding to increase the noise. Aneline DOLET Chapter 2. Multispectral photoacoustic imaging: systems, imaged samples and dataset

Phantom manufacturing

As shown in Figure 2.9, the dierent ink absorbances follow the same tendency as the multispectral photoacoustic signals (Figure 2.9, solid lines and * markers, respectively). The blue ink, however, shows a larger discrepancy. RMSE is calculated for all the inks to quantify this discrepancy and are given for each ink on Figure 2.9. The best similarity between both measurements is for brown ink (RMSE of 0.09) and the worse one is for blue ink (RMSE of 7.39). Indeed, for the blue ink, the absorbance peak is much larger for the PA signal than for the spectrophotometry signal. The red ink also gives a high RMSE (5.58) which comes from a large discrepancy for the smallest wavelengths. With the yellow ink, there is a dierence in the decay rate of the absorbance between both of these modalities. However, for all the media, the agreement is largely comparable to the results reported in literature [Cai 2011], with both of these modalities giving similar trends in absorbance for the same medium.

It is of major importance to accurately know the spectral properties of optical absorbers used to produce calibrated phantom for photoacoustic. Indeed, the selected optical absorbers for setting a PA designed to discriminate or quantify media depend on: (1) the spectral properties of optical absorbers and

(2) the wavelength range used for PA imaging acquisition. The interest of this study is that a medium absorbance evolution over wavelength is often known, or easily measurable, while the multispectral PA evolution is harder to get. This study conrms then that optical absorbers, for photoacoustic phantom design, can now be chosen regarding their absorbance evolution over wavelengths [Dolet 2017].

Multi-colored phantoms production

The processings developed during the present thesis aim at discriminating dierent media. Their validations need to be done on phantoms created with dierent optical absorbers. The easiest way to realize that is to produce a multi-colored phantom (e.g., bi-colored phantoms are shown Figure 2.10).

This type of phantoms is created in three steps (explained here for bi-colored phantoms): (1) a rst phantom of a single color is done, (2) it is cut in two parts and one of these is put again in the same phantom mold, and (3) a 4% agar mixture is colored with a second color and put in the empty part of the mold. After a few hours in the fridge, a bi-colored phantom is ready with no air between both colors to avoid any problem for ultrasound propagation. Firstly, I created spherical inclusions with a mold already used to create phantoms for US studies (Figure 2.11(a)). This mold allows the creation of spherical inclusions with various diameters: 20mm, 15mm, 10mm and 5mm. To make these inclusions, colored 4% agar mixture is made and injected in the mold with a needle (Figure 2.11(b)). When the inclusions are full, the mold is put in the fridge during 2 hours, which is enough as the inclusions are small. After unmolding the inclusions, they are xed inside a rectangular mold. To this end, some holes have been made in the mold to insert a nylon wire which holds the inclusions (Figure 2.11(c)). Uncolored 4% agar is made to ll the mold all around the inclusions.

After stiening in the fridge, the phantom is ready ( The spherical inclusion production is long and needs a lot of attention. A single mold is available, only one inclusion of each size can then be done at the same time. Each color has to be done at two dierent times, with hours in the fridge between. Also, the inclusions are small, it can then be a hard and long work to x them. It happens thus that the wire cut the inclusion, principally the smallest ones.

Aneline DOLET Chapter 2. Multispectral photoacoustic imaging: systems, imaged samples and dataset

To go further to validate the processings developed during this PhD, dierent shapes and smaller inclusions need to be created. Indeed, the laser beam is ≈ 5mm diameter large which depends on the selected wavelengths. Imaging two spherical inclusions of 5mm diameter is then impossible. Moreover, as the smallest inclusions are the hardest to hold, other sizes are most of the time used, and it is hard to excite a region where both spherical inclusions are located.

Improvements

I designed a new mold to improve each of the points mentioned above. 8 rectangular inclusions of each size can be done allowing the creation of dierent colored inclusions of the same size at the same time. It allows the creation of rectangular shape inclusions of 10mm, 8mm, 5mm and 3mm width and 35mm high (Figure 2.12(b)). Some big inclusions (from 5mm to 10mm) are kept available as the 5mm are already hard to x in the mold. The 3mm ones were only a test to know if this size could be used.

The mold had to support ≈ 80 • C, to put the hot agar inside, but also, the temperature of a freezer.

Indeed, for further photoacoustic developments, PVA could be later used and this material needs to go to a freezer during some hours. The mold has then been done in acrylic resin material. The DAS algorithm already developed on the scanner has been adapted to this context.

3D real-time imaging

The acquisition of the new ULA-OP scanner at CREATIS also makes 3D photoacoustic imaging possible. A 3D photoacoustic real-time mode was then implemented on the ULA-OP256 scanner. As the 2D mode was already developed, it was used and adapted to create the 3D one. A matrix prototype array probe (Vermon, Tours, France) is used. The probe is composed of 8 × 24 elements with 3.84M Hz center frequency and a bandwidth of 65.1% at -6dB. The ULA-OP256 scanner is programmed to acquire ultrasound images using diverging waves. For passive and active acquisitions 192 lines are used. The images are reconstructed with a 3D DAS algorithm.

Technological transfer to CREATIS

The hardware developments presented above have been conducted in the MSDLab (Florence, Italy).

It was after installed on the experimental platform in CREATIS. On this platform, the laser cannot be triggered by the ultrasound scanner so the inverse is done. For ULA-OP64, a generator is used to trig the ULA-OP64 with the laser. The same has to be done for the ULA-OP256 scanner.

For the real-time mode, the required triggered signal is presented in red in

2D real-time imaging

To validate the developed 2D real-time PA mode, the ULA-OP256 system is connected to the linear probe LA523E. The imaged phantom is shown Figure 2.17(a). It is an uncolored 4% agar phantom with two included graphite mines. The yellow rectangle on Figure 2.17(a) highlights the imaged region and the orange ellipse, the laser beam size. The phantom is excited by a laser beam at 1064nm (Figure 2.17(b), the visible green light is due to a residual laser component at 532nm). The phantom was translated from left to right by maintaining the plane of the graphite mines aligned with the imaged region. The photoacoustic signal is in blue while the active B-mode image is in black and white (Figure 2.17(c)).

Some artefacts are visible on Figure 2.17(c), particularly at the second phantom position (t1). They are due to the photoacoustic signal multiple echoes. Because the ultrasounds do not well pass through the graphite mines, the photoacoustic signal of the bottom graphite mine is lower than the one of the other mine (Figure 2.17(c)-t1). The validation is done using the Vermon matrix array and a gold particle (Figure 2.18(a)) inserted in an uncolored 4% agar phantom. The gold particle has a cylindrical shape of 1mm diameter by 3mm length. Also in 3D, US and photoacoustic images are superimposed in real-time. The 3mm axis is shown 3D real-time imaging could be improved to cancel the artefact that can be seen in Figure 2.18(c-d).

Regarding now the resolution with a 6dB dynamic (Figure 2.18(c-d)), the resolution of PA imaging is a little better than the one of US imaging. The lateral resolution of US imaging is of 3.5mm and 1.7mm

in the length and diameter planes, respectively. The length and diameter planes lateral resolutions are, for PA imaging, 3.2mm and 1.4mm, respectively. The axial resolution is quite similar between both modalities: 0.7mm for both techniques in the length plane and 1.7mm and 1.4mm in the diameter plane for US and PA images, respectively. The 3D resolution is then consistent with the particle dimension.

Multispectral photoacoustic dataset

To validate the developed processings in this thesis, various dataset will be used. They could be entirely composed of optical absorbers, to have a large ROI, or partially, with small ROI like inclusions, which can be spaced or close. A large number of dataset has been used during these three years, but only the ones used to present the results in this manuscript are presented.

Optical absorbers uniformly distributed in the phantom

The acquisition system used to acquire the dataset presented in this section is the Vevo LAZR.

Acquisitions are done using the full available range of wavelengths ([680nm ; 970nm] with 1nm steps) as it does not require too much acquisition time (≈ 1min). Classicaly in multispectral PA imaging, between 5 and 10 wavelengths are used as it is enough to discriminate or quantify media with an acceptable processing time. 8 of the acquired wavelengths are then selected (from 680nm to 820nm, 20nm steps) to compose our dataset. The selected wavelengths were chosen in a range where the used optical absorbers (blue and green inks) can be discriminated as their spectral evolutions are dierent.

The colored phantoms used to acquire the dataset are composed of three dierent 4% colored agar parts. The left one is done using only blue ink (400g of water, 16g of agar and 380µL of ink) and the right part, only green (400g of water, 16g of agar and 950µL of ink). Both parts are considered as relative ink concentrations equal to 1 (i.e., pure media). Blue and green ink quantities are not equal because both inks do not exhibit the same maximum photoacoustic signal amplitude. Based on the previous presented study, the ink quantities have been chosen to obtain the same maximum photoacoustic signal amplitude for both pure media.

The central part of the rst described phantom is composed of a 0.53 dilution of blue ink relatively to the pure considered concentration presented before. The concentration is not perfectly equal to 0.5 because of the 10µL pipette accuracy used to measure the ink dose (with a minimal possible dose of 20µL).

It is composed of 100g of water, 4g of agar, 50µL of blue ink. This phantom is shown Figure 2.19(a) and is referred as B-Bdil-G(Vevo) for the rest of the manuscript.

Aneline DOLET Chapter 2. Multispectral photoacoustic imaging: systems, imaged samples and dataset 2.5.2.2 Inclusions close to each other

The dataset now presented has been acquired on the commercial system Vevo LAZR. The imaged phantom is a small 4% agar phantom done in the rectangular inclusion mold. The inclusion size is 10mm

width. The acquisitions were done putting the phantom in water.

Two 5mm rectangular inclusions have rst been done, one blue and one green, using the concentrations corresponding to the pure media described before. These inclusions are put inside the 10mm inclusion mold at two opposite corners (Figure 2.22). The rest of the inclusion is completed with a mix of 0.74 blue and 0.29 green (100g of water, 4g of agar, 70µL of blue ink and 70µL of green ink) relatively to the considered pure media. This dataset is referred as IncB-mix-G(Vevo) in the following Chapters.

3D imaged region dataset

Only two, out of the three systems used during the present thesis, allow 3D ROI imaging. CREATIS experimental platform with the matrix array probe Vermon and the Vevo LAZR commercial system using a mechanical sweeping-translation of the LZ400 probe. The dataset presented in this section was acquired on the Vevo LAZR commercial system with 0.1mm mechanical sweeping-translation steps.

38

Aneline DOLET

Conclusion

The imaged phantom is an uncolored 4% agar phantom with two 5mm diameter spherical inclusions, one fully blue and one fully green, using the same relative full concentrations as for the other phantoms (Figure 2.23). This dataset is called 3D(Vevo) for its use in the next Chapter.

Conclusion

In this Chapter, the three photoacoustic systems used during this thesis have been presented. They are based on two dierent photoacoustic acquisition techniques: PAT and PACT. The acquired images have then dierent resolutions and the possible imaged ROI sizes can vary. To validate the developed processing presented in the next Chapters, phantoms and inclusions of various sizes and shapes needed to be produced. The developments done for their production have been highlighted. As the laser beam of CREATIS photoacoustic experimental platform is narrow with low optical energy in the interesting wavelength range, technological developments have been done to well image the ROI. These developments of real-time 2D and 3D PA/US modes and their installation to facilitate the acquisitions with the CREATIS system have been described here. Meanwhile, the narrow beam still limits the exploitation of data of this system.

Aneline DOLET

Introduction

Processing of medical imaging data usually requires a step called segmentation to identify regions of potential interest. In the case of multispectral photoacoustic imaging, this step is called clustering. The main goal of clustering is to group with the same label all the pixels having the same spectral properties.

The unsupervised proposed method for multispectral PA clustering will be presented as well as the most used supervised methods of the literature.

Aneline DOLET Chapter 3. Clustering of multispectral photoacoustic dataset Each pixel of multispectral photoacoustic data x i located at a position s i has a spectrum A i . The N samples of the ROI, x i , are expressed as:

x i = [s i , A i ] ∈ X with s i ∈ R S , the spatial position A i ∈ R L , the spectral prole i ∈ [1 ; N ], the sample index (3.1)
The aim of the proposed pre-processing is to discriminate the background from the ROI. The background is determined by applying a threshold related to the gradient calculation from a Sobel ltering.

The Sobel threshold T h Sobel is computed as follows:

Methods of the literature

T h Sobel = 2 * G (3.2)
where G is the mean of the Sobel gradient magnitude [Pratt 1978]. The edges detected by the Sobel ltering are then only used to calculate G. T h Sobel is applied on the sum over wavelengths of the initial dataset to create a binary mask B of the same size as the PA data at a single wavelength: the background pixels are set to 0 while the ROI pixels are set to 1. The following clustering methods are only applied on the pixels out of the background.

Two strategies to deal with concentration dierences

Depending on the targeted application, it might be important to discriminate dierent concentrations of a single medium, e.g., to determine the concentration of a contrast agent in the body [Mienkina 2009], or to discriminate a single medium from all other media without considering its dilution, e.g., to determine the level of vascularization for calculation of the oxygenation rate [Mercep 2017]. To allow these two strategies, a normalization of the spectra is applied on each pixel of interest when no discrimination between medium concentrations is needed. Indeed, the spectral shape of a medium is independent of the concentrations of its constituents. To consider only the spectral shape without including the PA amplitude, the spectra of each pixel is normalized as follows. Let us consider the spectra, A i ∈ R L , of a given pixel, where a i,λ l is the discrete photoacoustic signal amplitude which is acquired at the wavelength λ l , l = 1 • • • L, chosen according to the media for discrimination. We rst integrate and normalize the amplitudes of A i , as:

a * i,λ l = l j=1 a i,λ j L j=1 a i,λ j (3.3)
where a * i,λ l is the normalized amplitude of a i,λ l at the wavelength λ l . This processing is applied to all of the pixel spectra of the dataset matrix A and at all of the wavelengths λ l of [λ 1 ; λ L ]. The normalized dataset has values from 0 to 1. Before any other processing, the dataset values are placed in the range [0; 255] with or without this normalization step. This normalization by integration over the wavelengths is preferred to normalization by the maximum value of the spectra. Indeed, the proposed normalization smooths the noise, whereas the maximum normalization does not. On one hand, applying this normalization step to the data allows dierent concentrations of a single medium to be merged as a single cluster. On the other hand, using non-normalized data allows discrimination between dierent concentrations as if they were associated with dierent absorbing media.

Methods of the literature

Segmentation or classication methods have been developed to discriminate between biological media using multispectral photoacoustic dataset. Supervised methods, such as, spectral-tting [Glatz 2011], approaches based on the least-square criterion [START_REF] Jansen | [END_REF]], and intra-class correlation [START_REF] Wang | [END_REF], Mallidi 2008] have been proposed in the literature. These methods imply to know a reference spectrum of each medium to identify during the clustering process.

Aneline DOLET Chapter 3. Clustering of multispectral photoacoustic dataset

Endmembers calculation

Supervised clustering methods are based on a similarity calculation between a reference spectrum and each dataset pixel spectrum. Most of the time, the reference spectrum of each constitutive medium is not known, particularly for biological tissues. Indeed, even if the absorption spectra of some biological components can be known (Figure 1.6), the real shape depends on the photoacoustic system used to acquire the dataset. So, a reference spectrum, also called endmember, can be computed directly on the acquired data for each constitutive medium, as the mean of its pixels spectra.

The k endmembers corresponding to the k constitutive components of the imaged area are calculated from the binary mask B (section 3.2.1). B is manually labeled with a single label for each constitutive medium. For each of the k labeled regions, the mean spectrum is calculated (endmember). The whole endmembers are gathered into a matrix E c of size k × L (k endmembers × L acquired wavelengths).

Spectral-tting method

The spectral-tting (SF) method is here used as a clustering method like described in [Glatz 2011] to discriminate ICG from Cy7 injected in biological tissues. SF aims at nding which of the endmember best ts the spectrum of a given pixel with a least-square criterion. The multispectral PA dataset matrix A of N × L measurements (N is the number of pixels in the imaged region and L the number of acquired wavelengths) and the endmember matrix E c are considered. The Moore-Penrose pseudoinverse matrix

E +
c is rst calculated as follows:

E + c = E T c (E c E T c ) -1 (3.4)
The coecient corresponding to the probability, in the spectral-tting sense, for a pixel to be part of a cluster is given by the SF matrix, calculated by the following matrix inversion:

SF = AE + c (3.5)
SF contains N × k coecients that are called abundance coecients as they represent the tting of the data with each of the endmember. To nally cluster the data, the maximum abundance coecient value over the k endmembers in SF , m i , as well as its position, e i , which corresponds to one of the endmembers is found for each pixel. In the clustering result matrix of the spectral-tting method, C SF of size N × 1, the i-th pixel is labeled by the value e i . In equation (3.6), SF ie is a coecient of SF with i referring to a pixel and e to an endmember. The least-square (LS) minimization algorithm used in the present thesis is inspired by the study [START_REF] Jansen | [END_REF]]. Using the same notations as before, the least-square distance ls ie between the i-th pixel spectra, A i , and the e-th endmember, E ce , is calculated as follows:

     [m i , e i ] = max e∈[1;k] (SF ie ) C SF i = e i
ls ie = L l=1 (A i -E ce ) 2 (3.7)
All the calculated least-square distances, ls ie with i ∈ [1 ; N ] and e ∈ [1 ; k], are got into the matrix LS of size N × k. Finally, the clustering result matrix C LS of size N × 1 is constituted of the endmember label e i which corresponds to the minimum ls ie coecient of the i-th pixel:

     [s i , e i ] = min e∈[1;k]
(LS ie )

C LS i = e i (3.8)
with s i the smallest value of the vector [ls i1 . . . ls ik ] and e i the corresponding endmember label.

Intra-class correlation method

The intra-class correlation (ICC) is a method which uses both the spatial neighborhood and the spectral feature of the pixels to compensate the tissue movements and the optical energy changes between photoacoustic acquisitions at dierent wavelengths [Mallidi 2008[START_REF] Wang | [END_REF]. As this method also uses spatial features, better clustering results were expected. ICC aims at calculating the strength of the linear relationship between the pixel spectra of the dataset and the endmembers. To take into account the spatial features, the considered spectrum corresponds to the mean of the pixel spectra, in a neighbor area of the i-th pixel A i , calculated as follows:

A i = 1 S S j=1 A j (3.9)
with S the number of pixels in the considering neighbor area and A j , j ∈ [1; S], these pixels spectra.

The intra-class correlation coecient icc ie between the i-th averaged pixel spectrum, A i , and the e-th endmember, E ce , is calculated as follows:

icc ie = 2 × L l=1 (A i -X)(E ce -X) L l=1 (A i -X) 2 + (E ce -X) 2 (3.10)
where X = A+Ec Even if the intra-class correlation used both spatial and spectral features, this method did not provide the expected performance improvement. The reference spectra calculation clearly impacts the results of the three previously presented methods. An unsupervised clustering method using the spatial and spectral features is then proposed. This, compared principally to intra-class correlation, should avoid the varying performances related with the reference spectra calculation.

To this end, the spatiotemporal mean-shift approach (STM-S) was adapted to our context. This algorithm demonstrated good performances for magnetic resonance image clustering by taking into account both spatial and temporal features [START_REF] Mure | [END_REF]].

Spatio-spectral regularization

Based on the clustering algorithm developed for magnetic resonance imaging [START_REF] Mure | [END_REF]], we proposed here a spatio-spectral mean-shift (SSM-S). It means that the spectra of multispectral PA data are handled instead of the temporal evolutions of longitudinal magnetic resonance images.

This algorithm requires two parameters, R S and R λ , that are now presented. For a given pixel x i , let us consider rst its neighboring pixels in the spatial dimension within a radial distance R S (Figure 3.2(a)).

These pixels x j = [s j , A j ] have to satisfy the equation:

1 R S 2 (s i -s j ) T (s i -s j ) ≤ 1 (3.12)
As the sample spectra of a same medium are close, R λ is xed as the maximum accepted distance between the spectra of two samples of the same medium:

1 R λ (A i -A j ) ∞ ≤ 1 (3.13)
In equation (3.13), the innity norm is the most appropriated norm, because it allows pixels with close spectral shapes to be merged even if the amplitudes at each wavelength are not exactly the same.

This also allows spectra with large dierences in amplitude at a single wavelength to be dierentiated, because the spectral shapes of these pixels are then considered dierent. However, dierent norms could be used, depending on the desired eect.

46

Aneline DOLET 3.4. Spatio-spectral mean-shift x

[t+1] i = N j=1 g s (s [t] i ; s [t] j ) • g λ (A [t] i ; A [t] j ) • x [t] j N j=1 g s (s [t] i ; s [t] j ) • g λ (A [t] i ; A [t] j ) (3.14)
where g s (s

[t] i ; s [t] j ) = 1 if equation (3.
12) is respected, or it is 0 otherwise, and g λ (A

[t] i ; A [t] j ) = 1 if equation (3.13
) is respected, 0 otherwise. These steps are applied iteratively to all pixels x i of the imaged region X until convergence of the procedure; i.e., stabilization of the featured value throughout the ltering process. At each iteration, the shift of each updated pixel during the regularization process is calculated. The algorithm stops when the sum of all these shifts becomes lower than a pre-dened threshold. The result of this SSM-S algorithm can easily be used to segment the imaged region. Even if the result contains too many clusters compared to the desired result, all of the pixels with close spectra are labeled together.

The main advantage of mean-shift spatial regularization is that two areas of the imaged region can be merged into the same cluster even if they are spatially distant, as long as their spectra are close.

Moreover, this helps two nearby pixels to be clustered together even if their spectra have small amplitude or shape dierences. Indeed, the same medium might be present at dierent locations within the imaged regions of interest, and two close pixels are likely to belong to the same medium.

The setting of both parameters R S and R λ is presented in section 3.5.5.

Aneline DOLET

Cluster number reduction

With SSM-S, a medium is represented by more than one label because of light attenuation and ultrasonic dispersion impacting the amplitude spectra of a single medium. This is due to the strong amplitude dierences between these labeled spectra, which are obtained even if their shapes are similar.

Here, only the k largest labeled areas are kept, while k is the number of media to be discriminated.

The other pixels are kept unlabeled. Thus, the labeled areas correspond to the dierent media to be discriminated. To reduce the number of labels, spatial features post-processing is applied in the nearest neighbor manner. The sets of P labeled pixels and N -P unlabeled pixels can be referred to as Y = {y i } i=[1;P ] and Z = {z j } j=[1;N -P ] , respectively, with:

   y i = [s i ; p i ] z j = [s j ; p j = 0] (3.15)
where p i ∈ R is the associated label of a pixel x i . For each unlabeled pixel z j ∈ Z, the closest spatial pixel y i * in Y is found, and p j takes the value of its corresponding index p i * , as follows:

i * | j = argmin i=[1;P ] s j -s i then p j = p i * (3.16)
As all of the background pixels are equal to 0 at each wavelength since the pre-processing is applied, these pixels are easily clustered together by the SSM-S. However, with this method, some pixels of the ROI can be clustered with the background if, for all wavelength, the PA signal is always below R λ .

The SSM-S processing steps are summarized in Figure 3.3.

Figure 3.3: The spatio-spectral mean-shift proposed pipeline.

3.5. Results and performances 3.5 Results and performances

Quantitative validation method

The binary mask B manually labeled (section 3.3.1) is used as the ground truth reference for wellclustered pixels. A quantitative measure WCP is calculated as the percentage of well-clustered pixels.

The number of well-clustered pixels is counted after the end of the considered algorithm procedure. WCP is expressed relative to all of the pixels to cluster. To compute this criterion, a pixel merged with the background is considered as well clustered. With this strategy, WCP ranges from 0% to 100%; i.e., from "no well-clustered pixels" to "perfect discrimination result".

The aim of this evaluation is to determine whether or not the method can be used to discriminate between dierent media without considering the background. Consequently, the proposed validation method is preferred to classical Dice similarity coecient (DSC) [START_REF] Dice | Measures of the Amount of Ecologic Association Between Species[END_REF]] validation because DSC is more restrictive than WCP. Indeed, the pixels clustered with the background make the performance decrease with DSC which might appears less appropriate for the targeted application.

Medium dilution

The dierent steps of all the clustering methods are highlighted for the Blood-Ink(PAT) dataset.

First, to remind the dataset, the photoacoustic signal acquired at 730nm is shown on Figure 3. 4(a). The sum over wavelengths of the noise removed by the pre-processing threshold can be seen in Figure 3.

Blood-Ink(PAT) without normalization step

First, the clustering is done on the dataset without normalization. Three media have to be discriminated: pure blood, diluted blood and ink (k = 3). The SF, LS and ICC results are presented in Figure 3.5.

The rst three lines correspond to diluted blood, pure blood and ink abundance maps, respectively. The last line of images shows the clustering result maps followed by the W CP performance values. Each column presents a method results: from left to right, SF, LS and ICC. These methods have similar performances (W CP SF = 76%, W CP LS = 71% and W CP ICC = 78%). ICC result is given after the optimization of the average area size S. This area is a square of R S-ICC side size (S = R S-ICC 2 ). Here, R S-ICC = 8pixels. In this case, ICC and its use of spatial feature seems to improve the discrimination result.

Aneline DOLET

Blood-Ink(PAT) with normalization step

The normalization step can also be applied aiming to cluster both blood inclusions altogether. Only two media have here to be discriminated: blood and ink (k = 2). The ground truth B labeling is presented in Figure 3.4(g) and highlights the desired discrimination result. Figures 3.7 and 3.8 present the literature method and SSM-S results when the normalization step is done, respectively. The SSM-S parameters are tuned with the same values as without normalization (R S = 15pixels and R λ = 19a.u.)

and R S-ICC is optimized with the value of 1pixel which means that the spatial features is not taken into account as it decreases the performance in this context. When the normalization is applied, R S-ICC is always optimized with the value of 1pixel.

SF, LS and ICC give really similar results (Figure 3.7). Indeed, the resulting maps cannot be distinguished with eyes and the W CP values are close (W CP SF = 85%, W CP LS = 87% and W CP ICC = 87%).

As R S-ICC = 1pixel, and compared to the results without normalization, it is normal that the ICC does not give better performance than SF and LS because the spatial feature is not taken into account. Aneline DOLET Chapter 3. Clustering of multispectral photoacoustic dataset

For Blood-Ink(PAT) dataset, the SSM-S method is the one reaching the best performances with and without applying the normalization step (see Table 3.1 which summarizes the used parameter values and the method performances). Almost all the pixels of interest are well clustered, only a few pixels at the blood inclusions boundaries are grouped with the ink cluster. Regarding the dataset (Figure 3.4(a)), all the inclusion boundary pixels exhibit really low photoacoustic signal. As the ink spectra are the ones with the lowest amplitudes (Figure 3.4(e)), it is not surprising that some of the boundaries pixels are clustered with this medium.

The dierent clustering steps and method principles have been highlighted on this rst dataset. For all the others, only the endmembers and nal result map for each clustering method is presented and discussed in this Chapter.

Blood-Ink(PAT) Normalization Parameters W CP (%) k R S-ICC R S R λ SF LS

B-Bdil-G(Vevo) dataset

The results presented in Figure 3.9 come from the application of all the clustering methods to the B-Bdil-G(Vevo) dataset. The dierent processings have been done using the optimized parameters: when no normalization is applied (R S-ICC = 3pixels, R S = 30pixels and R λ = 2a.u.), and with the normalization step (R S-ICC = 1pixel, R S = 30pixels and R λ = 4a.u.). These values are summarized in the Table 3.2. Figure 3.9(a-b) shows the endmembers used as references for the supervised clustering methods. It can be seen that without normalization, diluted blue and green reference spectra are really close as well as blue and green spectra with normalization. It highlights how the clusterization is challenging. As the normalization step aimed at clustering the diluted medium with the pure one, only two media (blue and green) have to be discriminated when the normalization is applied and only two endmembers are shown in Figure 3.9(b).

The results are presented as follows, Figure 3.9(c-d) shows the map results without and with normalization respectively. From up to bottom, the method result maps are presented: ground truth, SF, LS, ICC and SSM-S results. The performances are highlighted next to the result maps and summarized in Table 3.2. When no normalization is applied, the supervised methods assess similar performances, Both dataset results presented in this section test the behaviour of the algorithms on data with diluted medium. It has been highlighted that SSM-S give always the best performances in this context, with a W CP value superior to 86% and the result maps clearly allow to distinguish and localize the dierent imaged media. The pre-processing normalization step has also been tested and it well permits the merging of pure and diluted medium altogether. The proposed method has then been validated on diluted media for two dierent acquisition systems.

only

Mix of media B-mix-G(Vevo) dataset

The methods are now tested on phantoms with mix of media. This rst phantom is entirely composed of optical absorbers. Figure 3.10 presents the clustering results with the optimized parameters summarized in Table 3.3. When no normalization is done, R S-ICC = R S = 30pixels and R λ = 10a.u.

With the application of the normalization step, only R S-ICC value is modied and takes the value of 1pixel. In this context, the normalization aims at merging the mix with one of the pure media. As the green medium gives photoacoustic signal for every used wavelengths and the mix is principally composed of green (0.67), it should be clustered with the pure green medium.

The results are presented in Figure 3.10. Without normalization, the best performance is assessed by SSM-S algorithm (W CP SSM -S = 94%). ICC also well performs (W CP ICC is of 89%) and its result map (Figure 3.10(c)) clearly highlights the three dierent phantom parts. The two other methods have lower performances (W CP SF = 69% and W CP LS = 62%) as the mix and green medium parts are not well discriminated (Figure 3.10(c)). Regarding the green and mix endmembers (Figure 3.10(a)) which are really similar, the results of SF and LS are not surprising. This case shows the performance improvements using the spatial features with the high W CP values of ICC and SSM-S.

When the normalization is applied and the mix grouped with the pure green, all the supervised methods allows a well discrimination (W CP ≈ 93%). However, here the proposed method misclusters the mix and grouped it with the blue phantom part which leads to a really low W CP (38%).

IncB-mix-G(Vevo) dataset

This dataset is of particular interest as there are two dierent media in the axial direction, which means that the light is absorbed by some imaged media and also other media have to be imaged deeper.

It leads to important light attenuation, so photoacoustic signal decreases over depth. Moreover, there are high photoacoustic signals at the horizontal inclusion boundaries (Figure 2.22) that do not correspond to a particular absorbing medium. To test the clustering methods, these boundaries are manually removed.

Aneline DOLET The mix is here considered as a particular medium that has to be clustered alone, whether the normalization is applied or not. The results are presented in Figure 3.12. The used parameters are reported in Table 3.4. The supervised method results are similar without or with normalization application (W CP ≈ 52% and W CP ≈ 57%, respectively). The pure media are well clustered but the mix seems harder to discriminate. Regarding the endmembers and normalized endmembers (Figure 2.22(a-c)), the blue and mix reference spectra are really close which explains the algorithm diculties to discriminate these two media.

The proposed SSM-S method performs the best discrimination (W CP = 71% without normalization, W CP = 68% with normalization). Pure blue, pure green and one of the two mix parts are well clustered whether the normalization is applied or not (Figure 2.22(b-d)). However, the other mix part is clustered with the blue medium. It can be explained by the fact that the mix is principally composed of blue (0.7 of blue and 0.3 of green).

With this particular light attenuation, both mix parts could be considered as dierent absorbing media as they are not excited by the same optical energies. Using k = 4 and four endmembers in the processing chains clearly shows that it is the light attenuation which impacts the clustering performances.

Indeed, Figure 3.13(a-c) highlights the similarity of the blue and up mix parts spectra as well as the one of the green and bottom mix parts in both cases of normalization. This explained the diculties of the supervised methods to perform on this dataset (W CP ≈ 60% with or without normalization). Anyway, the SSM-S method assesses high performances in this context which highlights the interest of the proposed method (Figure 3.13(b-d)). It can be seen that with very small R λ (i.e., where too many clusters are created) or very large R λ (i.e., where the dierent media are clustered together), the performances decrease. However, the range where R λ is well adapted to the data is large enough (i.e., from 5 a.u. to 15 a.u. in the present case), and in this range, R S does not particularly aect the performance. Moreover, when R λ is above this range, the reduction of R S provided increased performance. It is explained by the fact that it helps the agglomeration of close pixels that are probably of the same medium. Finally, even if these settings change a little for each dataset, choosing R λ inferior to 20a.u. with R S in the range of 10pixels to 45pixels

IncB-mix-G(Vevo) Normalization Parameters W CP (%) k R S-ICC R S R λ SF LS ICC SSM-S (pixels) (pixels) (a.
provides good performances for all of the dataset, which makes the proposed algorithm easy to use.

Conclusion

This Chapter presents a comparative study of three dierent supervised clustering methods (SF, LS and ICC) and the unsupervised proposed method: SSM-S. The clustering is tested on dataset acquired on two dierent acquisition systems: the PAT and Vevo LAZR commercial systems. All the performances are summarized Table 3.5 where it can be seen that SSM-S gives the best result in each case except one particular, only because of the spatial post-processing. Indeed, even if the SSM-S procedure gives three main clusters well corresponding to the three phantom parts, the two main clusters kept by the post-processing are the ones of the mix and green parts and not the ones of both pure media.

62

Aneline DOLET

Conclusion

The proposed SSM-S method improves then the discrimination of multispectral photoacoustic dataset.

It has been validated on diluted media as well as on mix of media and particular light attenuation (light passing through two dierent absorbing media). Also, this Chapter highlights that using the spatial features, in addition with the spectral one, already improves the discrimination. Compared with ICC method, unsupervised clustering better performs as the endmember calculation dependence is overcome.

Considering 3D imaging, the SSM-S algorithm should be adapted to process the volume. This, instead of applying the procedure slice by slice, could reduce the misclustered pixels. Indeed, using the R S parameter in the 3 directions would allow to group more voxels with close spectra altogether.

Normalization Aneline DOLET

Parameters W CP (%) k R S-ICC R S R λ SF LS ICC SSM-S (pixels) (pixels) (a.u.) Blood-

Introduction

In photoacoustic imaging dataset, each pixel can be considered as a weighted sum, i.e a mix, of the imaged pure media where the weights represent the media concentrations. The pure medium spectra are called endmembers and the weights, the abundance coecients. The study of multispectral photoacoustic dataset unmixing, where the objective is to access these concentrations, is presented in this Chapter.

To unmix the photoacoustic data, regions of interest are discriminated from the background with the processing method presented in section 3.2.1. Unmixing methods are applied only to these regions of interest. In the previous Chapter 3, we have shown that supervised methods conducted to the calculation of reference spectra, which has a large impact on the results. Only unsupervised strategies will then be considered in this Chapter. First strategy is the use of unsupervised unmixing methods, from photoacoustic literature, that gives the abundance coecients. Second, unsupervised methods that only extract the endmembers are exploited and combined to a supervised but robust unmixing method, used with the endmembers previously extracted.

Aneline DOLET

Principal component analysis

As already used in [Glatz 2011], principal component analysis (PCA) can be used to unmix dierent media using multispectral photoacoustic imaging. This method is a blind source unmixing technique which assumes that the source components are statistically uncorrelated. It means that it is an unsupervised unmixing method which can be used only if the dierent imaged medium spectra are considered as statistically uncorrelated. We assume this hypothesis in our context.

PCA is a linear orthogonal transformation from the initial dataset to a new coordinate system where each axis corresponds to an uncorrelated principal component. The PCA algorithm gives: (1) the extracted uncorrelated principal component matrix which could correspond to the endmembers matrix and is then referred as E P CA and (2) the transformation matrix which corresponds to the abundance matrix U P CA . This is summarized by the following equation:

A = U P CA E P CA (4.1)
PCA extracts more components than the number of pure media in the imaged region. The endmembers kept in E P CA are the rst k extracted components which correspond to the k larger eigen values.

The E P CA matrix is then of size k × L. The U P CA abundance matrix considered is restricted to the k abundance coecients corresponding to the endmembers kept (U P CA of size N × k).

Independent component analysis

Independent component analysis method (ICA) was also used by [Glatz 2011] and is another blind source separation technique. The assumption considered here is that the sources are statistically independent, which is a stronger assumption than the one considered for PCA. We assume this hypothesis in our context.

The ICA method is a transformation from the initial dataset to a new coordinate system where each axis corresponds to an independent component. The ICA algorithm gives: (1) the extracted independent components matrix E ICA and (2) the abundance matrix U ICA . This is summarized by the following equation:

A = U ICA E ICA (4.2)
The independent components extracted with the ICA algorithm are more numerous than the number of pure media in the imaged region. The k larger independent components are kept in E ICA (size of k × L). The U ICA abundance matrix is restricted to the k abundance coecients corresponding to the endmembers kept (U ICA of size N × k).

In remote sensing eld 4.2.2.1 Rationale for the study

The methods of the photoacoustic imaging literature previously presented do not assess abundance matrix which allows an accurate calculation of medium concentrations. Mainly, the abundance coecient matrix U can have some negative values which is contrary to satisfactory concentration values. Indeed, an endmember concentration is equal to zero if the medium is absent in the considered pixel, equal to one if it is a pure medium pixel and in ]0 ; 1[ if the pixel is composed of the medium; e.g. diluted medium

or mix of media. The abundance matrix should then be constrained to the following requirements:

(1) abundance coecient values should be in the range [0 ; 1] and (2) for a single pixel, the sum of all the k coecients, corresponding to the k endmembers, should be inferior (diluted medium) or equal (pure medium or mix of media) to 1.

These kind of unmixing constraints are found in remote sensing eld where the imaged area contains, for example, road, building, sea, or forest that correspond to the media to unmix. The abundance, or concentration, is then the percentage area of a pixel containing a medium. Indeed, as the resolution is of some meters, dierent media can be found inside a single pixel. Both previously presented constraints are then also required in this eld.

As remote sensing unmixing is a huge and already deeply explored domain, numerous algorithms have already been developed and tested, which is not the case in our photoacoustic context. Moreover, the hyperspectral images dataset in remote sensing eld are close to the multispectral photoacoustic ones.

Both techniques image 2D or 3D region of interest at dierent wavelengths (from 2 to 20 in multispectral PA and more than 100 in remote sensing). These similarities motivated the study and adaptation of dierent unmixing methods, used in remote sensing, for their application to multispectral PA unmixing to accurately access medium concentrations.

Linear mixing model

Data unmixing requires solving a mixing model which can be linear or non-linear. In multispectral imaging, the interesting area is imaged at dierent wavelengths. Each pixel in the image is characterized by a spectrum which is the collection of acquired intensity values at each wavelength. A pixel's spectrum is either pure and considered as an endmember, or mixed (i.e., composed of a mixture of the endmembers).

The linear mixing model (LMM) [Keshava 2002] has then to be considered in our context considering a mixed pixel as a convex combination of the endmembers. More formally, it can be dened by the following equation:

A i = k j=1 u ji E j + g i , ∀ i ∈ [1; N ] (4.3)
where A i ∈ R L is the L-dimensional spectrum of the i-th pixel, k denotes the number of endmembers, u ji is the abundance of the j-th endmember in the i-th pixel, i.e. coecient of the matrix U , E j is the L-dimensional spectrum of the j-th endmember, g i is a vector of Gaussian white noise accounting for sensor noise and error of the model. All vectors are column vectors. The abundances, because they are contributions, must be positive and their sum has to be equal to one:

Aneline DOLET 67 Chapter 4. Unmixing of multispectral photoacoutsic dataset

u ji ≥ 0, k j=1 u ji = 1 (4.4)
The LMM is a simple but very representative model which was extensively studied in remote sensing, see for example the survey in [Keshava 2002]. For medical applications, the LMM can be a powerful tool for quantifying media in order to accurately assess the medium concentrations. In this context, a spectrum is assumed to be either fully concentrated, diluted or a mix of several endmembers. The LMM can be naturally extended to this scenario, where the spectra of fully concentrated media are regarded as pure spectra and their concentrations in diluted spectra are regarded as abundances. However, the sum-to-one constraint should be relaxed given that a mixed spectrum can be a diluted version of the endmember; e.g., it can be equal to 50% of the endmember. This can be done by simply adding a zero endmember in the LMM, also known as a shadow endmember.

Unmixing pipeline

In remote sensing literature, dierent unsupervised algorithms are presented to extract the endmembers and/or to estimate the abundance matrix. Unsupervised methods, such as Group lasso with unit sum and positivity constraints (GLUP) [START_REF] Ammanouil | [END_REF]] or vertex component analysis (VCA) [Nascimento 2004], can do both in the same time like the previously presented methods PCA and ICA but considering the constraints of equation (4.4). Other unsupervised methods, such as N-FINDR [START_REF] Winter | [END_REF]], only extract the endmembers. It then requires the use of a supervised algorithm to calculate the abundances, again considering the constraints of equation (4.4). A supervised method, called fully constrained least-square (FCLS) [Heinz 2001], has reached interesting performances in remote sensing eld when endmembers are accurately extracted. It has been used together with GLUP and VCA (considering only the extracted endmembers) and N-FINDR. It could also be used with endmembers calculated with other strategies.

Because of the FCLS performances, this method is used to calculate the abundance matrix in this comparison study. The other algorithms (GLUP, VCA and N-FINDR) are compared to nd the most appropriate one to extract endmembers from multispectral photoacoustic dataset.

Endmember extraction methods

4.2.3.1 Group lasso with unit sum and positivity constraints GLUP assumes that the endmembers are unknown but present in the image, among the observations [START_REF] Ammanouil | [END_REF]]. Given this assumption, and without loss of generality, the linear mixing model (equation (4.3)) can be reformulated as follows:

A i = N j=1 u Gji A j + g i , ∀ i ∈ [1; N ] (4.5)

Abundance calculation

Similarly as in section 4.2.2.2, u Gji is the abundance of A j in A i . On one hand, if A j is an endmember, the row U G j of the matrix U G , which is the abundance matrix calculated with GLUP algorithm, has non-zero entries and represents the corresponding abundance map. On the other hand, if A j is a mixed pixel, U G j has all its elements equal to zero. As a consequence, U G admits N -k rows of zero, the other rows being equal to rows of U . The premise in GLUP is that U G allows the identication of the endmembers in A through its non-zero rows. This property is exploited in GLUP in order to nd the endmembers among the observations. The unmixing problem under investigation requires that U G only has a few rows dierent from zero, in addition to the non-negativity and sum-to-one constraints which leads to the following convex optimization problem:

min U G ( 1 2 N j=1 (A -AU G ) 2 F + µ N j=1 U G j 2 ) subject to    u G ji ≥ 0 ∀ i, j N j=1 u Gji = 1 ∀ i (4.6)
with µ ≥ 0 a regularization parameter, and A = [A 1 , . . . , A N ] the dataset matrix. The rst term in equation ( 4.6) ensures that the observations match the model of equation ( 4.5), the second term is the Group Lasso regularization which induces sparsity by possibly driving several rows of U G to zero [Yuan 2006]. The minimization is constrained to ensure that the abundances obey the positivity and the sum-to-one constraints. The resulting optimization is solved using a primal dual method (see [START_REF] Ammanouil | [END_REF]). In conclusion, GLUP allows to identify the endmembers in A by identifying the nonzero rows in U G . Note that GLUP also provides the estimated abundances corresponding to the non-zero rows in the estimated matrix U G . Similarly to [START_REF] Ammanouil | [END_REF]], given the endmembers estimated by GLUP, FCLS is used to estimate the abundances.

Vertex component analysis

VCA also assumes the presence of pure pixels in the data [Nascimento 2004]. The principle of this algorithm is to project data onto a direction orthogonal to the subspace constituted from the endmembers already extracted. The new extracted endmember is the farthest signal in this projection. A new subspace, considering this endmember, is calculated and the same procedure is iteratively done until the extraction of the asked number of endmembers, k.

The rst step of this algorithm is to nd the initial considered subspace. It can be calculated using two dierent methods depending on the dataset SNR. If the SNR is superior to SNR th (equation (4.7)), this rst subspace is calculated using singular value decomposition (SVD) algorithm [Scharf 1991]. Otherwise, the considered subspace is constructed with the k -1 rst axes extracted by PCA. SNR th is dened as follows:

SN R th = 15 + 10 log 10 (k) (4.7)

Let us named the z-th subspace SU B z of the z-th iteration. The vector v z orthonormal to SU B z is calculated as follows:

v z = r z -SU B z SU B + z r z r z -SU B z SU B + z r z (4.8)
Aneline DOLET with r z is a zero-mean random Gaussian vector and SU B + z the pseudoinverse matrix of SU B z . The following equation then describes the projection of the dataset A onto v z .

f z = v T z A (4.9)
The extreme data of f z is found and considered as an endmember. This new endmember is then added to the already extracted ones to create the subspace considered in the next iteration. When the k endmembers are extracted, they are put into their matrix E V . The unmixing matrix containing the abundance coecient U V is calculated projecting the dataset A onto E V . In order to compare each of the method to extract the endmembers, only E V is considered to use FCLS for the abundance estimation.

N-FINDR

N-FINDR is a method which also assumes that a pure pixel for each medium to unmix is present in the dataset [START_REF] Winter | [END_REF], Plaza 2005]. The rst step is the random generation of a set of k endmembers producing then the matrix E N 0 . At each iteration z, the following volume is calculated:

V z = |det(Etemp z )| (k-1)! with E tempz = 1 . . . 1 E Nz 1 . . . E Nz k (4.10)
where E Nz p is the p-th endmember contained in E Nz (which is the matrix of endmembers at the zth iteration). All the N dataset pixels of A (x i , i ∈ [1 ; N ]) iteratively takes place in the E tempz matrix at the position of each E Nz p vector, p ∈ [1 ; k]. The volume V z+1 is again calculated. For example:

V z+1 = det 1 . . . 1 E Nz 1 . . . x i . . . E Nz k (k -1)! (4.11)
If the new calculated volume is greater than the previous one, V z and E tempz are updated with the new values. A pixel x i of the dataset A is then considered as an endmember and inserts in the matrix E tempz as follows:

if V z+1 > V z , V z = V z+1 and E Nz = E Nz 1 . . . x i . . . E Nz k (4.12)
When all the N pixels of the dataset have been tested, the k calculated endmembers are in the nal matrix E N . No abundance coecients are calculated with this method, the FCLS is used to this end.

Spatio-spectral mean-shift

The SSM-S clustering method has been previously presented in Chapter 3 as a multispectral PA clustering method. The spatio-spectral regularization allows the access to the clustering result and to the representative spectrum of each cluster (the average spectrum of each cluster).

70 Aneline DOLET

Results and performances

As already explained in Chapter 3, the spatio-spectral regularization creates more cluster than the number of media to cluster. To reduce the number of clusters, only the k largest labeled areas are kept. The representatives spectra of these k clusters, calculated by the SSM-S, can be considered as the endmembers to use the FCLS. These spectra are got into the matrix E S of the endmembers extracted by the SSM-S.

Abundance estimation

To estimate the abundances knowing now the endmembers, an unmixing method, which demonstrated interesting performances in remote sensing unmixing, called fully constrained least-square (FCLS), is used. This method takes into account the constraints of equation (4.4) which are required in our context.

FCLS solves the following equation:

min U ( 1 2 N j=1 (A -U E) 2 F ) with    u ji ≥ 0 ∀ i, j N j=1 u ji = 1 ∀ i (4.13)
which is a similar problem to the GLUP one, presented in equation ( 4.6), with µ = 0 and using the estimated endmembers matrix E rather than the dataset matrix A. FCLS is better conditioned than the problem of equation ( 4.6) because the endmembers are now known. It should then lead to well abundance assessments.

Depending on the method used to extract the endmembers, the considered estimated endmembers matrix E could be: E G , E V , E N or E S if GLUP, VCA, N-FINDR or SSM-S are used, respectively. It generates then the matrix U G , U V , U N or U S if GLUP, VCA, N-FINDR or SSM-S are used, respectively.

Results and performances

Method of validation

All the unmixing methods presented above calculate abundance maps corresponding to each medium to unmix in the considering dataset. To measure the performances of the dierent methods, these abundance maps are compared to the ground truth.

To do so, regarding a single abundance map, corresponding to an endmember, an average concentration value is calculated on each part of the dataset (using the mask B and its labeling, presented in the previous Chapter). For example, on a phantom composed of a pure medium and a dilution of it, two average values are calculated: one on the pure medium part and one on the dilution region. As we assume that the pure medium, corresponding to the endmember of interest, is present in the imaged region, these average values are normalized by the maximum. If the unmixing is well done, the maximum value corresponds to the pure medium of the considered endmember and is then set to 1. The other normalized values correspond then to the concentration relative to the pure medium.

This calculation is done on all the abundance maps given by the dierent unmixing methods. The normalized average values are compared to the ground truth to measure the performances of the dierent algorithms.

Aneline DOLET

On medium dilution

Blood-Ink(PAT)

First, the results of the unmixing methods of the literature in multispectral photoacoustic imaging (PCA and ICA) are presented for the Blood-Ink(PAT) dataset (Figure 4.1 and Table 4.1). The Figure highlights the abundance maps of PCA (a) and ICA (b). The described ground truth can be seen in Figure 4.2, rst column. Table 4.1 gives the results of PCA and ICA using the validation method presented before, as well as the ground truth. The result did not provide a satisfactory unmixing as the pure ink part is not accurately localized. Indeed, on the ink abundance map, the maximum average value is found in the pure blood inclusion (highlighted in blue in the Table ). Also for ICA, the maximum average value in the blood abundance map is found in the diluted blood inclusion (highlighted in green in the Table ). The PCA and ICA abundance maps contained negative values that do not t a possible concentration value. This is mostly the reason why the results cannot be exploited to access medium concentrations. Indeed, the uncorrelated principal components, and the independent components, could not be considered as an accurate representation of pure media.

The impossible use of PCA and ICA algorithms to access accurate medium concentrations motivates the use of FCLS which is under positivity and sum-to-one constraints. Indeed, these constraints allow the access to medium concentrations and can be exploited to their assessment. For the other presented dataset, PCA and ICA results are no more presented. Only the comparison on GLUP, VCA, N-FINDR and SSM-S combined with FCLS is now done. . The shadow endmember map highlights that the diluted blood inclusion (left-up) is the one of unpure medium as it is mostly present in this map. Regarding both other maps, it can be linked to the blood medium and then considered as a dilution of it because it is also really present in the blood endmember abundance map. For all these methods, the pure media are well localized but the calculated average concentrations are higher than the ground truth. The closest, to the ground truth, calculated dilution is reached with VCA algorithm and has a value of 0.75 instead of 0.5. Even if this value is the best one calculated on this dataset, it highlights too low performance. The concentrations that should be close to 0 are superior to 0.2 which is also too far from the ground truth to be considered as correct. Aneline DOLET Aneline DOLET The abundance maps are presented in Figure 4.10 and the average concentration values are summarized Table 4.4. The blue abundance maps calculated with the endmembers extracted with GLUP and VCA highlight the pure blue medium part (left) but the pure green part is dicult to localize. Indeed, the abundance values are similar in central and right parts. However, on the N-FINDR and SSM-S blue abundance maps the three parts can be distinguished. On the SSM-S abundance map, the green part of the phantom has really low concentration which is close to the ground truth as it can also be seen in Table 4.4 with an average value of 0.1. For the green abundance maps, all the methods allow the discrimination of the blue part (left) from the mix part (central). VCA abundance map does not permit the discrimination of the mix (central) and green (right) parts but the discrimination is possible with the other methods.

Results and performances

Regarding the average concentration values of Table 4.4, the rst remarkable value is the blue calculated concentration of the mix assess with VCA with the value of 0.42 which is exactly the ground truth value (highlighted in blue in the Table ). However, with VCA algorithm, the green concentration in the mix as well as the value of the pure media in the opposite abundance maps are higher than the ground Aneline DOLET Aneline DOLET 83 Chapter 4. Unmixing of multispectral photoacoutsic dataset truth (0.89 for the green concentration in the mix instead of 0.67 ; 0.27 instead of 0 for the blue concentration in pure green part ; 0.20 instead of 0 for the green concentration in pure blue part). Considering now the result given by the SSM-S/FCLS pipeline (in red in the Table 4.4). The calculated mixture is composed of 0.35 of blue and 0.73 of green instead of 0.42 and 0.67, respectively. It corresponds to an error of less than 7% which is already interesting result to target biological application. Also, the abundance parts where the concentration should be 0 have calculated values lower than 0.1 using SSM-S which is really satisfactory and better than any other compared methods. The pipeline SSM-S/FCLS is here validated for the unmixing of mix of media.

Conclusion

Dierent quantication methods of multispectral photoacoustic dataset have been compared to assess the medium concentrations. PCA and ICA algorithms were already used for another application in the photoacoustic imaging eld [Glatz 2011] but do not allow the calculation of concentrations as their abundance maps contain negative values.

Considering the other compared methods to extract the endmembers (GLUP, VCA, N-FINDR and SSM-S), the results obtained on the Blood-Ink(PAT) dataset are far from the ground truth. The dilutions calculated with GLUP, N-FINDR and SSM-S are all of ≈ 0.9 instead of 0.5. Even the lower average dilution value calculated, with VCA, is still of 0.75 which corresponds to an error of 25%. This asks an important question about the blood dilution ground truth as the measurement method is unknown.

Using the Vevo LAZR, we have demonstrated that SSM-S procedure combined with FCLS algorithm allows the accurate assessment of medium concentration for dilution and mix of media.

Lots of acquisition systems are used in the photoacoustic eld. Meanwhile, the multispectral PA images are largely impacted by the laser energy, the acquisition geometry (optical excitation vs US probe) and the targeted applications. The assessment of medium concentrations is then challenging because of this variability. The regions of interest of the imaged area rst need to be extracted. To this end, we suggest to use a threshold calculated from Sobel ltering to discriminate between the background and the ROI. The dierent media present in the ROI can accurately be discriminated using the SSM-S procedure to cluster the photoacoustic data. The proposed SSM-S algorithm also extracts the endmembers of the dataset, corresponding to the pure media, which can be used with the FCLS to correctly assess the medium concentrations of the imaged area.

In the present PhD, this pipeline has been proposed and validated as a robust procedure to unmix the multispectral photoacoustic data in order to assess the imaged media concentrations. The sum-to-one constraint imposed in the FCLS can however be discussed. Indeed, it is a strong assumption particularly when diluted media are imaged. To relax this constraint, a shadow endmember has been introduced in this study but other choices could have been done. For example, other unmixing algorithms are developed in the remote sensing eld to deal dierently with it. Our proposed strategy has anyway permitted the accurate calculation of dilution medium and mix of media concentrations, with an error lower than 7% which is acceptable for biological applications and will be tested in the next Chapter.

Introduction

The proposed strategy, based on both clustering and unmixing methods (SSM-S and FCLS), is used on biological tissues to calculate the oxygen saturation rate (sO 2 ). This measurement is of great interest for various medical applications like the follow-up of tumors or the evaluation of tissue aging [Su 2012]. This Chapter presents rstly, preliminary results on tumor tissues imaged with the Vevo LAZR and secondly, a study on the evaluation of tissue death, conducted on the experimental platform of CREATIS. These studies allow rst tests of the SSM-S and FCLS algorithms on in vivo and ex vivo biological tissues dataset.

For in vivo imaging, the sO 2 values calculated on tumors with our methods are compared to the values displayed by the Vevo LAZR. Indeed, such equipment provides a real-time measurement of the sO 2 concentration map that can be stored after the acquisitions. For the ex vivo acquisition, the evaluation of tissue death is studied on dierent types of meat. However, no ground truth is known in such case.

The photoacoustic signal and sO 2 changes over time are analysed with acquisitions at dierent times on the experimental platform.

Aneline DOLET Chapter 5. Application: evaluation of blood oxygen concentration 5.2 sO 2 calculation on tumors

Data acquisitions

The acquisitions have been done with the Vevo LAZR system on one MMTV-NeuT balb/c mouse.

This type of transgenic mouse spontaneously developed breast tumors which can be seen with the PA system after 17 weeks (between 5mm 3 and 10mm 3 size). The tumors can be seen with the naked eyes after around two more weeks. The studied mouse has been developed for another study conducted by the Dr PETRILLI (Centre de Recherche en Cancérologie de Lyon, France) who lends us the mouse for our acquisitions. These acquisitions have been done when the mouse had 20 weeks (Figure 5.1, the white arrows show some tumors). The mouse ROI were shaved the day before with commercial hair remover cream (VEET, Cream hair remover, Reckitt Benckiser, UK) to avoid any imaging noise coming from the mouse hairs. The mouse was anesthetized with a mix of around 3.5% of isourane and 96.5% of oxygen before the PA acquisitions. During the imaging, the mouse was kept asleep inhaling the same mix but with only around 1.5% of isourane inside the white box (Figure 5.1(a)) to minimize the movements between acquisitions at dierent wavelengths. The heartbeat and breathing are controlled by several detectors of the system at the mouse paws (highlighted by red arrows on Figure 5.1) and inside the white box, respectively. It allows the follow-up of the mouse health and monitors the acquisitions in order to acquire images at the same position. Heartbeat and breathing motions are so avoided as much as possible.

The dataset used to test our strategy is composed of 15 images acquired from 680nm to 960nm, 20nm steps. The results are compared to the Vevo LAZR sO 2 map. To this end, acquisitions on the same tumors are also done using the Vevo LAZR Oxy-Hemo mode.

Discussion

The sO 2 map of Vevo LAZR is smoother in comparison to our own results. It could come from dierent reasons and settings. For example, an average lter applied on the concentration maps or directly on the data before the calculation. Anyway, regarding the similarity of both endmembers extracted by the SSM-S (Figure 5.2(b)), the unmixing results on both tumors are not surprising. Considering as endmembers the average spectra of both larger SSM-S clusters is probably not a good strategy in this context. Indeed, in tumor tissues, pure Hb and HbO 2 are not the most present media. As it can be seen on the Vevo LAZR display, the calculated sO 2 values are principally between 25% and 75% which corresponds to mix of both pure media. The endmembers actually extracted represent then mix of media instead of pure media. Assuming that pure Hb and HbO 2 are not present in numerous pixels, dierent strategy can be tested to extract endmembers according to the present data.

The SSM-S can be used to extract more than 3 clusters (tested with 10) and the endmembers can be selected among these clusters (Figures 5.5 The high errors could come from (1) the bad endmember extractions or/and (2) the sum-to-one constraint of FCLS algorithm which is maybe a strong hypothesis in biological tissues, principally for medium dilutions. Both of these questions have been tested. First, to test the endmember extraction, the FCLS has been used with the theoretical spectra of Hb and HbO 2 . The considered endmembers and the calculated sO 2 and error maps are shown (Figures 5.5(b) and 5.6(b)) for both tumors. The average error values are clearly reduced with this method, principally for the rst presented tumor (18.3% instead of 39.17% in the rst tested strategy, see Table 5.1 for others values). This means that the endmember extraction needs to be improved for biological tissue case. Second, the FCLS algorithm has been tested

with endmembers linked to the Vevo LAZR to test the sum-to-one hypothesis in biological tissues case.

To this end, the endmembers have been extracted from the data regarding the sO 2 map calculated using the Vevo LAZR concentration maps (Figures 5.3(b) and 5.4(b)). The maximal and minimal sO 2 values have been detected in this Vevo LAZR map and the corresponding data spectra have been extracted and considered as HbO 2 and Hb endmembers, respectively. The results are again improved with this strategy (Figure 5.5(c) and 5.6(c) and Table 5.1) which is not surprising as we used the considered ground truth to nd the endmembers. However, the average errors are still around 16%. The sum-to-one constraint should then probably be relaxed using another strategy than shadow endmember in the present case.

90

Aneline DOLET 5.2. sO 2 calculation on tumors In this study, the proposed method for sO 2 calculation has been compared to the sO 2 values calculated by the Vevo LAZR Oxy-Hemo mode on two dierent tumors. As the results were signicantly far from the Vevo LAZR ground truth, we investigate some other strategies testing the potential obstacles. The endmember extraction and sum-to-one constraint relaxation needs to be deeply studied but the primarily results on biological tissues of our strategy (considering the presence of pure Hb and HbO 2 but in small quantity) is still encouraging. Moreover, the ground truth of oxygenation rate in the imaged tumors is not known. It is then dicult to know which method allows the assessment of the more accurate sO 2 values. Also, considering the presence of only Hb and HbO 2 media in the imaged area is maybe a too strong hypothesis. The water and other tissue absorbances should probably be taken into account.

Finally, other wavelength values could be tested to highlight their impact on the results. In the future, it will be interesting to compare our results with the optimal map provided in real-time by the Vevo LAZR.

However, for such comparison, the full processing and post-processing applied in the system have to be known.

Evolution of tissue death

Acquisitions and methods

Photoacoustic acquisitions have been conducted on the experimental platform of CREATIS using the CMUT probe for US detection (presented in section 2.2.2). Seven dierent meats have been studied: duck, chicken, pork, veal, tournedos, heifer and beef. All the meats were bought in a supermarket without knowing the animal death time. The evolution studied is then the sO 2 changes in the meat in ambient air (25 • C) from an initial time. The temperature inside the meat has been carefully checked during all the acquisitions and kept constant to be sure that none evolution could come from temperature changes.

All the acquisitions have been done from 470nm to 690nm with 10nm steps. For each meat, 20 acquisitions have been conducted at dierent times: from initial time to 4h and from 19h to 20h, one every 30 minutes. The sO 2 is calculated for each meat and for each acquisition time. To this end, the FCLS is used to calculate the concentration maps of Hb and HbO 2 . Equation (5.1) is then used for sO 2 calculation. As the optical excitation region is really small on the CREATIS experimental platform, we cannot assume that Hb and HbO 2 are present, as pure media, in the imaged area. The theoretical multispectral evolutions of these media (Figure 5.7) are then considered as endmembers instead of using SSM-S to their extraction. Because this strategy has given satisfactory result on tumor tissue, it seemed to be the best choice in this context.

Results

The results are presented only for two dierent meats but the relative evolutions of the photoacoustic signal and of the sO 2 over time are similar for all the tested meats. Figures 5.8 and 5.9 present the result for the beef and pork meat, respectively. The photoacoustic signal increases over time keeping a shape relatively similar (Figures 5.8(a) and 5.9(a)). Only the acquisitions after 19h have dierent shapes for the wavelengths under 550nm. Comparing both meats, the beef gives higher photoacoustic signal than the pork (maximum values of 4.2a.u. and 3a.u., respectively). This was expected as it is principally blood that gives photoacoustic signal and beef is a bloodier meat compare to pork.

Aneline DOLET signal and the decrease of sO 2 is particularly impressive. Indeed, a rise in photoacoustic signal is clearly linked to a decrease of sO 2 . For example, on pork (Figure 5.9), the sO 2 is not exactly decreasing over time: sO 2 = 78.61% at 19h and sO 2 = 78.64% at 19h30. Meanwhile, the photoacoustic signal at 19h is higher (in purple) than the one at 19h30 (in gray) which is not logical compare to other photoacoustic signal evolutions over time, but it is according to the sO 2 evolution.

Discussion

The photoacoustic signal evolution over time after death has been clearly highlighted in this study as it was already done in [Su 2012]. Here, the photoacoustic signal evolutions and the sO 2 changes over time after tissue death have been linked. However, no sO 2 ground truth exists in this study, only the decrease of relative sO 2 values can be linked to the increase of photoacoustic signal. The sO 2 values interpretation is then complex. The cause of these evolutions is however not explained in this study. Of course, the evolution comes from oxygenation changes in the tissue but it is probably linked to lots of other changes. For example, it should also come from the decrease of water amount inside the meat. The meat becomes drier over time and it impacts the photoacoustic signal as water absorbs the optical energy in the studied wavelength range. Stiening, elasticity changes or other evolutions impact for sure the acquired photoacoustic signal but the contribution of these changes in the photoacoustic signal evolution is not quantied here. Only the meat temperature has been checked during all the procedure so, at least, it should not come from it.

Conclusion

Oxygenation rate (sO 2 ) has been calculated with the proposed quantication pipeline in two dierent cases: (1) tumor oxygenation and (2) evolution of tissues after death. In both studies, the proposed method gives satisfactory preliminary results even if further study needs to be done. Indeed, it allows the calculation of satisfactory sO 2 values compared to the ones given by the Vevo LAZR Oxy-Hemo mode.

Also, satisfactory relative sO 2 values have been calculated from photoacoustic signal acquired on the experimental platform of CREATIS.

However, no sO 2 ground truths were available for both studies which leads to dicult accurate conclusion and only relative values can then be interpreted. It could be of great interest to compare the calculated values to real sO 2 . To this end, comparable studies could be conducted removing some blood during the acquisitions and calculating the sO 2 with blood test. One solution could be to image, on the Vevo LAZR, a known solution of Hb and HbO 2 to evaluate the accuracy of both proposed strategies.

Similar acquisition could also be conducted on the CREATIS platform.

Aneline DOLET

Conclusion and perspectives

Photoacoustic (PA) imaging is a promising medical imaging technique which is of interest for several applications. The access to functional information, combined with the structural information of ultrasound imaging, makes photoacoustic imaging a powerful diagnostic technique. Since this technique is relatively recent, dierent PA processing methods have been developed but few of them are enough robust to be applied to data acquired on various experimental systems.

In this context, the objective of the present thesis was to develop a processing chain for the discrimination and quantication of media using multispectral photoacoustic images. An unsupervised clustering method, called spatio-spectral mean-shift (SSM-S), has been validated to discriminate the PA signal from the imaged media. This algorithm, based on spatial and spectral characteristics, has been shown able to cluster the investigated media. However, this algorithm also extracts the reference spectra of the pure imaged media to allow the use of a supervised unmixing algorithm, called fully constrained least-square (FCLS), to calculate the concentration of each pure medium in all of the pixels. Such algorithm, with constraints linked to the positivity of concentrations and the concentrations sum lower than or equal to 1, outperforms the literature strategies. These developments were tested on dierent acquisition systems: a photoacoustic tomography (PAT) set-up, a photoacoustic computed tomography (PACT) experimental platform and a commercial system. This pipeline (SSM-S/FCLS) has been validated and shown as a powerful quantication method to access imaged media concentrations usable on dierent acquisition systems.

During this PhD, the PACT experimental platform of CREATIS has been improved with the cosupervision laboratory. Two real-time photoacoustic modes have been developed and installed in the platform allowing real-time imaging in 2D and 3D. This improvement makes the acquisitions easier and could allow to image moving media in further studies. For example, within the frame of a recent collaboration with a laboratory in Rennes, France, the follow-up of nanoparticles in the body using photoacoustic imaging can be studied [Ciancone 2018]. Indeed, it could be of interest to release drugs in a specic location in the body. By conning the drug in the nanoparticles, PA imaging could allow following-up their movement. Furthermore, the photoacoustic eect could also be used to break the nanoparticles in order to release the drug in the desired area.

In vivo and ex vivo studies were conducted to validate all the developments of this project. The calculation of the blood oxygenation rate (sO 2 ) in biological tissue, quantifying the oxygenated (HbO 2 ) and deoxygenated (Hb) blood in photoacoustic dataset, has been validated in two dierent cases. First, the performance of the proposed pipeline has been compared to the one of a commercial photoacoustic acquisition system (Vevo LAZR) for the calculation of the sO 2 in mouse breast tumors.

Second, the evolution of biological tissues after death has been studied on the experimental platform of CREATIS. The increase of photoacoustic signal over time has been linked to the decrease of the sO 2 .

Further studies should however be conducted to quantify the importance of oxygenation changes in the photoacoustic signal increase. Indeed, it is also linked to other changes like drying or stiening of tissue and elasticity changes. The contribution of all these changes in the photoacoustic evolutions could be deeper investigated to characterize the tissue evolution after death. The assessment of tissue death time could then be studied with multispectral photoacoustic imaging in further studies.

4. 1

 1 abbreviations meaning). (b) PAT acquisition system highlighting a arc-shape ultrasound probe [Brecht 2009]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.9 The optical excitation and US detector are focalized to obtain a thin resolution [Moore 2016]. 12 1.10 PACT set-up where optical source and US probe are (a) separated [Montilla 2013] or (b) close to each other [Needles 2013]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.11 Thesis objectives and contributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

  Figure bottom. The image axes are in mm. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Theoretical multispectral absorption evolutions of Hb and HbO 2 used as endmembers in unmixing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Study on beef. (a) Photoacoustic spectra for the dierent acquisitions and (b) sO 2 calculated from the photoacoustic signal acquired at each acquisition time. . . . . . . . . . . . 5.9 Study on pork. (a) Photoacoustic spectra for the dierent acquisitions and (b) sO 2 calculated from the photoacoustic signal acquired at each acquisition time. . . . . . . . . . . . .

Y

  time of the US wave from s to s q , i.e., number of media to cluster/unmix W = [λ 1 ; λ L ] Used wavelength range, λ l is one of the used wavelengths with l ∈ [1 ; L]SAbsorbance measured by a spectrophotometer S P A Normalized absorbance measurement comparable to multispectral PA signal A Dataset matrix of size N × L a λ l Discrete photoacoustic amplitude of A at the wavelength λ l a * λ l Normalized value of A at the wavelength λ l A i Photoacoustic spectrum of the i-th pixel, A i ∈ R L s i Spatial position of the i-th pixel, s i ∈ R S x i i-th dataset pixel, i ∈ [1 ; N ] = {y i } i=[1;P ] Set of P labeled pixels Z = {z j } j=[1;N -P ] Set of N -P unlabeled pixels Introduction During the last century, improvements in the medical domain permitted the reduction of the death rate linked to many diseases. This decrease comes from treatment enhancements, better diagnosis and earlier illness detection. New technologies and techniques, principally in the medical imaging eld, have favored these improvements. Non-invasive and non-ionizing imaging techniques are particularly interesting. Transportable systems, which can be moved to the patient, are also preferred. Ultrasound (US) imaging techniques bring all these requirements together and a lot of imaging techniques have been developed around this modality during the last decades. For nearly thirty years, a new imaging modality, close to US imaging, has been developed: photoacoustic (PA) imaging, which combines optical and ultrasound imaging. It adds to the US imaging the access to functional information. This characteristic is of great interest for diagnostic and follow-up of diseases. Photoacoustic imaging requires developing new acquisition systems and image processing algorithms to accurately access functional information. The present thesis takes place in this context and aims at contributing to the photoacoustic imaging eld. Principally, the large number of photoacoustic acquisition systems makes dicult the development of robust image processing algorithms adapted to each one. The main objectives of this thesis are to develop, test and validate a processing chain for photoacoustic dataset in order to measure with high accuracy dierent concentrations of an imaged medium. Validation was made on dataset coming from dierent acquisition systems to validate the robustness of the proposed methods. Also, since real-time imaging is a remarkable characteristic of ultrasound imaging, we have developed a real-time photoacoustic imaging mode on the ULA-OP256 research scanner, during my co-supervision stay in the MSDLab in Florence.

Figure 1 . 1 :

 11 Figure 1.1: Optical and ultrasound imaging resolution over the imaging depth (based on [Dubois 2007]).

Figure 1 . 2 :

 12 Figure 1.2: Scattering phenomenon in US imaging [Zahnd 2007].

Figure 1 . 3 :

 13 Figure 1.3: Photoacoustic principle: (a) optical illumination, (b) thermal expansion of the optical absorbers, (c) ultrasound wave generation, (d) acquisition and reconstruction of a photoacoustic image.

  Ultrasound signals acquired on each US probe transducer need to be beamformed to reconstruct the photoacoustic images (Figure 1.3(d)). These signals, called radiofrequency signals (RF signals) or radiofrequency data (RF data), are typically reconstructed using DAS algorithm (Figure 1.4). A pixel at the spatial position s in the photoacoustic image I is reconstructed by delaying and summing the signals as follows:

Figure 1 . 4 :

 14 Figure 1.4: Reconstruction of the photoacoustic images: DAS algorithm considering the delays (τ q ) and weighting factors (W q ).

4 .

 4 Interest of photoacoustic imagingA given medium can conduct to dierent photoacoustic signals because (1) light attenuation impacts the response of tissue as a function of depth (Figure1.5), (2) various concentrations change the optical absorption and (3) dierent target sizes create ultrasound waves at dierent frequencies. For many studies in photoacoustic imaging, light or ultrasonic attenuations are not taken into account before applying image processing method.

Figure 1

 1 Figure 1.5: Photoacoustic image of an homogeneous absorbing medium. The amplitude of the PA image decreases as a function of depth because of light attenuation. The image axes are in mm.

Figure 1

 1 Figure 1.6, green vertical line). Then, HbO 2 and Hb exhibit dierent PA signals, which allows their dierentiation.

  2 and Hb can be dierentiated by their PA spectral evolutions in the range [600nm ; 900nm] [Hill 2016], which is much more dicult in the range [250nm ; 550nm] (see Figure 1.6). The range [650nm ; 900nm] is largely used for biological tissue imaging because, in this range, the optical absorption of water is low. As biological tissues are principally composed of water, this wavelength range allows better imaging depth: it is called the diagnostic window (in gray in the Figure 1.6).

Figure 1 . 6 :

 16 Figure 1.6: Absorption coecients of dierent media as a function of the wavelengths. (a.u.; arbitrary units) Green vertical line highlights Hb and HbO 2 optical absorption at 650nm. The diagnostic window is shown in gray.

  Figure 1.7: Resolution over imaging depth of PAT, PAM and PACT [Wang 2016].

Figure 1

 1 Figure 1.8: (a) PAT acquisition system where the imaged sample rotates inside the Nd:YAG laser/US probe system [Liu 2016] (see the List of abbreviations and symbols for the Figure abbreviations meaning). (b) PAT acquisition system highlighting a arc-shape ultrasound probe [Brecht 2009].

Figure 1 . 9 :

 19 Figure 1.9: The optical excitation and US detector are focalized to obtain a thin resolution [Moore 2016].

Figure 1 .

 1 Figure 1.10: PACT set-up where optical source and US probe are (a) separated [Montilla 2013] or (b) close to each other [Needles 2013].

Figure 1 .

 1 Figure 1.11: Thesis objectives and contributions.

  the entire thesis, thanks to a collaboration with the University of Macau (People's Republic of Macau, China). The acquisition set-up is presented Figure 2.1.

Figure 2 . 1 :

 21 Figure 2.1: Experimental set-up used for the PAT acquisitions.

[

  Vallet 2015], Figure2.2. The multispectral optical excitation is a Nd:YAG pulsed laser (Quanta-Ray INDI Series, Spectra-Physics, USA), with a pulse duration of 6ns and 10Hz repetition rate, parametric oscillator (OPO -versaScan, SpectraPhysics, USA). With this system, it is possible to acquire images with a large excitation wavelength range, from 410nm to 1250nm. There are dierent laser beam outputs that are reported on Figure2.2(b). These outputs give dierent wavelength ranges. Only one output is used at the same time so the others are stopped by dierent beam stops. This laser is the property of CEA in Grenoble, France. With the help of Jean-Marc DINTEN of LISA-CEA-LETI, which has provided a part of this photoacoustic system equipment, it has been possible to realize photoacoustic acquisitions in CREATIS.On this platform, dierent ultrasound scanners can be used: SonixMDP (Ultrasonix Corp., Vancouver, BC, Canada) or ULA-OP64 (Ultrasound Advanced Open Platform -64 transmission/reception channels, MSDLab, Florence, Italy -Figure2.3(a)). As this thesis is in co-supervision with the MSDLab, the use of this last ultrasound scanner has been preferred. For 2D imaging, two dierent US linear array probes are used: LA523E, Acuson -192 elements, 7.8M Hz center frequency, bandwidth 92% at -6dB (Figure 2.3(b)) or CMUT probe, HF3, ACULAB -192 elements, 10M Hz center frequency, bandwidth 100% at -6dB (Figure 2.3(c)). For 3D imaging, a matrix array probe is used: Vermon -8 × 24 elements, 3.84M Hz center frequency, bandwidth 65.1% at -6dB (Figure 2.3(d)).

Figure 2

 2 Figure 2.2: (a) PACT platform at CREATIS. (b) Schematic of this experimental set-up.

Figure 2 . 4 :

 24 Figure 2.4: Linear (a) and orthogonal (b) laser beam/imaged sample/ultrasound probe geometries.

  Figure 2.5(a)). The particularity of this system is that an optical ber is directly inserted inside the US probe, the optical sources are then close to the ultrasound transducers producing an hybrid US/PA probe (Figure 2.5(b)). The PA signals are then acquired by the ultrasound transducers on the same side as the optical excitation (Figure 2.5(c)).The optical source of the Vevo LAZR is a Nd:YAG pulsed laser with pulse duration of 5ns and 20Hz repetition rate coupled with an OPO to access various wavelengths[Arthuis 2017]. The wavelength range can be programmed from 680nm to 970nm. The PA probe is the LZ400, composed of 256 elements, that acquires ultrasounds in the frequency range from 18M Hz to 38M Hz. The imaging possible depth is ≈ 1.5cm and the imaged ROI is of ≈ 1cm large.The images are reconstructed with a DAS algorithm. RF data are not accessible and acquired dataset are impacted by the time gain compensation (TGC) and other post-processing algorithms. The interest of this platform is that ex vivo and in vivo acquisitions on small animals can be done as well as phantom imaging with an investigation depth around 1.5cm.

Figure 2

 2 Figure 2.5: (a) Multispectral photoacoustic commercial system: Vevo LAZR. (b) Hybrid US/PA probe. (c) Sample inside the system with the bimodal US/PA probe ready for acquisitions.

Figure 2

 2 Figure 2.6: Three dierent phantom molds. (a) White mold which is breakable to easily take out the phantom. (b) Two transparent molds with dierent heights.

  system that measures the contributions of both optical absorption and scattering (i.e., the absorbance) is used. Absorbances are measured in the wavelength range [400nm ; 1200nm], 1nm steps, with a Perkin Elmer Lambda900 spectrophotometer (Figure 2.7(a)) at the Institut Lumière Matière (UCB Lyon 1). The maximum absorbance that can be measured by the spectrophotometer made it impossible to perform the acquisitions directly on pure inks. The absorbance measurements are thus performed on diluted inks, with a dilution factor of 1:1601 (Figure 2.7(b)). These measurements have been done with the help of Elodie ROMEO and Thomas DEHOUX from the Institut Lumière Matière.

Figure 2

 2 Figure 2.7: (a) View of inside the spectrophotometer (diluted brown ink and water as reference in test tubes for the acquisition). (b) Diluted inks for the optical spectrophotometry measurements.

Figure 2

 2 Figure 2.8: (a) Blue phantom in the set-up. (b) Blue phantom during the acquisition of the multispectral photoacoustic signals.

Figure 2

 2 Figure 2.9: (a) Normalization ratio for spectrophotometry spectra measured with black ink. (b) PA spectra (A, * marker) and normalized absorbance spectra (S P A , solid line) for (b) brown, (c) yellow, (d) red, (e) blue and (f) green inks.

Figure 2 .

 2 Figure 2.10: Three bi-colored phantoms: (a) red/green, (b) blue/red and (c) black/red.

  Figure 2.11(d)). The wire used for the production can stay inside the phantom or be removed, experimental tests show that it gives negligible PA signal dierences.

Figure 2 .

 2 Figure 2.11: (a) Dierent sizes of molds for spherical inclusions, red arrows show the holes for lling. (b) Needle used to ll the inclusions. (c) Inclusions xed inside a phantom mold. (d) Uncolored phantom with colored inclusions (green and red).

Figure 2 . 4 .

 24 Figure 2.12: (a) Rectangular inclusion mold. (b) Rectangular inclusions of dierent colors and sizes. (c) Inclusions xed inside a phantom mold.

Figure 2 .

 2 Figure 2.13: ULA-OP256 research scanner.

Figure 2 .

 2 Figure 2.14: (a) Active and (b) passive ultrasound time of ight(s).

4 .

 4 photoacoustic modes are used with the same card, the dierent image reconstruction times in 2D or 3D just need to be taken into account.

Figure 2 .

 2 Figure 2.15: Scheme of the ULA-OP256 triggered by the laser pulse. The yellow arrows represent the trigger for passive acquisitions.

Figure 2 .

 2 Figure 2.16: External card used to make the trigger between the laser and the ULA-OP256 scanner.

Figure 2 .

 2 Figure 2.17: (a) Cross-sectional phantom image: the yellow rectangle highlights the region imaged with the active mode and the orange ellipse, the laser beam size. (b) Photoacoustic set-up. The visible green light is due to a residual laser component at 532nm. (c) Screenshots of the real-time display for dierent phantom positions. The photoacoustic signal is in blue.

Figure 2 .

 2 Figure 2.18(b)-right and the 1mm diameter circle plan Figure 2.18(b)-left. The photoacoustic signal is in blue and the US images are in black and white. As the 2D imaging validation, the 1064nm wavelength is used.

Figure 2 .

 2 Figure 2.18: 3D US/PA imaging. (a) Cylindrical gold particle (1mm diameter and 3mm length). (b) Screenshots of the real-time display at the center position. The photoacoustic signal is in blue. (c) The three central slice of the 3D active (US) signal with a 6dB dynamic. (d) The three central slice of the 3D passive (PA) signal with a 6dB dynamic. The resolution is given in each plane and is consistent with the particle dimension. The image axes are in mm.

Figure 2 .

 2 Figure 2.19: (a) Blue / 0.53 blue dilution / green colored phantom. (b) Spectra of each region of interest. (c) B-Bdil-G(Vevo) dataset. The image axes are in mm.

Figure 2 .

 2 Figure 2.21: (a) Blood and ink inclusions cylindrical phantom, top view. (b) Spectra of each region of interest. (c) Blood-Ink(PAT) dataset. The image axes are in mm.

Figure 2 .

 2 Figure 2.22: (a) Small phantom blue, green and both mix (0.74 blue and 0.29 green). (b) Spectra of each inclusion part. (c) IncB-mix-G(Vevo) dataset. The image axes are in mm.

Figure 2 .

 2 Figure 2.23: (a) Uncolored 4% agar phantom with blue and green 5mm diameter spherical inclusions. (b) Spectra of each inclusions. (c-e) Three dierent slices of the 3D(Vevo) dataset for ve wavelengths. The image axes are in mm.

For a given medium

  which should be clustered in a single label, the proposed methods have to face up to three main diculties: (1) the concentration within the medium can change spatially, (2) the light is attenuated as a function of depth and (3) US frequency and amplitude can change as a function of the optical absorber dimension. Consequently, media discrimination in multispectral photoacoustic imaging is a challenging task which requires robust clustering methods. A part of this Chapter has been published in EURASIP Journal on Advances in Signal between background and regions of interest When a region is imaged in PA imaging, some areas give photoacoustic signal, those with optical absorbers, and others do not. Multispectral PA images are then composed of ROI with multispectral PA signal (in purple Figure 3.1), and pixels containing only noise, called background, in green Figure 3.1.

Figure 3 . 1 :

 31 Figure 3.1: Representation of multispectral photoacoustic dataset. (a) Signals from the imaged region are acquired at dierent wavelengths with two pixels identied (in green and purple). (b) Spectra of two dierent pixels: in green, pixel from the background with low PA amplitudes at all wavelengths and in purple, pixel from the ROI with signicant PA amplitudes.

  background pixels of the matrix C SF are labeled to zero where the other pixels are labeled in the range [1 ; k].

2

  with A and E c the mean spectrum of the dataset and endmember matrix, respectively.All the intra-class correlation coecients are placed into the matrix ICC of size N × k. The clustering result matrix C ICC of size N × 1 corresponds to the endmember label e i which is the maximal icc ie coecient of the i-th pixel, m i , explained by the following equation:

Figure 3 . 2 :

 32 Figure 3.2: The spatio-spectral mean-shift principle. (a) Pixel spatial features at the rst iteration with the R S parameter, and (b) spectra of three pixels with the R λ parameter. The red pixel is the reference pixel.Although the green and blue pixels are both in the spatial area of the reference pixel, only the green pixel is used to update the features of the reference. Indeed, the spectrum of the blue pixel is not close enough to the one of the reference sample (by at least one wavelength; the blue and red spectra are separated by a distance greater than R λ ).

  Figure 3.4(f) because the normalization aims at merging both blood inclusions altogether which means that only two media (ink and blood) have to be discriminated.

Figure 3 . 4 :

 34 Figure 3.4: Pre-processings and endmembers calculation of Blood-Ink(PAT) dataset. (a) Photoacoustic image at 730nm. (b) Noise removed by the pre-processing threshold. (c) Binary mask B, the pixels of interest are in black. (d) Endmembers. (e) Labeling of the binary mask B which corresponds to the ground truth when no normalization is done. (f) Normalized endmembers. (g) Labeling of the binary mask B when the normalization step is applied. (h) Normalized data. The image axes are in mm.

Figure 3

 3 Figure 3.6: Blood-Ink(PAT) dataset SSM-S method without normalization. (a) Result of SSM-S procedure with numerous clusters. (b) The k = 3 biggest clusters are kept (ink in black, pure ink in red and diluted blood in blue), the unlabeled pixels are in gray. (c) SSM-S nal result, the purple arrow shows the misclustered pixels at the inclusion boundary. The image axes are in mm.

  Figure 3.7: Blood-Ink(PAT) dataset clustering method results with normalization. First two lines: pure blood and ink abundance maps. Last line: clustering map results (blood in red and ink in black) with W CP performance values. Columns, from left to right: SF, LS and ICC method results. The image axes are in mm.

Figure 3

 3 Figure 3.8: Blood-Ink(PAT) dataset SSM-S method with normalization. (a) Result of SSM-S procedure which shows the reduced number of clusters, as both media clusters are already well done. (b) The k = 2 biggest clusters are kept (ink in black and blood in red), the unlabeled pixels are in gray at the bottom inclusion boundary (highlighted by the yellow arrow). (c) SSM-S result after post-processing, purple arrows show the misclustered pixels at the blood inclusion boundaries. The image axes are in mm.

  Figure 3.9: B-Bdil-G(Vevo) dataset results without (a-c) and with (b-d) normalization, respectively. (a) Endmembers used as references for the supervised methods when no normalization is applied. (b) Endmembers calculated when the normalization step is processed. (c) Method result maps without normalization, from up to bottom: ground truth, SF, LS, ICC and SSM-S results (the diluted blue cluster is represented in cyan). The performances are also highlighted next to the results. (d) Same as (c) but with the application of the normalization step. The image axes are in mm.

  Figure 3.10: B-mix-G(Vevo) dataset results without (a-c) and with (b-d) normalization, respectively. Endmembers used as references for the supervised methods when (a) no normalization or (b) the normalization step is processed. (c) Method result maps without normalization, from up to bottom: ground truth, SF, LS, ICC and SSM-S results (the mix is represented in orange). The performances are also highlighted next to the result maps. (d) Same as (c) but with the application of the normalization step. The image axes are in mm.

Figure 3 .

 3 Figure 3.11: SSM-S steps with B-mix-G(Vevo) dataset and with the normalization applied. (a) SSM-S procedure result with numerous clusters. (b) Unlabeling step. Only the k = 2 biggest clusters are kept. The unlabeled pixels are in gray. The image axes are in mm.

4 :

 4 Figure 3.12: IncB-mix-G(Vevo) dataset results without (a-c) and with (b-d) normalization, respectively. Endmembers used as references for the supervised methods when (a) no normalization or (b) the normalization step is processed. (c) Method result maps without normalization, from up to bottom: ground truth, SF, LS, ICC and SSM-S results (the mix is in orange). The performances are also highlighted next to the result maps. (d) Same as (c) but with the application of the normalization step. The image axes are in mm.

Figure 3 .

 3 Figure 3.15: Setting of R S and R λ . Success map for varying R S and R λ with the Blood-Ink(PAT) dataset, when the normalization step is included. (Scale in %)

Figure 4 .

 4 Figure 4.1: PCA and ICA results on the Blood-Ink(PAT) dataset. (a) PCA and (b) ICA algorithms abundance maps. The image axes are in mm.

Figure 4 .

 4 Figure 4.2 presents the results of the unmixing methods coming from remote sensing elds (GLUP, VCA and N-FINDR), as well as the proposed SSM-S. All these methods are used to extract the endmembers from the dataset. The FCLS is after used to calculate the abundance maps and nally measure the average concentration in each ROI. All the endmembers extraction are rst presented to highlight where the pure media are found depending on the considered method.

  Blood-Ink(PAT) dataset unmixing results summary for PCA and ICA methods.Considering all the pixels of interest to apply the GLUP requires too much computation time. In this way, 200 pixels are randomly chosen inside the ROI and the endmember extraction is done only considering them. These pixels, as well as their spectra, are presented in, respectively.The GLUP extracted endmembers are chosen in both pure media (Figure 4.3(c), red arrow shows the position of the pure blood extracted endmember and dark green arrow the one of the pure ink endmember). The corresponding spectral evolutions are presented in Figure 4.3(d) as well as the added shadow endmember (only composed of zero values). The abundance maps calculated by the FCLS with these three endmembers are highlighted Figure 4.2, second column. The rst two lines validate the well unmixing of both pure inclusions (right-up for the pure ink one and the bottom one for pure blood) as they are entirely present in the corresponding abundance map with high concentration values (and low concentration values in the other abundance map)

Figure 4 . 2 :

 42 Figure 4.2: Remote sensing unmixing method results on the Blood-Ink(PAT) dataset. The rst lineshows the abundance maps corresponding to the blood endmember, the second to the ink endmember and third the abundance of the shadow endmember. Each column gives the results of a method combined with FCLS: ground truth, GLUP, VCA, N-FINDR and SSM-S (from left to right). The image axes are in mm.

Figure 4 .

 4 Figure 4.3: GLUP endmember extraction. (a) 200 randomly chosen pixels. (b) 200 spectra of the randomly chosen pixels. (c) Both extracted endmembers positions. (d) Both extracted endmembers as well as the added shadow endmember. The image axes are in mm.

Figure 4 . 4 :

 44 Figure 4.4: VCA endmember extraction. (a) Both extracted endmembers positions that are at the boundaries of the pure blood inclusion. The ink extracted endmember is then not really of the pure medium. (b) Both extracted endmembers as well as the added shadow endmember. The image axes are in mm.

Figure 4 .

 4 Figure 4.5: N-FINDR endmember extraction. (a) Both extracted endmembers positions that are in the same inclusion. The ink extracted endmember does not come from the pure medium. (b) Both extracted endmembers as well as the added shadow endmember. The image axes are in mm.

Figure 4 .

 4 Figure 4.6: SSM-S endmember extraction. (a) Class kept during the post-processing. The average spectra of each class corresponding to both pure media (in dark green for ink and red for blood) are the extracted endmembers. (b) Both extracted endmembers as well as the added shadow endmember. The image axes are in mm.

Figure 4 . 8 :

 48 Figure 4.8: Remote sensing unmixing method results on the B-Bdil-G(Vevo) dataset. The rst column shows the abundance maps corresponding to the blue endmember, the second to the green one. Each line gives the results of a method combined with FCLS: ground truth, GLUP, VCA, N-FINDR and SSM-S, from up to bottom. The image axes are in mm.

Figure 4 . 9 :

 49 Figure 4.9: Endmembers extracted with (a) GLUP, (b) VCA, (c) N-FINDR and (d) SSM-S.

Figure 4 .

 4 Figure 4.10: Remote sensing unmixing method results on the B-mix-G(Vevo) dataset. The rst column shows the abundance maps corresponding to the blue endmember, the second to the green one. Each line gives the results of a method combined with FCLS: ground truth, GLUP, VCA, N-FINDR and SSM-S, from up to bottom. The image axes are in mm.

  4: B-mix-G(Vevo) dataset unmixing results summary for GLUP, VCA, N-FINDR and SSM-S combined with the FCLS algorithm.

Figure 5

 5 Figure 5.1: Mouse with breast tumors in the Vevo LAZR system. White and red arrows highlight tumors and heartbeat detectors, respectively.

Figure 5

 5 Figure 5.2: SSM-S procedure to extract the endmembers. (a) 3 main clusters kept after the SSM-S procedure (red, yellow and white). The black pixels are the unlabeled ones. (b) Extracted endmembers corresponding to HbO 2 (yellow pixels in (a)) and Hb (white pixels in (a)). The image axes are in mm.

Figure 5

 5 Figure 5.3: Hb and HbO 2 concentration maps and sO 2 map (in %) calculated with (a) SSM-S, FCLS and equation (5.1) and (b) Vevo LAZR concentration maps and equation (5.1), from left to right, respectively. (c) Vevo LAZR display: left, US image of the tumor and right, sO 2 map. (d) sO 2 calculation errors (in %) between our method and the Vevo LAZR one. The image axes are in mm.

  (a) and 5.6(a)). The selection is done by comparison with the theoretical absorbance spectra of Hb and HbO 2 . The closest extracted endmember, to the theoretical spectrum, is considered for the FCLS calculation for both media. The sO 2 map calculated with FCLS and the resulting errors map are presented on the same Figure (center and right, respectively) as well as the average error value. On Figure5.6, the average error value calculated on the reduce ROI is highlighted in green. With such strategy, the average error values are reduced (see

Figure 5 . 5 :

 55 Figure 5.5: Test of new strategies for sO 2 calculation with FCLS on the rst presented tumor: (a) considering pure Hb and HbO 2 present in the data but in small quantity, (b) using the theoretical absorbance spectra of Hb and HbO 2 and (c) extracting the endmembers from the Vevo LAZR sO 2 map which is considered as the ground truth and presented at the Figure bottom. The image axes are in mm.

Figure 5

 5 Figure 5.7: Theoretical multispectral absorption evolutions of Hb and HbO 2 used as endmembers in unmixing.

  Figure 5.8: Study on beef. (a) Photoacoustic spectra for the dierent acquisitions and (b) sO 2 calculated from the photoacoustic signal acquired at each acquisition time.

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Combination of optical and ultrasound imaging . . . . . . . . . . . . . . . . . . . 4 1.2.1 Optical imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.2 Ultrasound imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Ultrasound wave generation with pulsed light illumination . . . . . . . . . . . . 6

	Chapter 1
	Photoacoustic imaging
	Contents
	1.1

1.3.1 Physical concept of the photoacoustic eect . . . . . . . . . . . . . . . . . . . . . . . 6 1.3.2 Mathematical representations of photoacoustic signal . . . . . . . . . . . . . . . . . . 6 1.3.3 Photoacoustic image reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Interest of photoacoustic imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4.1 Photoacoustic characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4.2 Interest of multispectral illumination in photoacoustic imaging . . . . . . . . . . . . 9 1.5 Current researches in multispectral photoacoustic imaging . . . . . . . . . . . . 10 1.5.1 Photoacoustic imaging systems and probes . . . . . . . . . . . . . . . . . . . . . . . 10 1.5.2 Examples of biological and medical multispectral photoacoustic applications . . . . . 12 1.5.3 Image processing methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.6 Motivations of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.6.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.6.2 Targeted applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.7 Manuscript organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

  State of the artMulti-colored phantoms are produced to do preliminary tests but more complex phantoms are required to go further. The objective is to put small optical absorbing objects in an uncolored 4% agar phantom to acquire dataset with punctual photoacoustic signals. These objects are called inclusions and are also made with colored 4% agar.

	2.3. Phantom manufacturing
	2.3.4 Inclusion production
	2.3.4.1

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Pre-processings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

	Chapter 3
	Clustering of multispectral photoacoustic
	dataset
	Contents
	3.1

3.2.1 Discrimination between background and regions of interest . . . . . . . . . . . . . . 42 3.2.2 Two strategies to deal with concentration dierences . . . . . . . . . . . . . . . . . . 43 3.3 Methods of the literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3.1 Endmembers calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3.2 Spectral-tting method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3.3 Least-square minimization algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.4 Intra-class correlation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4 Spatio-spectral mean-shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.1 Rationale of the study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.2 Spatio-spectral regularization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.3 Cluster number reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.5 Results and performances . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5.1 Quantitative validation method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5.2 Medium dilution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5.3 Mix of media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.5.4 3D dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5.5 Setting of R S and R λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

Table 3

 3

			(pixels) (pixels) (a.u.)		ICC SSM-S
	without with	3 2	8 1	15	19	76 71 85 87	78 87	99 98

.1: Blood-Ink(PAT) dataset summary of the used parameter values and the clustering method performances. The best results are highlighted in red.

Table 3

 3

	Ink(PAT)	without with	3 2	8 1	15	19	76 71 85 87	78 87	99 98
	B-Bdil-G(Vevo)	without with	3 2	3 1	30	2 4	65 63 71 70	54 70	86 90
	B-mix-G(Vevo)	without with	3 2	30 1	30	10	69 62 93 92 93 89	94 38
	IncB-mix-G(Vevo)	without with without with	3 4	3 1 1	40	2 10 2 10	51 51 52 59 57 61 55 68	54 59 58 68	71 68 97 95

.5: Clustering parameters and results summary. The best performance for each dataset and strategy are highlighted in red.

Table 4

 4 .2 summarizes these concentration measurements.Compared with PCA and ICA calculation, both pure media are well localized. However, the calculated blood dilution is of 0.91 for a ground truth of 0.5. The concentration calculated on the ROI of pure media, in the other pure medium abundance map, are of ≈ 0.86 for both instead of 0. Even if the localization of pure medium regions of interest is satisfactory, the calculated concentrations are high compared to the ground truth.

	Aneline DOLET

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.2 sO 2 calculation on tumors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2.1 Data acquisitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.2.2 sO 2 calculation pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

	Chapter 5
	Application: evaluation of blood oxygen
	concentration
	Contents
	5.1

5.3 Evolution of tissue death . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.1 Acquisitions and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

Table 5 .

 5 Table5.1 which summarized all these values) but still too high to consider the result exploitable for tumor follow-up. 1: Summary of the average sO 2 error values.

	Endmember extraction	First presented tumor	Second presented tumor All the ROI Reduced ROI
	3 main SSM-S clusters	39.17%	37.09%	23.42%
	10 main SSM-S clusters	31.76%	32.20%	19.93%
	Theoretical spectra	18.30%	26.69%	16.84%
	From Vevo LAZR sO 2 map	17.40%	23.48%	15.04%

Aneline DOLET

Acknowledgments

I would like to begin these acknowledgments thanking my jury. I am greatfull to my reviewers, who have accepted to read and evaluate this manuscript. I also want to thank all the jury for their presences

3D dataset

The 3D(Vevo) dataset is processed only with the SSM-S method. The objective is to highlight that the parameters R S and R λ have to be tuned according to the imaged media and not the region of interest geometry. Indeed, it could be thought that the R S parameter really inuences the result while the rst feature which has to be taken into account is the spectral one. The dataset used here is composed of two pure media (k = 2).

Figure 3.14 shows the result maps with the SSM-S tuned with R S = 30pixels and R λ = 10a.u.

The algorithm is used only in 2D. The data is processed slice by slice (Figure 3.14(a-b-c)) and then reconstructed in 3D ). As the regions of interest are spherical inclusions, the algorithm is tested with dierent distances between both pure media using the same R S parameter. The unmixing methods coming from the remote sensing eld as well as the SSM-S are now compared on the B-Bdil-G(Vevo) dataset. As the extraction method steps have been largely presented on the Blood-Ink(PAT) dataset, the focus is now only on the results: the extracted endmembers and the abundance maps.

The endmembers extracted with GLUP, VCA, N-FINDR and SSM-S are presented in Figure 4.7(ad), respectively. All the extracted endmembers are dierent and allows the discrimination of both pure media (here, blue and green 4% agar). The added shadow endmembers are also present in the Figure. 4.3. For all the methods, the pure media are well localized and the green abundance maps clearly highlight the green part. The best one is assessed with SSM-S as really low green abundances are calculated for both blue phantom parts. Indeed, the corresponding average concentrations are 0.05 and 0.07 for blue and diluted blue parts, respectively. The blue abundance maps have, on the contrary, higher values for the green part. However, the SSM-S blue abundance map have the lowest value (0.26).

Meanwhile, SSM-S allows the calculation of a dilution value of 0.51 instead of 0.53 for the ground truth.

For this dataset, the unmixing pipeline SSM-S/FCLS reached the best performances (highlighted in red in the Table 4.3) and allows the access to an accurate blue dilution factor. 5.2.2 sO 2 calculation pipeline

On Vevo LAZR system

The Oxy-Hemo mode acquires two photoacoustic images at 750nm and 850nm. With both images, the concentrations of Hb and HbO 2 (expressed by [Hb] and [HbO 2 ], respectively) are measured [START_REF] Visualsonics | Fujilm Visualsonics. Oxy-Hemo Mode[END_REF]]. The saturation in oxygen is then calculated with the expression:

(5.1) This is the classically used ratio in biomedical applications to access oxygenation information [Li 2009, Deán-Ben 2014]. Both concentration maps are saved, as well as the Vevo LAZR sO 2 displayed map for further comparisons.

SSM-S unmixing strategy

The pre-processing threshold (section 3.2.1) is rst applied on the dataset. The SSM-S is after processed to extract the endmembers corresponding to Hb and HbO 2 media. Finally, the abundance maps of both media are calculated with the FCLS algorithm. The sO 2 of each pixel is assessed using the abundance maps and the equation (5.1).

Results

The abundance maps, corresponding to the concentrations assessed, extracted with our strategy can be compared to the ones saved with the Vevo LAZR. The sO 2 maps are calculated and can be compared Aneline DOLET

sO 2 calculation on tumors

To measure quantitatively our method performances, errors between the two sO 2 maps for which sO 2 values are accessible (Figure 5.3(a-b), right) are computed. To this end, the absolute dierence between both measurements is calculated for each pixel. The error is then in % and its map is presented in Figure 5.3(d). The ROI average error is of 39.17% which is high. Some pixels have really good sO 2 agreement and so low error but others have errors superior to 50% which is not convincing concentration assessment.

The validation has been done on another tumor (Figure 5.4). The same comments can be done and the average error value on all the ROI is of 37.09% which is again too high to be exploited (Figure 5.4(d)).

Here, the ROI contains lots of pixel from deeper tissues that do not correspond to the tumor. The average error needs then to be calculated on a smaller region (green rectangle on Aneline DOLET Chapter 5. Application: evaluation of blood oxygen concentration

During this work, various processing methods and acquisition systems have been considered for PA imaging. Some ideas have been excluded in the scope of the PhD but could deserve a deeper study for future work. First, as deep learning algorithms are now used in dierent elds with good performance [Leclerc 2017, Gasse 2017], it would be interesting to test it in the photoacoustic domain. However, a major drawback of such strategy is the database creation. Indeed, numerous data are required to train robust deep learning algorithms and the large variability of acquisition systems could hamper these developments. Deal with the tests of this type of algorithms in depth in photoacoustic elds could anyway be of interest as soon as enough data have been collected with the corresponding ground truths.

3D real-time photoacoustic imaging has also been developed during this PhD. It needs to be further studied. For example, CREATIS and LABTAU collaboration allows the use of 4 Verasonics ultrasound systems altogether to display and acquire high quality 3D images. One objective would be to test this unique platform (3 systems with similar performance in the community) in photoacoustic imaging. The increase of the number of elements of the probes would also lead to a higher contrast and better resolution.

Finally, the 3D real-time mode already developed could also be deeply used to develop, test and validate the proposed unmixing strategy on volume.

Finally, photoacoustic imaging is a powerful technique to access blood functional information. It

National journal

A. Dolet, F. Varray, E. Roméo, T. Dehoux, and D. Vray. Spectrophotometry and photoacoustic imaging: a comparative study. IRBM, vol. 38, no. 6, pages 352-356, 2017. https://doi.org/10.1016/j.irbm.2017.10.003.

International conferences