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 backprojeceigenfunc-tion 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], 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 implemen-tation and they can have a low compuimplemen-tational 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 reverrever-sal such as iterative time reverrever-sal 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.