Optoacoustic inversion via convolution kernel reconstruction in the paraxial approximation and beyond

In this article we address the numeric inversion of optoacoustic signals to initial stress profiles. Therefore we study a Volterra integral equation of the second kind that describes the shape transformation of propagating stress waves in the paraxial approximation of the underlying wave-equation. Expanding the optoacoustic convolution kernel in terms of a Fourier-series, a best fit to a pair of observed near-field and far-field signals allows to obtain a sequence of expansion coefficients that describe a given “apparative” setup. The resulting effective kernel is used to solve the optoacoustic source reconstruction problem using a Picard-Lindelöf correction scheme. We verify the validity of the proposed inversion protocol for synthetic input signals and explore the feasibility of our approach to also account for the shape transformation of signals beyond the paraxial approximation including the inversion of experimental data stemming from measurements on melanin doped PVA hydrogel tissue phantoms.


Introduction
The inverse optoacoustic (OA) problem is concerned with the reconstruction of "internal" medium properties from "external" measurements of acoustic pressure signals. In contrast to the direct OA problem, referring to the calculation of a diffraction-transformed pressure signal at a desired field point for a given initial stress profile [1][2][3][4][5][6][7], one can distinguish two inverse OA problems: (I.1) the source reconstruction problem, where the aim is to invert measured OA signals to initial stress profiles upon knowledge of the mathematical model that mediates the underlying diffraction transformation [8,6,[9][10][11], and, (I.2) a kernel reconstruction problem, where the task is to find a convolution kernel that accounts for the apparent diffraction transformation shown by the OA signal. The latter arises quite naturally in a paraxial approximation wherein both signals can be related via a Volterra integral equation of the second kind [12]. Note that while problem I.1 is well established in the field of optoacoustics, we here make a first attempt at solving problem I.2, i.e. the kernel reconstruction problem, and demonstrate how it can be utilized for the reconstruction of initial stress profiles from observed OA signals.
Owing to its immediate relevance for medical applications [13][14][15][16][17][18][19], current progress in the field of inverse optoacoustics is spearheaded by OA tomography (OAT) and imaging applications in line with (I.1) [20][21][22][23][24][25], problem (I. 2) has not yet received much attention. However, quite similar kernel reconstruction problems are well studied in the context of inverse-scattering problems in quantum mechanics [26][27][28][29], and, from a more general point of view, are also studied in applied mathematics [30][31][32]. Note that analytic inversion formulae in OAT, aiming at reconstructing the full electromagnetic absorption distribution within the medium (see, e.g., Refs. [23][24][25]), assume that OA signals are detected from a full view, or, as in case of deconvolution reconstruction [33], still from a limited view of the object under scrutiny. If it is unfeasible to employ OAT techniques and one needs to resort to the inversion of data measured at a single point of the region of interest, due to either the inaccessibility of OAT inversion input or by other boundary conditions, kernel reconstruction in terms of (I.2) might provide an opportunity for OA inversion. However, note that the proposed approach does not evade the issue that reconstruction of data obtained from point detectors is, in general, not exact.
As a remedy, we here describe a numeric inversion scheme for problem (I.2), applicable to OA signals observed at a single field point, allowing to solve for a 1D absorption depth profile. Our aim is not to propose a competitive image reconstruction method for OAT applications that would require the reconstruction of full 3D domains from OA signals recorded at numerous detection angles. More precisely, in the presented article, we focus on the kernel reconstruction problem in the paraxial approximation to the optoacoustic wave-equation, where we suggest a Fourier-expansion approach to construct an approximate optoacoustic convolution kernel. We show that once (I.2) is solved for a given "apparative" setup, this then allows to subsequently solve (I.1) for different signals obtained using an identical apparative setup. A central and reasonable assumption of our approach is that the influence of the stress wave propagator on the shape change of the OA signal is negligible above a certain cut-off distance. After developing and testing the numerical procedure in the paraxial approximation, we assess how well the inversion protocol carries over to more prevalent optoacoustic problem instances, featuring the reconstruction for: (i) the full OA wave-equation, (ii) non Gaussian irradiation source profiles, and, (iii) measured signals exhibiting noise.
The article is organized as follows: in Section 2 we recap the theoretical framework of OA signal generation in the paraxial approximation, in Section 3 we discuss our approach to the OA kernel reconstruction problem, in Section 4 we illustrate the subsequent solution of the source reconstruction problem via the obtained approximate convolution kernel, and allude to the challenging problem of OA signal inversion beyond the paraxial approximation in Section 5. In Section 6, we conclude with a summary.

The direct OA problem in the paraxial approximation
The dominant microscopic mechanism contributing to the generation of acoustic stress waves is expansion due to photothermal heating [34]. In the remainder we assume a pulsed photothermal source with pulse duration short enough to ensure thermal and stress confinement [8]. Then, in case of a purely absorbing material exposed to an irradiation source profile with beam axis along the z-direction of an associated coordinate system, a Gaussian profile in the transverse coordinates r and nonzero depth dependent absorption coefficient μ a (z), limited to z ≥ 0 and varying only along the z-direction, the initial acoustic stress response to photothermal heating takes the form Therein Γ, f 0 and a B signify the Grüneisen parameter, the intensity of the irradiation source along the beam axis and the 1/e-width of the beam profile orthogonal to the beam axis, respectively. Given the above initial instantaneous acoustic stress field p r ( ) 0 , the scalar excess pressure field p r t ( , ) at time t and field point r can be obtained by solving the inhomogeneous OA wave equation [4,8]  x y z (satisfying the homogeneous wave equation) under scrutiny [35], it is possible to identify frequencies that correspond to solutions + p r t ( , ) and p r t ( , ) that travel in positive and negative zdirection, respectively. This allows to derive a more simple propagation equation than Eq. (2) that forms an adequate approximation to p + , only. I.e., employing a first order expansion of the dispersion relation in that corresponds to a parabolic approximation of Eq. (2). Introducing time-retarded coordinates t → τ = t + z D /c the equivalent partial differential equation takes the form [4,12,35] Along propagation directions close to e z a solution p r t ( , ) to Eq. (3) yields a reasonable approximation to Therein the Volterra operator features a convolution kernel K mediating the diffraction transformation of the propagating stress waves [12]. A derivation of the exponential convolution kernel in terms of the transfer function method working in the spectral domain is detailed in Ref. [4]. for a problem setup that resembles the experimental setup reported in Ref. [36]. For a comparison of the prediction of OA signals in terms of the effectively 1D approach provided by Eq. (4) as opposed to the full 3D wave equation please refer to appendix B of Ref. [7].

Reconstruction of the OA convolution kernel
Note that the solution of the direct problem and inverse problem (I.1) in terms of Eq. (4) is feasible using standard numerical schemes based on, e.g., a trapezoidal approximation of the Volterra operator for a generic kernel [37], or highly efficient inversion schemes for the particular form of the above convolution kernel [11]. As pointed out earlier, considering inverse problem (I.2), we here suggest a Fourierexpansion of the convolution kernel involving a sequence of N expansion coefficients and a cut-off distance R above which the resulting effective kernel is assumed to be zero, i.e.
The expansion functions k ℓ (x ; R) are given by and Θ(·) signifies the Heavyside step-function. Then, for a suitable sequence a, the Fourier approximation to the Volterra integral equation, Eq. (4), reads with auxiliary functions Now, consider a given set of input data (p 0 , p D ) for known apparative parameters p sys , both in a discretized setting with constant mesh interval Δ, mesh points t where t 0 = 0, t i = t i−1 + Δ, and t M large enough to ensure a reasonable measurement depth. Then, bearing in mind that τ i = t i + z D /c, the optimal expansion coefficient sequence a ★ can be obtained by minimizing the sum of the squared residuals (SSR) In the above optimization formulation of inverse problem (I.2), we considered a trapezoidal rule to numerically evaluate the integrals that enter via the functions Φ ℓ (τ i ; R). In an attempt to construct an effective Volterra convolution kernel K(x ; a, R) for a controlled setup with a priori known parameters p sys , one might use the high-precision "Gaussian-beam" esti- to obtain an initial sequence a 0 of expansion coefficients by means of which a least-squares routine for the minimization of Eq. (9) might be started. In a situation where, say, a B is only known approximately or the assumption of a Gaussian beam profile is violated, one has to rely on a rather low-precision coefficient estimate obtained by roughly estimating the apparative parameters and resorting on the above "Gaussian-beam" estimate.
An exemplary kernel reconstruction procedure is shown in Fig. 2 , where the OA signal p D at p sys = (1 cm/s, 0.1 cm, − 0.5 cm), i.e. D ≈ 3.75, is first obtained by solving the direct OA problem for Eq. (4) for an absorbing layer with μ a = 24 cm −1 in the range z = 0 −0.1 cm, see black (p 0 ) and blue (p D ) curves in Fig. 2 (a). The set (p 0 , p D ) is then used as inversion input to compute the effective Volterra kernel for various sets of reconstruction parameters p rec = (N, R). In particular, considering N = 51, the minimal value of s(a ★ , R ★ ) ≈ 1.47 is attained at R ★ = 0.06 cm, see the inset of Fig. 2 (b). As evident from the main plot of Fig. 2 (b), the effective Volterra kernel for p rec = (51, R ★ ) follows the exact stress wave propagator for almost two orders of magnitude up to cΔτ ≈ 0.05 cm. Beyond that limit, the noticeable deviation between both does not seem to affect the overall SSR s(a, R) too much. In this regard, note that the kernel approximated for the (non optimal) choice p rec = (51, 0.04 cm) exhibits a worse SSR.

The inverse OA problem -source profile reconstruction
Note that the above Fourier-expansion approximation might be interpreted as a gauge procedure to adjust an effective Volterra kernel K(x ; a ★ , R) for an (possibly unknown) apparative setup p sys , here indirectly accessible through the diffraction transformation of the OA signal p D relative to p 0 . That is, once the kernel reconstruction (I.2) is accomplished for a set of reference curves p p ( , ) 0 D ref under p sys , the source reconstruction problem (I.1) might subsequently be tackled also for all other OA signals measured under p sys by solving the OA Volterra integral equation Eq. (4) in terms of a Picard-Lindelöf "correction" scheme [38]. The latter is based on the continued refinement of a putative solution, starting off from a properly guessed "predictor" p ( ) From a practical point of view we terminated the iterative correction scheme as soon as the max-norm of two successive solutions decreases below c n ≤ 10 −6 . We here refer to the final estimate simply as p PL . Note that, attempting a solution of (I.1) in the acoustic near-field, a high-precision predictor can be obtained by using the initial guess p p PL (0) D . This is a reasonable choice since one might expect the change of the OA near-field signal due to diffraction to be still quite small. Further, source reconstruction in the acoustic far-field might be started using a high-precision predictor obtained by integrating the OA signal p D in the far-field approximation [11]. In contrast to this, low-precision predictors for both cases can be obtained by setting p c PL (0) 0 , where, e.g., c 0 = 0. The solution of the source reconstruction problem for the OA signal p D used in the approximation of the Volterra kernel for the above setting p sys = (1 cm/s, 0.1 cm, − 0.5 cm) is shown in Fig. 2 (a). We assessed the scaling of the speed of convergence, measured using the number of steps n max taken by the Picard-Lindelöf correction scheme, with increasing number of expansion coefficients N, finding n max ∝ N 1.3 . Note, however, that the computational burden of the Picard-Lindelöf correction scheme is inferior to the minimization of the SSR according to Eq. (9). The apparent agreement of the data curves p PL for p rec = (51, R ★ ) and p 0 does not come as a surprise since p D was used for the gauge procedure in the first place. As a remedy we attempt a source reconstruction for a second independent OA signal, simulated for the same apparative setting only with two absorbing layers μ a,1 = 24 cm −1 from z = 0 −0.05 cm and μ a,2 = 12 cm −1 from z = 0.05 − 0.12 cm. As evident from Fig. 2 (c), inversion using the effective Volterra kernel from the previous gauge procedure yields a reconstructed stress profile p PL in excellent agreement with the underlying exact initial stress profile p 0 .

Inversion beyond the paraxial approximation
Given the apparent feasibility of the kernel reconstruction routine as a gauge procedure to model the diffraction transformation of OA signals in terms of an effective stress wave propagator in the framework of the OA Volterra integral equation, we next address the inversion of OA signals to initial stress profiles beyond the paraxial approximation. Therefore, we first consider a borderline far-field signal for a top-hat irradiation source The position of the peak of the effective kernel seems due to the finite extension ρ 0 of the employed top-hat beam profile (bear in mind that originally, the underlying 1D approach assumes a Gaussian beam profile). We find that, the larger the respective top-hat width ρ 0 , the further out that peak occurs. Consequently, the cut-off distance R above which the kernel is assumed to vanish needs to be chosen large enough to still enclose the peak of the kernel. The excellent agreement of the stress profiles p 0 and p PL suggests that the kernel reconstruction routine also applies to a more general OA setting, based on the full OA wave equation. Finally, we consider an OA signal resulting from an actual measurement on PVA hydrogel based tissue phantoms [10]. In this case we carefully estimated the apparative parameters p sys = (150000 cm/s, 0.054 cm, 0.081cm/s, − 0.3 cm) as well as μ a = 11 cm −1 in the range z = 0 −0.095 cm, i.e. D ≈ 6.73, in order to create a set of synthetic input data by means of which an appropriate kernel gauge procedure can be carried out. The result of the procedure using p rec = (51, 0.1 cm) is shown in Fig. 3 (b). So as to perform the source reconstruction for the experimental signal p E , we considered data within the interval cτ = [0, 0.15] cm, only. As evident from the figure, the reconstructed stress profile p PL fits the signal p 0 used in the gauge procedure remarkably well.

Conclusions
In the presented article we have introduced and discussed the kernel reconstruction problem in the paraxial approximation of the optoacoustic wave equation for both, synthetic input data and experimental data resulting from controlled measurements on melanin doped PVA hydrogel tissue phantoms. We suggested a Fourier-expansion approach to approximate the convolution kernel which takes a central role in the theoretical framework. The developed approach proved useful as gauge procedure by means of which the diffraction transformation experienced by OA signals can effectively be modeled, allowing to subsequently solve the source reconstruction problem in the underlying apparative setting. From this numerical study we found that the developed approach extends beyond the framework of the paraxial approximation and also allows for the inversion of OA signals described by the full OA wave equation in the acoustic far field. It would be tempting to explore other kernel expansions in terms of generalized Fourier series and to assess the presented method with regard to different signal-to-noise ratios in the input data. Also, the effects of acoustic attenuation, the impulse response of the employed transducer and uncertainty in the system parameters might be studied in detail so as to probe the limits of the proposed inversion scheme. Such investigations are currently in progress with the aim to shed some more light on this intriguing inverse problem in the field of optoacoustics and to facilitate a complementary approach to conventional OA signal inversion.
A Python implementation of our research-code for the solution of inverse problems (I.1) and (I.2), along with all scripts needed to reproduce the figures is publicly available on one of the authors figshare profile, see Ref. [39].
Finally, note that here we discussed the problem of optoacoustic inversion in the limit of unscattered transmission. However, in general, the propagation of light in biological tissue is goverened by a scattering coefficient μ s with finite value [40]. In the limit μ s ≫ μ a this renders the transmission process effectively diffusive. As a consequence, the light beam may behave as Gaussian only for depths > 1 mm. In this depth range, optical-resolution photoacoustic microscopy (OR-PAM) [41] might be of interest.

Conflict of interests
None. Fig. 3. Inversion of OA signals to initial stress profiles beyond the paraxial approximation. Both figures illustrate the kernel and source reconstruction procedures for (a) inversion of an OA signal featuring a top-hat irradiation source profile (see text). The main plot shows the input (p 0 , p D ) to the inversion procedure (solid black and blue lines, respectively) as well as the reconstructed initial stress profile p PL (dashed red line), and, (b) inversion of an OA signal resulting from an actual measurement [10]. The main plot shows the synthetic initial stress profile p 0 (solid black line) used during the gauge procedure as well as the inversion input p E (orange line) for which the reconstructed initial stress profile p PL (dashed red line) is obtained. In both figures, the inset illustrates the effective Volterra kernel resulting from the Fourier-approximation.