*uef.fi*

**PUBLICATIONS OF **

**THE UNIVERSITY OF EASTERN FINLAND**
*Dissertations in Forestry and Natural Sciences*

ISBN 978-952-61-3268-6 ISSN 1798-5668

*Dissertations in Forestry and * *Natural Sciences*

**DISSERTATIONS | JENNI TICK | IMAGE RECONSTRUCTION AND MODELLING OF UNCERTAINTIES IN… | No 364**

**JENNI TICK**

**IMAGE RECONSTRUCTION AND MODELLING OF ** **UNCERTAINTIES IN PHOTOACOUSTIC TOMOGRAPHY**

**PUBLICATIONS OF**

**THE UNIVERSITY OF EASTERN FINLAND**

*Photoacoustic tomography is a biomedical *
*imaging modality, where an image of an *
*initial pressure distribution is formed as an *
*inverse problem based on ultrasonic boundary *

*measurements. This thesis introduces new *
*computational methods to solve the inverse *
*problem using Bayesian formalism. The new *

*methods provide quantitative information of *
*the estimated initial pressure, its uncertainties, *

*and improve the accuracy of the estimates by *
*modelling uncertainties in the speed of sound.*

**JENNI TICK**

PUBLICATIONS OF THE UNIVERSITY OF EASTERN FINLAND DISSERTATIONS IN FORESTRY AND NATURAL SCIENCES

N:o 364

**Jenni Tick**

## IMAGE RECONSTRUCTION AND MODELLING OF UNCERTAINTIES IN

## PHOTOACOUSTIC TOMOGRAPHY

ACADEMIC DISSERTATION

To be presented by the permission of the Faculty of Science and Forestry for public examination in the Auditorium SN200 in Snellmania Building at the University of Eastern Finland, Kuopio, on December 20th, 2019, at 12 o’clock.

University of Eastern Finland Department of Applied Physics

Kuopio 2019

Grano Oy Jyväskylä, 2019

Editors: Pertti Pasanen, Jukka Tuomela, Matti Tedre, and Raine Kortet

Distribution:

University of Eastern Finland Library / Sales of publications julkaisumyynti@uef.fi

http://www.uef.fi/kirjasto

ISBN: 978-952-61-3268-6 (print) ISSNL: 1798-5668

ISSN: 1798-5668 ISBN: 978-952-61-3269-3 (pdf)

ISSNL: 1798-5668 ISSN: 1798-5676

Author’s address: University of Eastern Finland Department of Applied Physics P.O. Box 1627

70211 Kuopio, FINLAND email: jenni.tick@uef.fi

Supervisors: Associate Professor Tanja Tarvainen University of Eastern Finland Department of Applied Physics P.O. Box 1627

70211 Kuopio, FINLAND email: tanja.tarvainen@uef.fi Docent Aki Pulkkinen University of Eastern Finland Department of Applied Physics P.O. Box 1627

70211 Kuopio, FINLAND email: aki.pulkkinen@uef.fi Professor Jari Kaipio University of Auckland Department of Mathematics Private Bag 92019

Auckland 1142, NEW ZEALAND email: jari@math.auckland.ac.nz Reviewers: Associate Professor Georg Stadler

New York University

Courant Institute of Mathematical Sciences 251 Mercer Street

New York City NY 10012, USA

email: stadler@cims.nyu.edu Professor Roger Zemp University of Alberta

Department of Electrical and Computer Engineering 9107 – 116 Street

Edmonton

AB T6G 2V4, CANADA email: rzemp@ualberta.ca

Opponent: Research Professor Johanna Tamminen Finnish Meteorological Institute Space and Earth Observation Centre P.O. Box 503

00101 Helsinki, FINLAND email: johanna.tamminen@fmi.fi

Jenni Tick

Image reconstruction and modelling of uncertainties in photoacoustic tomography Kuopio: University of Eastern Finland, 2019

Publications of the University of Eastern Finland Dissertations in Forestry and Natural Sciences 2019; 364

ABSTRACT

Photoacoustic tomography (PAT) is an imaging technique with the advantages of high spatial resolution, strong optical contrast, non-invasiveness, and being based on nonionizing radiation. This hybrid imaging modality has many potential biomed- ical applications such as imaging tissue vasculature and small animal imaging.

In PAT, a short pulse of light is used to illuminate the target of interest. As light propagates in the target, it is absorbed by light absorbing molecules resulting in a localised increase in pressure that is called initial pressure distribution. This pressure propagates through the target as an acoustic wave, and can be measured on the boundary of the target using ultrasound sensors. In the inverse problem of PAT, the initial pressure distribution is estimated from these ultrasound measurements.

Estimation of the initial pressure is an ill-posed inverse problem. Due to the ill- posedness, even small errors or uncertainties in the measurements or modelling can result in significant errors in the solution of the inverse problem. Although many solution methods to the inverse problem of PAT exist, reconstructed photoacoustic images do not necessarily provide quantitative information of the parameters of interest or their uncertainty. In addition, taking into account measurement and modelling errors can be challenging.

In this thesis, a Bayesian approach to the inverse problem of PAT is developed.

With this approach, quantitative information of the parameters and their uncertainty can be achieved. In the Bayesian approach to PAT, the solution of the inverse prob- lem is a conditional probability distribution of the initial pressure that combines information obtained through the measurements, model and prior knowledge. Fur- thermore, in the Bayesian approach, it is possible to incorporate statistical informa- tion about uncertainties and inaccuracies of the models in the solution of the inverse problem. In this thesis, the Bayesian approximation error approach is utilised to model the uncertainties in the speed of sound.

All methods developed in this thesis were evaluated with simulations, and some of the methods were also validated with experimental data. The simulations and ex- periments show that the Bayesian approach is feasible for PAT producing accurate estimates of the initial pressure together with estimates for the reliability. Further- more, the simulations show that by utilising the approximation error approach, it is possible to at least partially compensate for the modelling errors due to uncer- tainties in the speed of sound and to improve the accuracy and feasibility of the solution.

**Universal Decimal Classification:**534.14, 534.22, 535.214, 517.956.27

**INSPEC Thesaurus:** imaging; biomedical imaging; tomography; photoacoustic effect; in-
verse problems; modelling; simulation; errors; error correction; Bayes methods; parameter
estimation; image reconstruction; ultrasonic velocity; uncertainty handling

**Yleinen suomalainen ontologia:**kuvantaminen; fotoakustinen kuvantaminen; tomografia;

inversio-ongelmat; mittaus; mittausmenetelmät; laskentamenetelmät; numeeriset menetelmät;

approksimointi; mallintaminen; simulointi; virheet; bayesilainen menetelmä; estimointi; ää- nennopeus; epävarmuus

ACKNOWLEDGEMENTS

This study was carried out in the Department of Applied Physics at the University of Eastern Finland during the years 2015-2019. The study was financially supported by the Academy of Finland, the doctoral school of the University of Eastern Fin- land, Alfred Kordelin Foundation, Instrumentarium Science Foundation, Finnish Foundation for Technology Promotion, Saastamoinen Foundation, Orion Research Foundation, and Kuopio University Foundation.

I want to express my deepest gratitude to my supervisors Associate Professor Tanja Tarvainen, Docent Aki Pulkkinen, and Professor Jari Kaipio. It has been a privilege to work under your guidance. You have always found me some time when I longed for help, advice or support. In addition, you have taught me a lot about scientific thinking and inspired me to research. Thanks to you, I am now a much better scientist. Furthermore, special thanks to Jari for the opportunity to visit the University of Auckland.

Many thanks to my co-authors Dr. Felix Lucka, Dr. Robert Ellwood, Senior Lecturer Ben Cox, Professor Simon Arridge, and Dr. Jarkko Leskinen. This thesis would not have been possible without your contribution. Special thanks to Felix, Robert, Ben and Simon for all the assistance, knowledge, and skills that I received during my visit to the University College London. Jarkko, in turn, deserves great thanks for introducing me to the wonderful world of photoacoustic measurements.

I am very thankful to the official pre-examiners Associate Professor Georg Stadler and Professor Roger Zemp for taking the time to review and assess my thesis. I greatly appreciate your professional comments and suggestions to improve this the- sis.

I also want to thank the past and present members of the Inverse Problems group, especially those who have shared an office with me at some point. You have provided a pleasant and inspiring working environment where both scientific and non-scientific conversations have been possible. Thanks for the fun time together.

Finally, I am grateful to my parents Anne and Jorma for all the support and encouragement you have given. In addition, thank you for inspiring me to get this highly educated. I also want to thank my sister Henna and my brother Jani for the support and interest towards my work. My warmest thanks go to my partner Martti for believing in me, constant support and motivating me when needed.

Kuopio, November 28, 2019 Jenni Tick

LIST OF PUBLICATIONS

This thesis consist of the present review of author’s work in the field of photoacous- tic tomography and following selection of the author’s publications:

**I** J. Tick, A. Pulkkinen, and T. Tarvainen, ”Image reconstruction with uncertainty
quantification in photoacoustic tomography”,The Journal of Acoustical Society
of America**139(4), 1951–1961 (2016).**

**II** J. Tick, A. Pulkkinen, F. Lucka, R. Ellwood, B. T. Cox, J. P. Kaipio, S. R. Arridge,
and T. Tarvainen, ”Three dimensional photoacoustic tomography in Bayesian
framework”,The Journal of Acoustical Society of America**144(4), 2061–2071 (2018).**

**III** J. Leskinen, A. Pulkkinen, J. Tick, and T. Tarvainen, ”Photoacoustic tomogra-
phy setup using LED illumination”,Proc. SPIE **11077, Opto-Acoustic Methods**
and Applications in Biophotonics IV, 110770Q (2019).

**IV** J. Tick, A. Pulkkinen, and T. Tarvainen, ”Modelling of errors due to speed of
sound variations in photoacoustic tomography using a Bayesian framework”,
Biomedical Physics & Engineering Express,**6(1), 015003 (2019).**

The original articles have been reproduced with permission of their copyright hold- ers. Throughout the overview, these papers will be referred to by Roman numerals.

AUTHOR’S CONTRIBUTION

The publications selected in this thesis are original research papers on photoacous-
tic tomography. The publications of this thesis are the result of joint work between
the author of this thesis and co-authors of the articles. In all publications, the au-
thor implemented the numerical computations and carried out the solution of the
inverse problem. However, in the numerical solution of the wave equation, a k-wave
MATLAB toolbox was used, and the adjoint operator of the PAT forward problem
was developed by co-author Felix Lucka. The author was the main writer of publi-
cations**I,II**and**IV** and the co-writer of publication**III. The author participated in**
the measurements in publication **III. The photoacoustic tomography measurement**
data in publication**II**originated from co-author Robert Ellwood.

### TABLE OF CONTENTS

1 Introduction **1**

2 Photoacoustic tomography **3**

2.1 Measurement setup... 3

2.2 Forward model... 5

2.3 Inverse problem... 5

2.4 Applications of PAI... 8

3 Bayesian inversion methods in PAT **11**
3.1 Bayesian approach to PAT... 11

3.1.1 Prior model... 13

3.1.2 Matrix-free implementation... 14

3.2 Bayesian approximation error approach for modelling of errors... 15

3.2.1 Uncertainties in the speed of sound... 15

4 Summary of the results **17**
4.1 Evaluating the Bayesian approach to PAT with simulations... 17

4.1.1 2D... 17

4.1.2 3D... 23

4.2 Experimental validation of the Bayesian approach to PAT... 26

4.2.1 Fabry-P´erot sensor based PAT setup... 26

4.2.2 LED-based PAT setup... 28

4.3 Compensation of modelling errors... 30

5 Discussion and conclusion **35**

BIBLIOGRAPHY **37**

ABBREVIATIONS

2D Two dimensional

3D Three dimensional

AFM Accurate forward model

CT Computed tomography

IFM Inaccurate forward model

IFM&EM Inaccurate forward model and modelling of errors

LED Light-emitting diode

MAP Maximum a posteriori

MCMC Markov Chain Monte Carlo

MRI Magnetic resonance imaging

PA Photoacoustic

PAI Photoacoustic imaging

PAM Photoacoustic microscopy

PAT Photoacoustic tomography

QPAT Quantitative photoacoustic tomography

TAT Thermoacoustic tomography

NOMENCLATURE

|| · || Norm

*∂s**∂*,_{∂s}^{∂}^{2}_{2} Derivatives

∇^{2} Laplacian

∝ Directly proportional to

arg max Maximum point

·^{T} Transpose

SYMBOLS

*β* Thermal coefficient of volume expansion

c Speed of sound

ec Fixed speed of sound

Cp Isobaric specific heat capacity

CV Specific heat capacity at constant volume

*δ* Dirac delta function

e Noise

ek Unit vector

Ep_{0} Relative error

*e* Approximation error

*ρ* Mass density

Γ Covariance matrix

Γˆ Grüneisen parameter

H Absorbed optical energy

I Identity matrix

*κ* Isothermal compressibility

K Forward operator

l Characteristic length scale

*µ*_{a} Absorption coefficient

M Number of measurements

n Total error

N Number of elements

Ns Number of samples

N Gaussian distribution

*η* Mean

p Acoustic pressure

p_{0} Initial pressure

p_{0,MAP} Maximum a posterioriestimate of the initial pressure

pt Measured pressure signals

Φ Photon fluence

*π*(·) Probability density

r Position vector

**R** Real numbers

*σ* Standard deviation

*σ*^{2} Variance

t Time

T Temperature

### 1 Introduction

In biomedical imaging, one is interested in obtaining images of the internal struc- ture of the biological tissue non-invasively. The most desirable imaging techniques would be those that have both high contrast and high resolution. Various imaging modalities such as computed tomography (CT), magnetic resonance imaging (MRI) and ultrasound have been widely used for biomedical imaging. However, there is still a need for novel imaging techniques that can provide new information about the imaged target. During past few decades, several hybrid or coupled physics modalities have been introduced. Among these is photoacoustic imaging (PAI).

PAI is based on photoacoustic (PA) effect that refers to the generation of acous- tic waves by an optical excitation [1–7]. Although the underlying physics of PA effect was described by Alexander Graham Bell already in 1880 [8], it took more than a century that PAI was investigated for biomedical imaging and the first PAI images were produced [9–17]. The major PAI modalities can be divided into PA tomography (PAT) and PA microscopy (PAM). PAT refers to techniques which re- quire reconstruction algorithms to obtain PA images whereas in PAM, a PA image is directly obtained by mechanically scanning. Depending on whether acoustic fo- cusing or optical focusing is used in the scanning, PAM can further be divided into acoustic-resolution and optical-resolution PAM. Generally, PAT allows deep-tissue imaging with a capability of imaging entire organs while PAM is a relatively super- ficial imaging technique that can be used to image small sections of tissue. In this thesis, the focus is on PAT.

A PAT measurement typically starts by illuminating an imaged object with a short laser pulse. As light propagates in the object, it is partially absorbed by either endogenous chromophores, such as oxygenated and deoxygenated haemoglobin, melanin, lipids, and water, or exogenous contrast agents, such as organic dyes, nanoparticles and fluorescence proteins [18–20]. The absorption generates a small temperature rise that results in a local pressure rise through thermoelastic expan- sion. This initial pressure propagates through the object as an acoustic wave, and it is measured by ultrasound sensors on the boundary of the object. In the inverse problem of PAT, an image of the initial pressure distribution is reconstructed from these boundary measurements which thus gives information on the spatial distribu- tion of chromophores.

Due to PA effect, PAT can combine benefits of optical and acoustic techniques.

That is, PAT has the rich contrast of optical methods and the high spatial resolu- tion of ultrasound techniques. One feature of PAT is the scalability of the spatial resolution and the maximum imaging depth. That is, the spatial resolution can be improved at the expense of the imaging depth and vice versa. At imaging depths of approximately 3 to 7 cm, resolutions of 0.5 to 5 mm can be achieved [21]. Further- more, PAT is a non-ionizing and non-invasive imaging technique. These features make PAT an attractive imaging modality with great potential for biomedical imag- ing where it has found applications for example in breast [22–25], vascular [26, 27]

and small animal imaging [5].

The inverse problem of PAT refers to the problem of forming an estimate of

the initial pressure distribution of an imaged target based on the boundary mea- surements of acoustic pressure waves. This problem has been widely studied and several different reconstruction algorithms have been proposed [28, 29]. Commonly used reconstruction methods include backprojection algorithm [30–33], eigenfunc- tion expansion method [34, 35], time reversal [36–39] and penalized least squares approaches [40–51]. Nevertheless, there are still challenges in solving the inverse problem. Estimation of the initial pressure is generally a moderately ill-posed prob- lem, but it becomes more ill-posed if limited view sensor geometries (sensors do not enclose the imaged target) are used. Ill-posedness means that a problem is sen- sitive to errors and uncertainties both in measurements and in modelling. Thus, in order to produce as good estimates as possible, a measurement situation has to be modelled with sufficient accuracy. In practice, this can be difficult because inherent uncertainties in parameters that are required for modelling, such as mate- rial properties, sensor properties, are always present. In addition, taking modelling uncertainties or errors into account in the solution of the inverse problem can be challenging. Furthermore, a single estimated image is typically obtained as the solution of the inverse problem, so assessing the reliability and uncertainty of the solution can be difficult. To overcome these challenges, improvements in solving of the inverse problem are needed.

This thesis develops Bayesian estimation and uncertainty quantification methods for PAT. The Bayesian approach is a statistical estimation method in which mea- surements, model and prior information are used to infer information about the parameters of primary interest. The aims of this thesis are

1. Formulate a Bayesian framework and computational methods for the inverse problem of PAT. In the Bayesian approach to PAT, one is aiming at obtaining quantitative estimates of the initial pressure distribution and its uncertainties.

2. Develop a Bayesian approximation error modelling based approach to com- pensate for errors caused by an unknown speed of sound in PAT.

3. Evaluate the feasibility of the developed methods with simulations and exper- imental data.

This thesis is organized so that Chapter 2 describes the forward and inverse problem of PAT and reviews the research. In Chapter 3, the Bayesian framework for PAT and its utilisation in modelling of errors are described. The results of the publications are reviewed in Chapter 4. Finally, Chapter 5 gives the discussion of the results and final conclusions.

**Figure 2.1:** Principle of PAT. a) Illumination of an imaged object with a light pulse.

b) Absorption of light and a corresponding localised pressure increase. c) Measure- ment of the pressure wave with ultrasound sensors on the surface of the object.

### 2 Photoacoustic tomography

In this chapter, the principle of PAI and the forward and inverse problems of PAT are described. Furthermore, solution methods and challenges related to the inverse problem are reviewed. Finally, applications of PAI are introduced. Further infor- mation about history of PAT, physical principles, measurement systems, solution methods, and applications in biomedicine can be found in several comprehensive reviews, see e.g. [1–7, 28, 52].

2.1 MEASUREMENT SETUP

The procedure of a PAT measurement situation illustrated in Figure 2.1 is as follows.

First, a short (nanosecond scale) pulse of light is used to illuminate an imaged tar- get. By tuning the wavelength of light, the penetration depth inside the target and absorption by chromophores of interest can be altered. Typically, light in the visible or near-infrared region is utilised [3]. In this spectral region, light has its maximum depth of penetration in tissue due to the relatively weak absorption of chromophores in tissue. In the majority of the PAT measurement setups, the illumination is imple- mented with solid state lasers, but also laser diodes [53–56], light-emitting diodes (LEDs) [57–60] and broadband xenon lamps [61] have been used as a light source.

From the light source, light is guided to the target using for example optical lenses, mirrors, optical fibers or fiber bundles.

Light propagation in a random media such as biological tissue can be described using analytic theory and transport theory [62]. In the analytical theory, Maxwell’s

equations are used to describe the transportation of particles. The transport theory can be derived using stochastic approach, which models individual particle interac- tions [63], or deterministic approach, which describes the transportation of particles using partial differential equations such as the radiative transport equation and its approximations [62, 64]. As the light pulse propagates within the imaged target, some of its energy is absorbed by chromophores. The absorbed energy at pointris given by

H(r) =*µ*a(r)_{Φ}(r)_{,} _{(2.1)}
where *µ*a(r) is the absorption coefficient and Φ(r) is the photon fluence. The ab-
sorbed energy is converted into heat and the temperature of the target increases
slightly, which temperature increase is given by

T(r) = ^{H}(r)

*ρC*V , (2.2)

where*ρ*is the mass density andCV is the specific heat capacity at constant volume.

This temperature rise further causes thermoelastic expansion of the target, which in turn induces an excess pressure. If the duration of the light pulse is less than both the thermal and stress confinement times, the local fractional volume expansion and thermal diffusion are negligible, and the pressure increase can be considered to be instantaneous. This initial pressure can be written as

p0(r) = * ^{βT}*(r)

*κ*

= * ^{βH}*(r)

*ρC*V

*κ*

= ^{βc}

2H(r)

Cp

= ΓH^{ˆ} (r), (2.3)

where *β* is the thermal coefficient of volume expansion, *κ* is the isothermal com-
pressibility,cis the speed of sound,Cpis the isobaric specific heat capacity, and ˆΓis
the Grüneisen parameter that indicates the efficiency of conversion of the absorbed
optical energy into pressure [2, 7]. The generated pressure relaxes and propagates
as an acoustic wave that is known as PA wave. Generally, the amplitude of the wave
is relatively low but the spectral content of the wave is broad (tens of MHz).

The PA signals are acquired by ultrasound sensors at multiple locations around the object. Sensors may be in contact to the surface of the imaged object or away when they need to be coupled. Three most commonly used sensor geometries are spherical [23,65,66], cylindrical [67–72] and planar detection geometry [73–80] where the ultrasonic detection follows a (truncated) spherical surface, cylindrical surface, or plane, respectively. While spherical or cylindrical sensor geometries requires ac- cess to all sides of the object, planar sensor geometries can be utilised in a wider range of imaged targets. However, the image quality provided by spherical and cylindrical geometries is superior to the image quality provided by a planar de- tection geometry due to their larger apertures. The sensor geometries can be im- plemented by using multiple sensors or by moving a single sensor. Multielement arrays have also been applied. Moving a single sensor around the imaged object can be time-consuming. Therefore, multiple sensors or arrays are preferred if a fast data acquisition is required. Commonly used sensors can be divided into piezoelec- tric [81, 82] and optical sensors [75, 83–86].

2.2 FORWARD MODEL

The forward problem of PAT is to solve the time-varying pressure at the sensors surrounding the object when the acoustical properties and initial pressure distribu- tion are given. In PAT, the propagation of the PA wave can generally be modelled using equations of linear acoustics, because PA waves are typically of small ampli- tude. For soft biological tissue, it is generally assumed that the medium is isotropic and quiescent and that the shear waves can be neglected. Further, if the medium is assumed to be non-attenuating, the time evolution of the induced PA wave fields can be modeled using a wave equation

*∂*^{2}

*∂t*^{2} −c(r)^{2}∇^{2}^{}p(r,t) =0
p(r,t=0) =p0(r)

*∂t**∂*p(r,t=0) =0

(2.4)

wherep(r,t)is the acoustic pressure at pointrand timet,p_{0}(r)is the initial pressure
distribution, andc(r)is the speed of sound [2, 6].

However, in many practical applications, PA waves travel through a complex heterogeneous biological medium that cannot be considered non-attenuating. A model of ultrasound propagation that accounts acoustic attenuation has been used in [39,87–91]. Furthermore, neglecting the shear waves is not valid assumption when the imaged target includes for example bone. Therefore, an elastic wave equation model has been considered [92, 93]. In this thesis, the wave equation (2.4) is used to model the wave propagation.

The numerical solution of the wave equation (2.4) can be obtained using for example finite difference or finite element methods [94]. These methods are su- perb for many applications since they are simple, adaptable and easily parallelised.

However, these methods can become computationally expensive when broadband ultrasound fields in the time domain are considered [95]. Alternatively, a k-space pseudospectral method can be utilised [96]. This method enables a computation- ally efficient solving of the wave equation by combining the spectral calculation of spatial derivatives with a temporal propagator expressed in the spatial frequency domain ork-space. In this thesis, thek-space method implemented with the k-Wave MATLAB Toolbox [95] is used for the numerical solution of the wave equation.

2.3 INVERSE PROBLEM

In PAT, estimates of the initial pressure distribution, also called photoacoustic im- ages, are of principal interest. Estimation of the initial pressure from the measured PA waves is the inverse problem of PAT. This problem has an unique solution al- most for all practical acoustic detection surfaces, but the solution can sometimes be unstable [29]. Thus, the problem is ill-posed. However, the inverse problem is only mildly ill-posed in full view sensor geometries in which the imaged object is fully surrounded by sensors on a closed surface. In this case, qualitative accurate esti- mates of the initial pressure can be achieved. Sensor geometries that do not enclose the object are called limited view geometries. In these sensor geometries, only a part of the wavefront is recorded which can result in images that suffer from blurry edges, loss of details and reduced image quality. More specifically, sharp boundaries of the object can be reconstructed accurately if they are facing the detection surface

**Figure 2.2:** Two examples of limited view imaging scenarios with sensors on an
arc (left) and on a line (right). The region enclosed by the sensors (gray shaded
area), the inclusion boundaries that are reconstructed accurately (solid line), and the
inclusion boundaries that are blurred (dashed line). Image is adapted from [6, 97].

but suffer from blurring if they are perpendicular to the detection surface [6, 97].

This is illustrated in Figure 2.2.

The inverse problem of PAT has been widely studied and a large number of solution methods are available. The approaches can be classified either analytical or numerical methods. The former aim to solve the problem analytically and the latter utilise the numerical solution of the problem.

The analytical methods include backprojection algorithm [30–33] and eigenfunc- tion expansion method [34, 35]. The backprojection approach is based on the in- verse circular/spherical Radon transform, and the initial pressure is reconstructed by summing up the backprojected measured pressure signals with appropriate time delays. In the eigenfuction expansion method, the initial pressure is obtained as the series solution and series coefficients are calculated from the measured pressure signals.

The numerical methods include time reversal [36–39], penalized least squares
approaches [40–51], and Bayesian approach [98], **I**and **II. In the time reversal al-**
gorithms, an image is reconstructed by running a numerical model of the forward
problem backwards in time. That is, the measured acoustic pressure signals are
backpropagated into the domain in the time-reversed order. The approaches based
on penalized least squares perform an image reconstruction by minimizing the
sum of the least square error between the measured signal and signals predicted
by the photoacoustic forward model and a regularizing penalty functional such as
Tikhonov [40, 42, 46] or total variation [41, 43, 45, 46, 48]. The Bayesian approach is a
statistical inversion method and it is reviewed more in detail in Section 3.1.

Generally speaking, the analytical methods can provide an ease of implemen- tation and they can have a low computational cost. However, they are limited to specific geometries such as spherical, cylindrical, and planar acoustic detection sur- faces and require that the PA signals are densely sampled on the detection surface.

In contrast, the numerical methods are more versatile, albeit computationally more intensive. In addition, the above described methods, apart from the Bayesian ap- proach, produce a single estimated image as the solution to the inverse problem.

Thus, the reliability, uncertainty, and errors of the solution can be challenging to assess.

Despite many solution methods, obtaining an accurate solution for the inverse problem of PAT is not a straightforward task in practice. This is partly due to ill- posed nature of the problem, which is why a careful modelling of the measurement situation is required. Next, challenges related to solving the inverse problem of PAT in practice are discussed.

In practice, limited view sensor geometries are often required to be used due to constraints associated with an experimental setup or an imaging application, although they suffer from the reduced image quality. Thus, methods for improv- ing the image quality have been developed. In the backprojection approach, the image quality of the limited view reconstruction has been improved by adding a weighting factor from a smoothing function to backprojection signals [99]. In the penalized least squares approaches, the artifacts resulting from the limited view sce- nario have been reduced by utilising prior information in form of a regularization term [100–102]. Especially, prior structural information has been shown to improve the quality of the reconstructed images significantly [103, 104]. In the time rever- sal approaches, an enhanced version of time reversal such as iterative time reversal can be used to improve the quality of the reconstructions in the limited view se- tups [45, 105, 106]. In the Bayesian approach, prior information can be incorporated into the image reconstruction procedure.

In the analytical reconstruction algorithms, it is assumed that the ultrasound sensors are ideal point-like sensors with an infinite bandwidth. However, real ul- trasound sensors have a finite aperture and a finite bandwidth. Both these char- acteristics of the sensor can cause blurring of the reconstructed images [107, 108].

However, in the numerical solution methods, the properties of sensors can be mod- elled. This has resulted in significantly reduced blurring and further the enhanced image quality in the penalized least squares approaches [109–111].

Reconstruction algorithms rely on having accurate knowledge of the speed of sound in the medium. This assumption is somewhat problematic in real applica- tions where the speed of sound is typically unknown. Basically, the speed of sound could be estimated from an adjunct measurement such as ultrasound tomography measurements [112–115] or jointly with the initial pressure distribution [116–120].

However, an implementation of adjunct imaging and non-uniqueness of the joint estimation problem cause practical challenges. Thus, a preassigned nominal value for the speed of sound is more commonly used, for example by setting it equal to the speed of sound in water or to average speed of sound in tissue. Using an ap- proximate speed of sound that inaccurately addresses the true speed of sound dis- tribution can yield estimates which contain severe artefacts [38, 121]. Some methods have been proposed to mitigate these artefacts, but these approaches compensate well mainly such acoustic heterogeneity-induced artifacts that originate from small speed of sound variations [122–126].

The most approaches for the inverse problem have been using the wave equation (2.4) as the forward model. In practice, this can be too simplified model for the wave propagation in tissue, which can result in significant errors in the solution of the inverse problem due to modelling errors. Thus, models incorporating more realistic tissue properties have been studied [39, 49, 88, 92, 127–129]. Although the utilisation of these models can reduce modelling errors and thus improve the accuracy of the solution of the inverse problem, it can also make solving of the inverse problem more challenging due to increased complexity of the model.

2.4 APPLICATIONS OF PAI

Photoacoustic imaging is an excellent modality for imaging blood vessels due to its high sensitivity to haemoglobin with certain wavelenghts [26, 27]. PAI can detect both tiny and large vessels so it can provide comprehensive imaging of the vascu- lature. In addition, a structure of vessels can be visualised using intravascular PA imaging catheters [130–133]. Since increased vascularity, high total blood concentra- tions and abnormally low haemoglobin oxygen saturation indicate the presence of tumours, PAI can potentially be used for tumour imaging [134].

PAI is well suited for molecular imaging which measures chemical and biological processes at the molecular and cellular level [18]. Molecular imaging is based on exogenous chromophores, which either change their optical properties if there are changes in the molecular level or have specific spectral signatures and target specific molecules.

In addition to structural and molecular imaging, PAI can also provide functional information. One important functional parameter that PAI is able to measure is blood oxygen saturation [135]. Since blood oxygen saturation is determined by the concentrations of oxyhaemoglobin and deoxyhaemoglobin, PAI is able to pro- duce an image of blood oxygen saturation. Other functional capability of PAI is the measurement of blood flow velocity [136]. In addition, PAI can provide functional information in form of temperature [137–139]. The Grüneisen coefficient changes with temperature, thus PAI can provide images of temperature distributions in tis- sue [139].

Another major application of PAI is small animal imaging, which enables study- ing human disease processes, monitoring treatment and developing new treatments [5]. The small size of the animal means that penetration depth requirements are generally modest and signals can usually be recorded around the whole body of the animal. Thus, high-resolution images can be obtained even deeper anatomical structures. In small animals, for example viscera [65, 68, 140–143], heart [144, 145], skin [74, 146, 147], and brain [67, 148–150] have been studied.

Breast imaging is a promising application of PAI [22–25]. In fact, breast imaging is one of the most advanced clinical applications of PAI and clinical measurement devices are available. PAI can provide comprehensive information of the vasculature and tumour properties. Thus, PAI is well suited to diagnostic imaging of the breast.

Another potential application of PAI is brain imaging [151]. Imaging of brain is challenging due to strong optical scattering and acoustic attenuation by the skull.

Thus, studies have been mainly performed using small animals due to their thin skull. In these studies, both brain structure and functionality have been imaged [67, 79, 149, 152]. In addition, brain tumours and disorders have been considered [153–

157]. Thus far, imaging of human brains have not been performed, but a few studies promotes its feasibility, especially if human neonates would be imaged subjects [72, 158].

PAI has also been applied in dermatologic imaging [74, 146, 147, 159, 160]. Super- ficial vasculature of small animals and humans have been imaged. In addition, skin cancers, burns and other skin lesions have been detected with PAI. Recently, PAI has been applied in diagnosis and monitoring of psoriasis [161].

The other possible biomedical applications include gynecological [162, 163], uro- logic [164, 165], endocrinological [166–168], and ocular [169] imaging. In addi- tion, PAI have been utilised in diagnosing joint diseases [170–172] and Crohn dis- ease [173].

In addition to imaging, PAI can be used as an assisting method. PAI can be utilised in planning and monitoring treatments and thus it could help to improve treatments and outcomes. In addition, PAI can be an useful tool for assessing the efficacy of treatments. Examples of applications are guidance of therapy, surgery and biopsy [134, 174–181].

### 3 Bayesian inversion methods in PAT

In this chapter, the Bayesian approach to the inverse problem of PAT is described
following publication **I. In addition, extending this approach to three dimensions**
(3D) is described following publication**II. Finally, the Bayesian approximation error**
approach and its utilisation in modelling uncertainties in the speed of sound is
described following publication**IV.**

3.1 BAYESIAN APPROACH TO PAT

In the inverse problem of PAT, an unobservable parameterp0(parameter of interest) is estimated based on the measurements pt. An observation model describes how these depend on each other. Since the practical measured pressure signals always contain noise that is commonly assumed to be additive, the discrete observation model of PAT corresponding to the forward model (2.4) can be written in form

pt=K(c)p0+e, (3.1)

where pt∈**R**^{M}is a vector composed of the acoustic pressure waves sampled at the
sensors at discrete time points, p0 ∈**R**^{N} is the discrete initial pressure distribution,
K(c)∈**R**^{M×N} is the linear operator which maps the initial pressure distribution to
the measurable data by discretizing the wave equation (2.4), ande∈**R**^{M}represents
the unknown noise. In practice, the operatorK(c)can be assembled by computing
the impulse response of a discrete system approximating (2.4) for each of the ele-
ments in the reconstructed area and placing the acoustic output on the columns of
K(c). In this thesis, the impulse responses were computed using the k-space time
domain method implemented in the k-Wave toolbox [95].

In the Bayesian framework [182, 183], all parameters are modelled as random variables, and information related to these parameters is expressed by probability distributions. The solution of the inverse problem is the posterior density that is a conditional probability density for the parameters of interest in each element of the domain. The posterior density extracts information and assesses uncertainty about the variables based on information obtained through the inference of measurements, model, and prior model for unknown parameters. In the case of PAT, the posterior density reflects information on the unknown initial pressure distribution p0 given the acoustic pressure measurements pt.

According the Bayes’ theorem, the posterior density follows proportionality
*π*(p0|pt)_{∝}*π*(p0)*π*(pt|p0). (3.2)
where *π*(p_{0}) is the prior probability density and *π*(pt|p_{0}) is the likelihood den-
sity. The prior density models what is known of the parameters of interest before
the measurements, while the likelihood describes the probabilities of measurement
outcomes over all choices of parameters .

Let us consider a linear observation model (3.1) and assume that p0 and e are mutually independent. Then, the likelihood becomes

*π*(pt|p0) =*π*e(pt−K(c)p0), (3.3)

where*π*eis the probability density of the noisee. The posterior density can now be
written as

*π*(p0|pt)∝*π*(p0)*π*_{e}(pt−K(c)p0). (3.4)
Let us further assume that the noiseeand initial pressurep0 are Gaussian dis-
tributed i.e. e ∼ N(*η*e,Γe) and p0 ∼ N(*η*p_{0},Γp_{0}) where *η*e and Γe are the mean
and covariance of the noise, respectively, and *η*p_{0} and Γp_{0} are the mean and co-
variance of the prior model, respectively. Now, in the case of a linear observation
model and Gaussian distributed noise and prior, the posterior density is a Gaussian
N(*η*_{p}_{0}_{|p}_{t},Γ_{p}_{0}|p_{t})and is described by the mean and covariance

*η*_{p}_{0}_{|p}_{t} = _{Γ}_{p}

0|pt

K(c)^{T}_{Γ}^{−1}_{e} (pt−*η*e) +_{Γ}^{−1}_{p}

0*η*p_{0}

(3.5)
Γ_{p}_{0}_{|p}_{t} = ^{}K(c)^{T}Γ^{−1}_{e} K(c) +Γ^{−1}_{p}_{0}−1

. (3.6)

It should be noted that the posterior density has this closed form solution only when an observation model is linear and noise and prior are Gaussian distributed.

In other situations, i.e. non-linear observation model and/or non-Gaussian dis- tributions, the posterior density can in principle be explored by employing Markov Chain Monte Carlo (MCMC) technique, which provides a representative set of sam- ples from the posterior density. However, in tomography, the posterior density is often high dimensional, and thus MCMC methods can be computationally pro- hibitively too expensive. Therefore, point estimates of the posterior density are typ- ically computed. The most commonly used point estimate is amaximum a posteriori (MAP) estimate that is defined as

p_{0,MAP}=arg max

p_{0}

*π*(p_{0}|pt). (3.7)

In Gaussian case, the MAP estimate is same as the mean of the posterior density i.e.

p_{0,MAP}=*η*_{p}_{0}_{|}_{p}_{t}.

In addition to obtaining an estimate of the parameter of interest, the posterior
density also provides information on the uncertainty of the estimate. This informa-
tion can be used to assess the accuracy of the estimate. From the posterior covariance
posterior variances *σ*_{p}^{2}

0|p_{t} = diag{Γ_{p}_{0}|p_{t}}and further standard deviations *σ*_{p}_{0}_{|p}_{t} for
the each unknown parameter can be extracted. In addition, the reliability of the
estimates can be assessed by computing credible intervals

h

*η*_{p}_{0}_{|p}_{t}_{,k}−*ασ*_{p}_{0}_{|p}_{t}_{,k}, *η*_{p}_{0}_{|}_{p}_{t}_{,k}+*ασ*_{p}_{0}_{|p}_{t}_{,k}

i, (3.8)

where *η*_{p}_{0}_{|p}_{t}_{,k} is the value of the posterior mean*η*_{p}_{0}_{|p}_{t} in thekth element,*σ*_{p}_{0}_{|p}_{t}_{,k} is
the value of *σ*_{p}_{0}_{|p}_{t} in thekth element and *α* = 1, 2, 3 corresponds to 68.3%, 95.5%

and 99.7% intervals, respectively. Also the posterior marginal densities can be in- formative for inspecting the uncertainty of the estimate. Since the posterior density is a Gaussian, all marginal densities are Gaussian. The marginal density of the kth element is defined as

p_{0,k}|pt∼ N(*η*_{p}_{0}_{|p}_{t}_{,k},Γp_{0}|p_{t},kk), (3.9)
where*η*_{p}_{0}_{|p}_{t}_{,k}is the value of*η*_{p}_{0}_{|p}_{t} in thekth element andΓ_{p}_{0}|p_{t},kk is the value of the
kth diagonal element of the matrixΓ_{p}_{0}|p_{t}.

**Figure 3.1:** Two sample draws of 2D spatial distributions (10×10 mm) from the
white noise (first column), Ornstein-Uhlenbeck (second column) and squared ex-
ponential (third column) priors with the mean*η*p0 =5 and the standard deviation
*σ*p0 = 10/6. The characteristic length scalel = 2 mm was used for the Ornstein-
Uhlenbeck and squared exponential prior.

3.1.1 Prior model

A prior model describes our pre-existing knowledge of the imaged target in form of a probability density. That is, it describes what values are probable for each parameter of interest and how these values are distributed. The prior model should be chosen in such a way that the imaged target is well supported by the model.

Commonly used prior models are Gaussian priors [182]. Three examples of
Gaussian priors are a white noise prior, Ornstein-Uhlenbeck prior and squared ex-
ponential prior [184,185]. Sample draws of 2D spatial distributions from these priors
are shown in Figure 3.1. Gaussian priors can be described by their means*η*_{p}_{0} and
covariance matrices Γp_{0}. If the range of the unknown parameter values is roughly
known, this information can be utilised to choose the mean*η*p_{0} of the prior. Corre-
spondingly, information related to the expected spatial distribution of the parame-
ters can be utilised to choose the covariance matrixΓp_{0}. If the estimated parameters
are assumed to be non-smooth (the parameters are independent of each other or
have no spatial correlation), the white noise prior can be an appropriate choice. In
the white noise prior, theΓp_{0} is defined as

Γp_{0},ij =*σ*_{p}^{2}_{0}*δ*_{ij}, (3.10)

where i and j are the element indices, *σ*p_{0} is the standard deviation, and*δ*ij is the
Dirac delta function. The prior knowledge on the range of the unknown parameter
values can be utilised to choose the standard deviation *σ*_{p}_{0}. On contrary, if the
parameters are assumed to be somewhat smooth, the squared exponential prior or
Ornstein-Uhlenbeck prior can be an useful choice. In the squared exponential prior,

theΓp_{0} is of form

Γp_{0},ij=*σ*_{p}^{2}_{0}exp −kr_{i}−r_{j}k^{2}
2l^{2}

!

(3.11)
whereas in the case of the Ornstein-Uhlenbeck prior theΓp_{0} is defined as

Γp_{0},ij=*σ*^{2}_{p}_{0}exp

−kr_{i}−r_{j}k
l

, (3.12)

where r_{i} and r_{j} are the element positions of the element indicesi and j and l is
the characteristic length scale which controls the spatial range of correlation. The
characteristic length scalelcan be chosen based on prior information about internal
structures of the imaged target and size of these structures.

In PAT, some spatial correlation in parameter values can be expected since pre- sumably properties of biological tissue types are relatively homogeneous. In addi- tion, it is possible that the target is composed of heterogeneities separated by sharper edges such as blood vessels. Thus, the squared exponential or Ornstein-Uhlenbeck prior can be regarded as more appropriate presumption for PAT than the white noise prior.

3.1.2 Matrix-free implementation

In a large dimensional inverse problem, the closed form matrix presentation i.e.

(3.5)-(3.6) might be impractical to evaluate and store. Therefore, a matrix-free method
for the computation of the posterior mean and covariance is needed. This can be
achieved by solving linear systems corresponding to the mean and covariance us-
ing iterative linear solvers, such as conjugate gradient methods, and evaluating the
multiplication by the forward modelKwith the forward operator [96], the multipli-
cation by the transpose of the forward modelK^{T} with the adjoint operator [45] and
the multiplication by the covariance of the prior Γp_{0} with a 3D convolution using
the fast Fourier transform as follows. The mean of the posterior density (3.5) can be
computed by solving the system of linear equations

Cη_{p}_{0}|p_{t} =d, (3.13)

where

C=_{Γ}_{p}_{0}K(c)^{T}_{Γ}^{−1}_{e} K(c) +I, (3.14)
d=Γp_{0}K(c)^{T}Γ^{−1}_{e} (pt−*η*_{e}) +*η*_{p}_{0} (3.15)
and Iis an identity matrix. Correspondingly, thekth column of the posterior covari-
ance (3.6) can be computed by solving the linear system

CΓ_{p}_{0}|pt,k =_{Γ}_{p}_{0}e_{k}, (3.16)
whereCis as in (3.14) ande_{k}is a unit vector with value one at thekth element and
zeros elsewhere. In this thesis, the linear systems (3.13) and (3.16) are solved using
a biconjugate gradient stabilized (l) method [186] built-in MATLAB.

3.2 BAYESIAN APPROXIMATION ERROR APPROACH FOR MOD- ELLING OF ERRORS

A forward model that describes the dependence of unknown parameters and mea- surements is an essential part of the solution of an inverse problem. However, the model may be inaccurate and it may contain parameters that are not well known. In addition, some of the unknown parameters may not be of direct interest and they are referred to as nuisance parameters. If the model is an approximation to the accurate physical model, this may affect unfavourably the solution of the inverse problem due to its ill-posedness. However, a consideration of the inverse problem utilising the Bayesian approximation error modelling [182] allows accounting for modelling errors.

The Bayesian approximation error modelling has been applied to handle many
kinds of approximation and modelling errors in a variety of inverse problems. Ex-
amples in optical and acoustic imaging modalities include quantitative PAT (QPAT)
where the Bayesian approximation error approach has been utilised in the marginal-
ization of scattering [187], mitigation of inaccuracies due to an acoustic solution
method [188] and compensation of discretisation errors [189]. In diffuse optical
tomography, a model reduction [190–193] and compensating uncertainties in mea-
surement geometry, boundary shape and model parameters has been considered
[194–196]. In full-waveform ultrasound tomography, errors due to domain trunca-
tion and approximate boundary models were treated [197]. In PAT, the Bayesian ap-
proximation error modelling has been used to model and compensate errors caused
by uncertainties in the speed of sound**IV**(Section 3.2.1) and ultrasound sensor lo-
cations [198].

3.2.1 Uncertainties in the speed of sound

In PAT, the speed of sound in the imaged object is commonly considered as a known nuisance parameter. However, the accurate speed of soundc, that can be constant or spatially varying, is typically not available in practical applications. Therefore, some nominal value ec is used in the inverse problem, which can cause modelling errors. However, these errors can be compensated by expressing the observation model utilising the Bayesian approximation error modelling [182] in the form

pt = K(_{e}c)p0+*e*+e

= K(_{e}c)p0+n, (3.17)

where *e* = K(c)p0−K(_{e}c)p0 is the approximation error describing the modelling
error, or discrepancy, between the exact forward model with the accurate speed of
soundcand the approximate forward model with an inexact speed of soundecand
n=*e*+eis the total error.

When the total errorn is approximated to be mutually independent of the un- known p0, the observation model (3.17) leads to a likelihood density

*π*(pt|p0) =*π*n(pt−K(_{e}c)p0), (3.18)
where*π*nis the probability density of the total errorn. Further, the posterior density
can be written as

*π*(p0|pt)∝*π*(p0)*π*n(pt−K(_{e}c)p0). (3.19)

Let us assume that noisee is Gaussian distributed e ∼ N(*η*e,Γe), and approx-
imate the modelling error *e* as a Gaussian *e* ∼ N(*η** _{e}*,Γ

*e*). Thus, the total error n is a Gaussian distribution n ∼ N(

*η*

_{n},Γn), where

*η*

_{n}=

*η*

_{e}+

*η*

*and Γn = Γe+Γ*

_{e}*e*

leading to a Gaussian approximation for the likelihood. With this likelihood, a
Gaussian prior, and a linear observation model, the posterior density is a Gaussian
N(*η*_{p}_{0}_{|p}_{t},Γ_{p}_{0}|p_{t})given by the mean and covariance

*η*_{p}_{0}_{|p}_{t} = _{Γ}_{p}_{0}_{|p}_{t}^{}K(_{e}c)^{T}_{Γ}^{−1}_{n} (pt−*η*n) +_{Γ}^{−1}_{p}_{0}*η*p_{0}

(3.20)

Γ_{p}_{0}|p_{t} = ^{}_{K}(_{e}_{c})^{T}_{Γ}^{−1}_{n} _{K}(_{e}_{c}) +_{Γ}^{−1}_{p}

0

−1

. (3.21)

The realisation of the approximation error*e*is unknown since its value depends
on the actual unknown p_{0}and inaccurately known nuisance parameterc. However,
samples for the modelling error can be simulated utilising ’teaching’ distributions
of unknowns, which are probability densities including the prior information about
the unknowns, as follows. First, a set of samples {p^{(l)}_{0} ,l = 1, ...,Ns}and {c^{(l)},l =
1, ...,Ns}are drawn from the teaching distribution of the initial pressure and speed
of sound, respectively. Then, samples of the approximation error are computed
using

*e*^{(l)}=K(c^{(l)})p^{(l)}_{0} −K(_{e}c)p_{0}^{(l)}. (3.22)
From these samples, the mean*η** _{e}*and covarianceΓ

*e*of the approximation error can be estimated as

*η** _{e}*=

^{1}Ns

N_{s}
l=1

### ∑

*e*^{(l)} (3.23)

Γ*e* = ^{1}
Ns−1

Ns

l=1

### ∑

(*e*^{(l)}−*η** _{e}*)(

*e*

^{(l)}−

*η*

*)*

_{e}^{T}. (3.24)

### 4 Summary of the results

This chapter summarizes the main results of the publications constituting this thesis.

First the Bayesian approach is validated with numerical simulations in Section 4.1.

These results are based on publications**I**and**II. Secondly, results of the experimen-**
tal validation of the approach are presented in Section 4.2. These results are based
on publications **II**and **III. Finally, in Section 4.3, results of numerical simulations**
on modelling of uncertainties in the speed of sound using the approximation error
modelling are presented. These results are based on publication**IV.**

4.1 EVALUATING THE BAYESIAN APPROACH TO PAT WITH SIM- ULATIONS

The feasibility of the Bayesian approach to PAT was first evaluated with numerical
simulations. In publication**I, simulations were done in 2D whereas in publicationII**
simulations were conducted in 3D. In the 2D simulations, the approach was inves-
tigated utilising two different distributions of a Gaussian prior in various imaging
situations including limited view measurement setups. In addition, different sensor
properties such as a size and bandwidth were considered. In these simulations, the
entire posterior density was determined and inspected. In addition to the poste-
rior standard deviation, the reliability of the estimates was assessed by computing
marginal densities and credible intervals. In the 3D simulations, the performance of
the approach was demonstrated in different sensor geometries, and point estimates
for the image reconstruction and uncertainty quantification were computed. In all
simulations, the quantitative accuracy of the posterior mean estimates was evalu-
ated by computing the relative errors of the estimates with respect to the true initial
pressure distribution using

Ep_{0} =100%· kp0−p0,MAPk

kp0k ^{,} ^{(4.1)}

where p0is the simulated initial pressure distribution and p0,MAP is the estimated value.

4.1.1 2D

In the 2D simulations, a square domain with a side length of 10 mm was used.

The medium was assumed to be non-attenuating with a constant speed of sound 1500 m/s. The simulated initial pressure distribution is shown in Figure 4.1. Sen- sors were arranged around the domain (full view), on two adjacent sides of the domain (two side) or on one side of the domain (one side). Ideal point-like sen- sors (ideal sensors), sensors of width 0.5 mm (finite sized sensors), and sensors with limited bandwidth (bandlimited sensors) were investigated. In the case of the ban- dlimited sensors, the frequency response of the sensors was modeled as a bandpass filter with−6 dB bandwidth of 8 MHz and central frequency of 6 MHz. The simu- lated measurement data was generated by solving the wave equation (2.4) using the

**Figure 4.1:** The simulated (true) initial pressure distribution given in arbitrary
units. The white square and circle indicate the locations where the marginal densi-
ties are plotted and the dashed line indicates the location where the credible interval
is plotted.

k-space time domain method and adding a 1% Gaussian noise to it. In the data sim-
ulation, the domain was discretised into 300×300 square pixels (pixel side length
of 33*µm) and 283 temporal samples (sampling frequency 20 MHz*)were used.

In the inverse problem, the mean and covariance of the posterior density were
computed using (3.5) and (3.6), respectively. In addition, credible intervals were
computed on a diagonal cross section of the domain (a dashed line in Figure 4.1) us-
ing (3.8). Furthermore, marginal densities were calculated in two locations (a square
and circle in Figure 4.1) using (3.9). In the computations, the domain was discretised
into 120×120 square pixels (pixel side length of 83*µm). The Ornstein-Uhlenbeck*
and white noise priors were used as the prior model for the initial pressure. For
both priors, the mean was set as the expected mean value of the initial pressure
*η*_{p}_{0} = 5 and the standard deviation was set as*σ*_{p}_{0} = 2.5, which means that 99.7%

of the initial pressure values were expected to be normally distributed within the range[−2.5, 12.5]. For the Ornstein-Uhlenbeck prior, the characteristic length scale l=1.25 mm was used. In addition, the accurate noise statistics were assumed.

The simulations show that the sensor geometry affects the initial pressure esti- mates and their uncertainties. In the full view sensor geometry, the inverse problem is only mildly ill-posed. Thus, the estimates of the initial pressure distribution are qualitatively similar with true one (first column of Figure 4.2). In addition, the es- timates are quantitatively accurate and the relative errors are smallest (Table 4.1).

Furthermore, the standard deviations are small indicating small uncertainty of the estimates (first column of Figure 4.2). In limited view sensor geometries, the inverse problem becomes more ill-posed. Thus, the estimates of the initial pressure dis- tribution suffer from artefacts and distortions (second and third column of Figure 4.2). The distortion of the estimates increases at larger distances from the sensors.

In addition, the quantitative accuracy of values in these distorted areas is reduced.

This results into higher relative errors (Table 4.1). However, the uncertainty of the estimates also increases, especially in the distorted areas, as the number of detection edges decreases. In all cases, the true value of the initial pressure is within the error limits obtained in the Bayesian approach (Figures 4.3 and 4.4). This indicates that the Bayesian approach can provide reliable uncertainty estimates.

The simulations also indicate that the prior has an influence on the solution of the inverse problem (Figure 4.2). In the full view sensor geometry case, the

**Figure 4.2:** The posterior mean (top block) and standard deviation (bottom block)
obtained using the Ornstein-Uhlenbeck prior (first row) and the white noise prior
(second row) in the case of the ideal sensors. The columns from left to right represent
the full view (first column), two side (second column) and one side (third column)
sensor geometries. The red dots in the first row images indicate the locations of the
sensors.

**Table 4.1:** The relative errors Ep0(%) of the estimated posterior mean obtained
using the full view, two side and one side sensor geometries in 2D. The ideal sen-
sors (ideal), finite sized sensors (finite sized), and bandlimited sensors (bandlim-
ited) were considered using the Ornstein-Uhlenbeck (OU) and/or white noise (WN)
prior.

Ideal Finite sized Bandlimited

OU WN OU OU

Full view 13 13 23 14

Two side 15 16 33 19

One side 34 36 56 51

**Figure 4.3:** The true initial pressure distribution (black solid line) and the posterior
mean (red dashed line) with 99.7% credible intervals (purple area) on a diagonal
cross-section shown in Figure 4.1 obtained using the Ornstein-Uhlenbeck prior (first
row) and the white noise prior (second row) in the case of the ideal sensors. The
columns from left to right represent the full view (first column), two side (second
column) and one side (third column) sensor geometries.

posterior mean estimates are almost identical and uncertainty estimates differ only
slightly. However, the significance of the prior increases in the case of the limited
view sensor geometries. In addition, the simulations with different noise levels in
publication**I**indicate that proper prior information is important in the case of noisy
measurements.

In addition, the simulations show that sensor properties such as a size and band- width influence the estimates of the initial pressure and their uncertainties (Figure 4.5). In the case of the finite sized sensors, the quality of the estimates is reduced and the uncertainty is increased. This can be seen very clearly in the results ob- tained using the limited view sensor geometries (second and third column of Figure 4.5). The reduced quality cause significant increase in the relative errors (Table 4.1).

Also, the use of the bandlimited sensors reduces the accuracy of the solution, but the effect is not that significant as for the finite sized sensors.

**Figure 4.4:**The marginal probability densities of the posterior densities at locations
denoted by a square (first row) and circle (second row) in Figure 4.1 in the case
of the ideal sensors. The columns presents results obtained using the full view
(first column), two side (second column), and the one side (third column) sensor
geometry. Shown in the graphs are the true initial pressure (vertical black line)
together with the marginal densities of the posterior density obtained using the
Ornstein-Uhlenbeck prior (blue solid line) and white noise prior (red dotted line).

**Figure 4.5:** The posterior mean (top block) and standard deviation (bottom block)
obtained using the finite sized sensors (first row) and bandlimited sensors (second
row). The columns from left to right represent the full view (first column), two side
(second column) and one side (third column) sensor geometries. The red dots in the
images indicate the locations of the sensors.