Image reconstruction and modelling of uncertainties in photoacoustic tomography

(1)

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

(2)

(3)

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

(4)

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

(5)

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

(6)

(7)

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 biomedical 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 problem 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 information 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 uncertainties 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; inverse problems; modelling; simulation; errors; error correction; Bayes methods; parameter estimation; image reconstruction; ultrasonic velocity; uncertainty handling

Yleinen suomalainen ontologia:kuvantaminen; fotoakustinen kuvantaminen; tomografia;

(8)

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

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

(9)

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 thesis.

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

(10)

(11)

LIST OF PUBLICATIONS

This thesis consist of the present review of author’s work in the field of photoacoustic 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 America139(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 America144(4), 2061–2071 (2018).

III J. Leskinen, A. Pulkkinen, J. Tick, and T. Tarvainen, ”Photoacoustic tomography 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 photoacoustic 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 author 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- cationsI,IIandIV and the co-writer of publicationIII. The author participated in the measurements in publication III. The photoacoustic tomography measurement data in publicationIIoriginated from co-author Robert Ellwood.

(12)

(13)

1 Introduction

In biomedical imaging, one is interested in obtaining images of the internal structure 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 acoustic 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 require reconstruction algorithms to obtain PA images whereas in PAM, a PA image is directly obtained by mechanically scanning. Depending on whether acoustic focusing 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 expansion. 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 distribution 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 resolution 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 imaging 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

(18)

the initial pressure distribution of an imaged target based on the boundary measurements 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 problem, 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 measurements, 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 compensate for errors caused by an unknown speed of sound in PAT.

3. Evaluate the feasibility of the developed methods with simulations and experimental 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.

(19)

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 information 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 target. 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 implemented 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

(20)

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 absorbed energy is converted into heat and the temperature of the target increases slightly, which temperature increase is given by

T(r) = ^H(r)

ρCV , (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) ρCVκ

= ^β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 detection geometry due to their larger apertures. The sensor geometries can be implemented 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].

(21)

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 distribution 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 amplitude. 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









 ∂²

∂t² −c(r)²∇²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₀(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 computationally 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 images, 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 almost 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 estimates 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

(22)

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 inverse 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], Iand 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 implementation and they can have a low computational cost. However, they are limited to specific geometries such as spherical, cylindrical, and planar acoustic detection surfaces 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 approach, 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.

(23)

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 reversal 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 setups [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 ultrasound 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 modelled. 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 applications 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 approximate speed of sound that inaccurately addresses the true speed of sound distribution 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.

(24)

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 vasculature. In addition, a structure of vessels can be visualised using intravascular PA imaging catheters [130–133]. Since increased vascularity, high total blood concentrations 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 produce 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 tissue [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 addition, PAI have been utilised in diagnosing joint diseases [170–172] and Crohn disease [173].

(25)

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].

(26)

(27)

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 publicationII. Finally, the Bayesian approximation error approach and its utilisation in modelling uncertainties in the speed of sound is described following publicationIV.

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^Mis 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^Mrepresents 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₀) is the prior probability density and π(pt|p₀) is the likelihood density. 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)

(28)

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 distributed i.e. e ∼ N(ηe,Γe) and p0 ∼ N(ηp₀,Γp₀) where ηe and Γe are the mean and covariance of the noise, respectively, and ηp₀ and Γp₀ are the mean and covariance 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₀_|p_t,Γ_p₀|p_t)and is described by the mean and covariance

η_p₀_|p_t = _Γ_p

0|pt

K(c)^T_Γ⁻¹_e (pt−ηe) +_Γ⁻¹_p

0ηp₀

(3.5) Γ_p₀_|p_t = K(c)^TΓ⁻¹_e K(c) +Γ⁻¹_p₀−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 distributions, the posterior density can in principle be explored by employing Markov Chain Monte Carlo (MCMC) technique, which provides a representative set of samples 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 typically computed. The most commonly used point estimate is amaximum a posteriori (MAP) estimate that is defined as

p_0,MAP=arg max

p₀

π(p₀|pt). (3.7)

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

p_0,MAP=η_p₀_|_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 information can be used to assess the accuracy of the estimate. From the posterior covariance posterior variances σ_p²

0|p_t = diag{Γ_p₀|p_t}and further standard deviations σ_p₀_|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₀_|p_t_,k−ασ_p₀_|p_t_,k, η_p₀_|_p_t_,k+ασ_p₀_|p_t_,k

i, (3.8)

where η_p₀_|p_t_,k is the value of the posterior meanη_p₀_|p_t in thekth element,σ_p₀_|p_t_,k is the value of σ_p₀_|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

(29)

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 exponential (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 exponential 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₀ and covariance matrices Γp₀. If the range of the unknown parameter values is roughly known, this information can be utilised to choose the meanηp₀ of the prior. Corre- spondingly, information related to the expected spatial distribution of the parameters can be utilised to choose the covariance matrixΓp₀. 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₀ is defined as

Γp₀,ij =σ_p²₀δ_ij, (3.10)

where i and j are the element indices, σp₀ 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₀. 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,

(30)

theΓp₀ is of form

Γp₀,ij=σ_p²₀exp −kr_i−r_jk² 2l²

!

(3.11) whereas in the case of the Ornstein-Uhlenbeck prior theΓp₀ is defined as

Γp₀,ij=σ²_p₀exp

−kr_i−r_jk 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 addition, 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 using iterative linear solvers, such as conjugate gradient methods, and evaluating the multiplication by the forward modelKwith the forward operator [96], the multiplication by the transpose of the forward modelK^T with the adjoint operator [45] and the multiplication by the covariance of the prior Γp₀ 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₀|p_t =d, (3.13)

where

C=_Γ_p₀K(c)^T_Γ⁻¹_e K(c) +I, (3.14) d=Γp₀K(c)^TΓ⁻¹_e (pt−η_e) +η_p₀ (3.15) and Iis an identity matrix. Correspondingly, thekth column of the posterior covariance (3.6) can be computed by solving the linear system

CΓ_p₀|pt,k =_Γ_p₀e_k, (3.16) whereCis as in (3.14) ande_kis 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.

(31)

3.2 BAYESIAN APPROXIMATION ERROR APPROACH FOR MOD- ELLING OF ERRORS

A forward model that describes the dependence of unknown parameters and measurements 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 measurement 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 approximation error modelling has been used to model and compensate errors caused by uncertainties in the speed of soundIV(Section 3.2.1) and ultrasound sensor locations [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(_ec)p0+e+e

= K(_ec)p0+n, (3.17)

where e = K(c)p0−K(_ec)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 unknown p0, the observation model (3.17) leads to a likelihood density

π(pt|p0) =πn(pt−K(_ec)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(_ec)p0). (3.19)

(32)

Let us assume that noisee is Gaussian distributed e ∼ N(ηe,Γe), and approximate 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+η_e and Γn = Γ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₀_|p_t,Γ_p₀|p_t)given by the mean and covariance

η_p₀_|p_t = _Γ_p₀_|p_tK(_ec)^T_Γ⁻¹_n (pt−ηn) +_Γ⁻¹_p₀ηp₀

(3.20)

Γ_p₀|p_t = _K(_e_c)^T_Γ⁻¹_n _K(_e_c) +_Γ⁻¹_p

0

−1

. (3.21)

The realisation of the approximation erroreis unknown since its value depends on the actual unknown p₀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)₀ ,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)₀ −K(_ec)p₀^(l). (3.22) From these samples, the meanη_eand covarianceΓeof the approximation error can be estimated as

η_e= ¹ Ns

N_s l=1

∑

e^(l) (3.23)

Γe = ¹ Ns−1

Ns

l=1

∑

(e^(l)−η_e)(e^(l)−η_e)^T. (3.24)

(33)

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 publicationsIandII. Secondly, results of the experimental validation of the approach are presented in Section 4.2. These results are based on publications IIand 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 publicationIV.

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 publicationI, simulations were done in 2D whereas in publicationII simulations were conducted in 3D. In the 2D simulations, the approach was investigated 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 posterior 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 evaluated by computing the relative errors of the estimates with respect to the true initial pressure distribution using

Ep₀ =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 sensors (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 bandlimited 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 simulated measurement data was generated by solving the wave equation (2.4) using the

(34)

Figure 4.1: The simulated (true) initial pressure distribution given in arbitrary units. The white square and circle indicate the locations where the marginal densities 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 simulation, 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) using (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₀ = 5 and the standard deviation was set asσ_p₀ = 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 estimates 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 estimates 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 distribution 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

(35)

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 sensors (ideal), finite sized sensors (finite sized), and bandlimited sensors (bandlimited) 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

(36)

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 publicationIindicate 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 bandwidth 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 obtained 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.

(37)

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).

(38)

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.

Image reconstruction and modelling of uncertainties in photoacoustic tomography

uef.fi

Dissertations in Forestry and Natural Sciences

JENNI TICK

IMAGE RECONSTRUCTION AND MODELLING OF UNCERTAINTIES IN PHOTOACOUSTIC TOMOGRAPHY

JENNI TICK

Jenni Tick

IMAGE RECONSTRUCTION AND MODELLING OF UNCERTAINTIES IN

PHOTOACOUSTIC TOMOGRAPHY

TABLE OF CONTENTS

1 Introduction

2 Photoacoustic tomography

3 Bayesian inversion methods in PAT

∑

∑

4 Summary of the results