Effective one-dimensional approach to the source reconstruction problem of three-dimensional inverse optoacoustics

The direct problem of optoacoustic signal generation in biological media consists of solving an inhomogeneous three-dimensional (3D) wave equation for an initial acoustic stress profile. In contrast, the more defiant inverse problem requires the reconstruction of the initial stress profile from a proper set of observed signals. In this article, we consider an effectively 1D approach, based on the assumption of a Gaussian transverse irradiation source profile and plane acoustic waves, in which the effects of acoustic diffraction are described in terms of a linear integral equation. The respective inverse problem along the beam axis can be cast into a Volterra integral equation of the second kind for which we explore here efficient numerical schemes in order to reconstruct initial stress profiles from observed signals, constituting a methodical progress of computational aspects of optoacoustics. In this regard, we explore the validity as well as the limits of the inversion scheme via numerical experiments, with parameters geared toward actual optoacoustic problem instances. The considered inversion input consists of synthetic data, obtained in terms of the effectively 1D approach, and, more generally, a solution of the 3D optoacoustic wave equation. Finally, we also analyze the effect of noise and different detector-to-sample distances on the optoacoustic signal and the reconstructed pressure profiles.

Figure 1

Illustration of an OA setup with plane-normal irradiation source profile $f\left(\right|{\vec{r}}_{\perp }\left|\right)$ and a source volume consisting of possibly multiple stacked absorbing layers. The irradiation source exhibits a $1/e$-intensity radius of $|{\vec{r}}_{\perp }|=a_0$, and the different layers possess absorption coefficients ${\mu }_a$ as indicated in the figure. Here, the figure shows a single absorbing layer of width $\mathrm{\Delta }z=0.1\mathrm{cm}$. The initial stress field $\propto p_0\left(\vec{r}\right)$ causes acoustic pressure waves that can be monitored as an OA signal at the detection point $\mathrm{D}$, here located on the beam axis, i.e., at $|{\vec{r}}_{\perp }|=0$, with distance $z_{\mathrm{D}}<0$ from the absorbing layer.

Figure 2

Forward calculation of OA signals $p\left(\tau \right)$ in the framework of the Volterra integral equation. In this figure, $p_0$ (black solid lines) indicates the on-axis profile of the initial stress according to Eq. (4). The figure illustrates the change in shape of the OA signals as perceived at detection points with increasing detector-to-layer distance $z_{\mathrm{D}}=-0.02\mathrm{cm}$ (magenta dotted line), $-0.1\mathrm{cm}$ (blue dash-dotted line), and $-1\mathrm{cm}$ (red dashed line), characterized by the diffraction parameters $D=0.17$, 1.0, and, 8.33, respectively.

Figure 3

Forward calculation of OA signals in a far-field approximation considering Eq. (11). In this figure, $p_0$ (black solid lines) indicates the on-axis profile of the initial acoustic stress. The figure illustrates the difference in the approximate FF formula, resulting in the signal $p_{\mathrm{FH}}$ (red dashed curve), relative to the exact solution $p$ (red dash-dotted curve) at $z_{\mathrm{D}}=-1\mathrm{cm}$, i.e., $D=8.33$.

Figure 4

Reconstruction of the initial acoustic stress ${\tilde{p}}_0$ from computed OA signals $p$ upon knowledge of the diffraction propagator $\mathcal{K}$ in the acoustic NF (${\tilde{p}}_{0,\mathrm{NF}}$) at $z_{\mathrm{D}}=-0.02\mathrm{cm}$, i.e., $D=0.17$, and the acoustic FF (${\tilde{p}}_{0,\mathrm{FF}}$) at $z_{\mathrm{D}}=-1.0\mathrm{cm}$, i.e., $D=8.33$. The OA source reconstruction is accomplished with the leap-frog algorithm introduced in Sec. 3b. In either case, the reconstructed stress profiles perfectly match the true initial stress profiles $p_0$.

Figure 5

Solution of the source reconstruction problem using the Volterra integral equation based inverse solver. The figure illustrates the reconstruction of the initial acoustic stress profile ${\tilde{p}}_0$ (dashed red line) from a FF signal $p$ (dash-dotted blue line), calculated at $z_{\mathrm{D}}=-1\mathrm{cm}$. The true initial stress profile is given by $p_0$ (solid black line). Reconstruction starting from a genuine OA signal $p$ that was superimposed by Gaussian white noise with a signal-to-noise ratio of 5.

Figure 6

Root-mean-square error (RMSE) between the true initial stress profile $p_0$ and the inverse estimate ${\tilde{p}}_0$ for increasingly narrow intervals $\mathrm{\Delta }z\equiv [c{\tau }_{-}:c{\tau }_{+}]$ along the $z$ and $\tau$ axes, respectively. Note that in a preprocessing step, $p_0$ was normalized and ${\tilde{p}}_0$ amplitude-adjusted, as detailed in the text.

Figure 7

Effect of a finite detector size on the OA signals and the reconstructed initial stress profiles. For the simulations, we considered detectors with a circular surface and different radii in the range $R_{\mathrm{D}}=0.1-4\mathrm{mm}$ (for further details, see Sec. 4c). (a) OA signals obtained by solving the Poisson-like integral Eq. (2) and integration over the detector surface. (b) Reconstructed initial stress profiles obtained by inversion via Eq. (10). Note the exact initial stress profile indicated by the black dashed line. (c) Reconstructed absorption coefficient profile ${\mu }_{\mathrm{rec}}\left(z\right)$. The exact absorption coefficient profile used to obtain the curves in (a) is indicated by the black dashed line. In either subfigure, the blue solid line emphasizes the curves for a detector radius $R_{\mathrm{D}}=0.5\mathrm{mm}$.

Figure 8

Comparison of two inversion procedures derived from two different linear integral equations describing the diffraction transformation of OA signals. The (fast) numerical procedure based on the inverse solution of a Volterra integral equation of the second kind is indicated by the blue dashed line (labeled V-inversion), and the procedure based on the solution of an inhomogeneous Fredholm equation of the first kind (here accomplished using Wiener deconvolution) is indicated by the red dotted line (labeled F-inversion); for further details see Sec. 5.