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 and nonzero depth dependent absorption coefficient μ(z), limited to z ≥ 0 and varying only along the z-direction, the initial acoustic stress response to photothermal heating takes the formTherein Γ, f0 and aB 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 , the scalar excess pressure field at time t and field point can be obtained by solving the inhomogeneous OA wave equation [4], [8]with c denoting the sonic speed in homogeneous media. Putting the dispersion relation of a harmonic wave solution with frequency ω and wave vector (satisfying the homogeneous wave equation) under scrutiny [35], it is possible to identify frequencies that correspond to solutions and that travel in positive and negative z-direction, 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 the transverse parameter yields a rational approximation that corresponds to a parabolic approximation of Eq. (2). Introducing time-retarded coordinates t → τ = t + zD/c the equivalent partial differential equation takes the form [4], [12], [35]Along propagation directions close to a solution to Eq. (3) yields a reasonable approximation to . In this paraxial approximation, the OA signal at a fixed field point along the beam axis can be related to the initial (t = 0) on-axis stress profile p0(τ) via a Volterra integral equation of the 2nd kind [12]Therein the Volterra operator features a convolution kernel K(τ − τ′) = ωD exp { − ωD(τ − τ′)}, 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]. The characteristic OA frequency effectively combines the defining parameters of the apparative setup psys ≡ (c, aB, zD). The acoustic near and far-field might be distinguished by means of the diffraction parameter , where near and far-field are characterized by D < 1 and D > 1, respectively. Subsequently we focus on OA signal detection in backward mode, i.e. zD < 0.

As an exemplary application, Fig. (1) illustrates the solution of the direct OA problem, i.e. the forward calculation of OA signals via Eq. (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].

Fig. 1

Illustration of OA signals obtained in the acoustic near-field (D = 0.22; red dashed curve) and far-field (D = 6.67; blue dash-dotted curve). The signals are obtained by solving the Volterra integral equation of the second kind, Eq. (4), that mediates the shape transformation of the initial acoustic stress profile (black solid curve) in the paraxial approximation. For the preparation of the initial stress profile it is assumed that μ = 5 mm−1 in the range [cτ−, cτ+] = [0, 0.3] mm and zero elsewhere. The convolution kernel reflects a Gaussian irradiation source profile with 1/e-width aB = 0.3 mm.

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 Fourier-expansion 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 byand Θ(·) signifies the Heavyside step-function. Then, for a suitable sequence a, the Fourier approximation to the Volterra integral equation, Eq. (4), readswith auxiliary functionsNow, consider a given set of input data (p0, pD) for known apparative parameters psys, both in a discretized setting with constant mesh interval Δ, mesh points where t0 = 0, t = t + Δ, and t large enough to ensure a reasonable measurement depth. Then, bearing in mind that τ = t + zD/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 Φℓ(τ ; R). In an attempt to construct an effective Volterra convolution kernel K(x ; a, R) for a controlled setup with a priori known parameters psys, one might use the high-precision “Gaussian-beam” estimator to obtain an initial sequence a0 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, aB 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 pD at psys = (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 μ = 24 cm−1 in the range z = 0 −0.1 cm, see black (p0) and blue (pD) curves in Fig. 2 (a). The set (p0, pD) is then used as inversion input to compute the effective Volterra kernel for various sets of reconstruction parameters prec = (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 prec = (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 prec = (51, 0.04 cm) exhibits a worse SSR.

Fig. 2

Kernel and source reconstruction within the paraxial approximation for system parameters psys = (c, aB, zD) ≡ (1 cm/s, 0.1 cm, − 0.5 cm). (a) Inversion input p0 (solid black line) and pD (solid blue line) used to derive effective kernel for N = 5, 11, and 51 Fourier-coefficients and cut-off parameter R = 0.06 cm. Solution of the respective source reconstruction problems yields the estimates pPL (dashed and dash-dotted red curves). (b) The main plot illustrates the effective kernel Keff(Δτ) ≡ K(Δτ ; a★, R) for two different cut-off distances R = 0.04 cm, and 0.06 cm. The inset shows the SSR s(R) ≡ s(a★, R) for N = 51 as function of the cut-off distance where the minimum is attained at R = 0.06 cm. (c) Solution pPL of the source reconstruction problem for a OA signal pD (solid blue line) resulting from a two-layer absorbing structure for the same system parameters as in (a). Source reconstruction is performed using the effective kernel for prec = (51, 0.06 cm) resulting from the gauge procedure.

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 psys, here indirectly accessible through the diffraction transformation of the OA signal pD relative to p0. That is, once the kernel reconstruction (I.2) is accomplished for a set of reference curves under psys, the source reconstruction problem (I.1) might subsequently be tackled also for all other OA signals measured under psys 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” , improved successively by solvingFrom a practical point of view we terminated the iterative correction scheme as soon as the max-norm of two successive solutions decreases below c ≤ 10−6. We here refer to the final estimate simply as pPL. 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 . 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 pD in the far-field approximation [11]. In contrast to this, low-precision predictors for both cases can be obtained by setting , where, e.g., c0 = 0.

The solution of the source reconstruction problem for the OA signal pD used in the approximation of the Volterra kernel for the above setting psys = (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 nmax taken by the Picard-Lindelöf correction scheme, with increasing number of expansion coefficients N, finding nmax ∝ N1.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 pPL for prec = (51, R★) and p0 does not come as a surprise since pD 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 μ = 24 cm−1 from z = 0 −0.05 cm and μ = 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 pPL in excellent agreement with the underlying exact initial stress profile p0.

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 sourcerecorded at the system parameters psys = (c, ρ0, aB, zD) = (1 cm/s, 0.1 cm, 0.1 cm, − 0.50 cm), and thus D = 2|zD|/(μa(aB + ρ0)) ≈ 1.04, obtained via an independent forward solver for the full OA wave equation designed for the solution of the OA Poisson integral for layered media [8], [10]. The inversion results are summarized in Fig. 3 (a), where the kernel reconstruction (inset) and source reconstruction (main plot) are shown for the parameter set prec = (41, R = 0.1 cm). 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 p0 and pPL suggests that the kernel reconstruction routine also applies to a more general OA setting, based on the full OA wave equation.

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 (p0, pD) to the inversion procedure (solid black and blue lines, respectively) as well as the reconstructed initial stress profile pPL (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 p0 (solid black line) used during the gauge procedure as well as the inversion input pE (orange line) for which the reconstructed initial stress profile pPL (dashed red line) is obtained. In both figures, the inset illustrates the effective Volterra kernel resulting from the Fourier-approximation.

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 psys = (150000 cm/s, 0.054 cm, 0.081cm/s, − 0.3 cm) as well as μ = 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 prec = (51, 0.1 cm) is shown in Fig. 3 (b). So as to perform the source reconstruction for the experimental signal pE, we considered data within the interval cτ = [0, 0.15] cm, only. As evident from the figure, the reconstructed stress profile pPL fits the signal p0 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 μ with finite value [40]. In the limit μ ≫ μ 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.