Optoacoustic inversion via Volterra kernel reconstruction
aa r X i v : . [ phy s i c s . c o m p - ph ] J un Optoacoustic inversion via Volterra kernel reconstruction
O. Melchert, ∗ M. Wollweber, and B. Roth
Hannover Centre for Optical Technologies (HOT),Leibniz Universit¨at Hannover, D-30167 Hannover, Germany (Dated: September 3, 2018)In this letter we address the numeric inversion of optoacoustic signals to initial stress profiles.Therefore we put under scrutiny the optoacoustic kernel reconstruction problem in the paraxialapproximation of the underlying wave-equation. We apply a Fourier-series expansion of the op-toacoustic Volterra kernel and obtain the respective expansion coefficients for a given “apparative”setup by performing a gauge procedure using synthetic input data. The resulting effective kernel issubsequently used to solve the optoacoustic source reconstruction problem for general signals. Weverify the validity of the proposed inversion protocol for synthetic signals and explore the feasibilityof our approach to also account for the diffraction transformation of signals beyond the paraxialapproximation.
The inverse optoacoustic (OA) problem is concernedwith the reconstruction of “internal” OA properties from“external” measurements of acoustic pressure signals. Incontrast to the direct
OA problem, referring to the cal-culation of a diffraction-transformed pressure signal at adesired field point for a given initial stress profile [1–4],one can distinguish two inverse OA problems: (I.1) the source reconstruction problem , where the aim is to in-vert measured OA signals to initial stress profiles uponknowledge of the mathematical model that mediates theunderlying diffraction transformation [4–6], and, (I.2) the kernel reconstruction problem , where the task is to recon-struct a proper OA stress-wave propagator to accountfor the apparent diffraction transformation shown by theOA signal. While, owing to its immediate relevance formedical applications [7], current progress in the field ofinverse optoacoustics is spearheaded by OA tomographyand imaging applications in line with (I.1) [8, 9], problem(I.2) has not yet received much attention (note that quitesimilar kernel reconstruction problems are well studiedin the context of inverse-scattering problems in quan-tum mechanics [10]). However, under ill-conditioned cir-cumstances that prohibit a consistent description of thestress-wave propagation or when the multitude of signalsthat form the inversion input to common backpropaga-tion approaches (see, e.g., Refs. [9]) are simply inacces-sible, kernel reconstruction in terms of (I.2) provides anopportunity to yield a reliable OA inversion protocol interms of single-shot measurements.As a remedy, we here describe a numerical approachto problem (I.2), appealing from a point of view of com-putational theoretical physics. More precisely, in thepresented letter, we focus on the kernel reconstructionproblem in the paraxial approximation to the optoacous-tic wave-equation, where we suggest a Fourier-expansionapproach to construct an approximate stress wave propa-gator. We show that once (I.2) is solved for a given “ap-parative” setup, this then allows to subsequently solve(I.1) for different signals obtained using an identical ap-parative setup. A central and reasonable assumption of our approach is that the influence of the stress wave prop-agator on the shape change of the OA signal is negligibleabove a certain cut-off distance. After developing andtesting the numerical procedure in the paraxial approxi-mation, we assess how well the inversion protocol carriesover 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 direct OA problem.
The dominant microscopicmechanism contributing to the generation of acousticstress waves is expansion due to photothermal heating[11]. In the remainder we assume a pulsed photothermalsource with pulse duration short enough to ensure ther-mal and stress confinement [5]. Then, in case of a purelyabsorbing material exposed to a irradiation source pro-file with beam axis along the z -direction of an associatedcoordinate system, a Gaussian profile in the transversecoordinates ~r ⊥ and nonzero depth dependent absorptioncoefficient µ a ( z ), limited to z ≥ z -direction, the initial acoustic stress response to pho-tothermal heating takes the form p ( ~r ) = f µ a ( z ) exp n − | ~r ⊥ | /a − Z z µ a ( z ′ ) d z ′ o . (1)Therein f and a B signify the intensity of the irradiationsource along the beam axis and the 1 /e -width of the beamprofile orthogonal to the beam axis, respectively. Giventhe above initial instantaneous acoustic stress field p ( ~r ),the scalar excess pressure field p ( ~r, t ) at time t and fieldpoint ~r can be obtained by solving the inhomogeneousOA wave equation [2, 5] (cid:2) ∂ t − c ∆ (cid:3) p ( ~r, t ) = p ( ~r ) ∂ t δ ( t ) , (2)with c denoting the sonic speed within the medium. Theacoustic near and far-field might be distinguished bymeans of the diffraction parameter D = 2 | z D | / ( µ a a ),where near and far-field are characterized by D <
D >
1, respectively.In the paraxial approximation where the full waveequation reduces to the parabolic diffraction equation[ ∂ τ ∂ z − ( c/ ⊥ ] p = 0 [2, 12], it can be shown that thetime-retarded ( τ = t + z D /c ) OA signal at a field pointalong the beam axis p D ( τ ) ≡ p ( ~r D , t ) can be related to theinitial ( t = 0) on-axis stress profile p ( τ ) ≡ p ( ~r ⊥ = 0 , z )via a Volterra integral equation of 2nd kind, reading [12] p D ( τ ) = p ( τ ) − Z τ −∞ K ( τ − τ ′ ) p ( τ ′ ) d τ ′ . (3)Therein the Volterra operator features a convolutionkernel K ( τ − τ ′ ) = ω D exp {− ω D ( τ − τ ′ ) } , mediating thediffraction transformation of the propagating stresswaves. The characteristic OA frequency ω D = 2 c | z D | /a effectively combines the defining parameters of the appa-rative setup p sys ≡ ( c, a B , z D ). Subsequently we focus onOA signal detection in backward mode, i.e. z D < The inverse OA kernel reconstruction problem.
Notethat the solution of the direct problem and inverse prob-lem (I.1) in terms of Eq. (3) is feasible using standard nu-merical schemes based on, e.g., a trapezoidal approxima-tion of the Volterra operator for a generic kernel [13], orhighly efficient memoization techniques for the particularform of the above convolution kernel [14]. As pointed outearlier, considering inverse problem (I.2), we here suggesta Fourier-expansion of the Volterra kernel involving a se-quence of N expansion coefficients a ≡ { a ℓ } ≤ ℓ 0 0.05 0.1 0.15 c τ (cm) p ( a r b . un it s ) z D = -0.50 cm R = 0.06 cm(a) p PL ( τ ) - N =51 p PL ( τ ) - N =11 p PL ( τ ) - N =5 p D ( τ ) p ( τ )10 -2 -1 0 0.025 0.05 c ∆τ (cm) K ( s ) N = 51(b) K ( ∆τ ) K eff ( ∆τ ) - R = 0.04 cm K eff ( ∆τ ) - R = 0.06 cm Rs(R) 0 0.05 0.1 0.15 c τ (cm) p ( a r b . un it s ) z D = -0.50 cm R = 0.06 cm (c) p PL ( τ ) - N =51 p D ( τ ) p ( τ ) FIG. 1. (Color online) Kernel and source reconstructionwithin the paraxial approximation for system parameters p sys = ( c, a B , z D ) ≡ (1 cm / s , . , − . p (solid black line) and p D (solid blue line) used to de-rive effective kernel for N = 5, 11, and 51 Fourier-coefficientsand cut-off parameter R = 0 . 06 cm. Solution of the respec-tive source reconstruction problems yields the estimates p PL (dashed and dash-dotted red curves). (b) The main plot il-lustrates the effective kernel K eff (∆ τ ) ≡ K (∆ τ ; a ⋆ , R ) for twodifferent cut-off distances R = 0 . 04 cm, and 0 . 06 cm. Theinset shows the SSR s ( R ) ≡ s ( a ⋆ , R ) for N = 51 as func-tion of the cut-off distance where the minimum is attained at R = 0 . 06 cm. (c) Solution p PL of the source reconstructionproblem for a OA signal p D (solid blue line) resulting from atwo-layer absorbing structure for the same system parametersas in (a). Source reconstruction is performed using the effec-tive kernel for p rec = (51 , . 06 cm) resulting from the gaugeprocedure. ensure a reasonable measurement depth. Then, bearingin mind that τ i = t i + z D /c , the optimal expansion co-efficient sequence a ⋆ can be obtained by minimizing thesum of the squared residuals (SSR) s ( a , R ) = M X i =0 h ( p ( τ i ) − p D ( τ i )) − N − X ℓ =0 a ℓ Φ ℓ ( τ i ; R ) i . (8)In the above optimization formulation of inverse prob-lem (I.2), we considered a trapezoidal rule to numer-ically evaluate the integrals that enter via the func-tions Φ ℓ ( τ i ; R ). In an attempt to construct an effec-tive Volterra kernel K ( x ; a , R ) for a controlled setupwith a priori known parameters p sys , one might usethe high-precision “Gaussian-beam” estimator a ℓ =(2 ω D /R ) R R k ℓ ( x ; R ) exp {− ω D x } d x to obtain an initialsequence a of expansion coefficients by means of whicha least-squares routine for the minimization of Eq. (8)might be started. In a situation where, say, a B is onlyknown approximately or the assumption of a Gaussianbeam profile is violated, one has to rely on a rather low-precision coefficient estimate obtained by roughly esti-mating the apparative parameters and resorting on theabove “Gaussian-beam” estimate.An exemplary kernel reconstruction procedure isshown in FIG. 1, where the OA signal p D at p sys =(1 cm / s , . , − . D ≈ . 75, is first ob-tained by solving the direct OA problem for Eq. (3)for an absorbing layer with µ a = 24 cm − in the range z = 0 − . p ) and blue ( p D ) curves inFIG. 1(a). The set ( p , p D ) is then used as inversioninput to compute the effective Volterra kernel for vari-ous sets of reconstruction parameters p rec = ( N, R ). Inparticular, considering N = 51, the minimal value of s ( a ⋆ , R ⋆ ) ≈ . 47 is attained at R ⋆ = 0 . 06 cm, see the insetof FIG. 1(b). As evident from the main plot of FIG. 1(b),the effective Volterra kernel for p rec = (51 , R ⋆ ) followsthe exact stress wave propagator for almost two ordersof magnitude up to c ∆ τ ≈ . 05 cm. Beyond that limit,the noticeable deviation between both does not seem toaffect the overall SSR s ( a , R ) too much. In this regard,note that the kernel approximated for the (non optimal)choice p rec = (51 , . 04 cm) exhibits a worse SSR. The inverse OA source reconstruction problem. Notethat the above Fourier-expansion approximation mightbe interpreted as a gauge procedure to adjust an ef-fective Volterra kernel K ( x ; a ⋆ , R ) for an (possibly un-known) apparative setup p sys , here indirectly accessiblethrough the diffraction transformation of the OA signal p D relative to p . That is, once the kernel reconstruc-tion (I.2) is accomplished for a set of reference curves( p , p D ) ref under p sys , the source reconstruction problem(I.1) might subsequently be tackled also for all other OAsignals measured under p sys by solving the OA Volterraintegral equation Eq. (3) in terms of a Picard-Lindel¨of“correction” scheme [16]. The latter is based on the con-tinued refinement of a putative solution, starting off from 0 0.05 0.1 0.15 c τ (cm) p ( a r b . un it s ) z D = -0.50 cm(a) p ( τ ) p D ( τ ) p PL ( τ ) -30-20-10 0 10 20 30 40 0 0.04 0.08 0.12 c ∆τ K eff ( ∆τ ) N =41 R = 0.1 cm 0 0.05 0.1 0.15 c τ (cm) p ( a r b . un it s ) (b) p PL ( τ ) p E ( τ ) p ( τ ) -50 0 50 100 150 0 0.04 0.08 0.12 c ∆τ K eff ( ∆τ ) N =51 R = 0.1 cm FIG. 2. (Color online) Inversion of OA signals to initial stressprofiles beyond the paraxial approximation. Both figures il-lustrate the kernel and source reconstruction procedures for(a) inversion of an OA signal featuring a top-hat irradia-tion source profile (see text). The main plot shows the in-put ( p , p D ) to the inversion procedure (solid black and bluelines, respectively) as well as the reconstructed initial stressprofile p PL (dashed red line), and, (b) inversion of an OA sig-nal resulting from an actual measurement [15]. The main plotshows the synthetic initial stress profile p (solid black line)used during the gauge procedure as well as the inversion in-put p E (orange line) for which the reconstructed initial stressprofile p PL (dashed red line) is obtained. In both figures, theinset illustrates the effective Volterra kernel resulting from theFourier-approximation. a properly guessed “predictor” p (0)PL ( τ ), improved succes-sively by solving p ( n +1)PL ( τ ) = p D ( τ )+ Z τ −∞ K ( τ − τ ′ ; a ⋆ , R ) p ( n )PL ( τ ′ ) d τ ′ . (9)From a practical point of view we terminated the iter-ative correction scheme as soon as the max-norm c n ≡k p ( n +1)PL ( τ ) − p ( n )PL ( τ ) k of two successive solutions decreasesbelow c n ≤ − . We here refer to the final estimate sim-ply as p PL . Note that, attempting a solution of (I.1) inthe acoustic near-field, a high-precision predictor can beobtained by using the initial guess p (0)PL ≡ p D . This isa reasonable choice since one might expect the changeof the OA near-field signal due to diffraction to be stillquite small. Further, source reconstruction in the acous-tic far-field might be started using a high-precision pre-dictor obtained by integrating the OA signal p D in thefar-field approximation [14]. In contrast to this, low-precision predictors for both cases can be obtained bysetting p (0)PL ≡ c , where, e.g., c = 0.The solution of the source reconstruction problemfor the OA signal p D used in the approximation ofthe Volterra kernel for the above setting p sys =(1 cm / s , . , − . p PL for p rec =(51 , R ⋆ ) and p does not come as a surprise since p D was used for the gauge procedure in the first place. Asa remedy we attempt a source reconstruction for a sec-ond independent OA signal, simulated for the same ap-parative setting only with two absorbing layers µ a, =24 cm − from z = 0 − . 05 cm and µ a, = 12 cm − from z = 0 . − . 12 cm. As evident from FIG. 1(c), inver-sion using the effective Volterra kernel from the previousgauge procedure yields a reconstructed stress profile p PL in excellent agreement with the underlying exact initialstress profile p . Inversion beyond the paraxial approximation. Giventhe apparent feasibility of the kernel reconstruction rou-tine as a gauge procedure to model the diffraction trans-formation of OA signals in terms of an effective stresswave propagator in the framework of the OA Volterraintegral equation, we next address the inversion of OAsignals to initial stress profiles beyond the paraxial ap-proximation. Therefore, we first consider a borderlinefar-field signal for a top-hat irradiation source f ( ~r ⊥ ) = ( , if | ~r ⊥ | ≤ ρ exp {− ( | ~r ⊥ | − ρ ) /a } , if | ~r ⊥ | > ρ , (10)recorded at the system parameters p sys =( c, ρ , a B , z D ) = (1 cm / s , . , . , − . 50 cm),and thus D = 2 | z D | / ( µ a ( a B + ρ )) ≈ . 04, obtainedvia an independent forward solver for the full OA waveequation designed for the solution of the OA Poissonintegral for layered media [5, 15]. The inversion resultsare summarized in FIG. 2(a), where the kernel recon-struction (inset) and source reconstruction (main plot)are shown for the parameter set p rec = (41 , . p and p PL suggests that the kernel reconstruction routinealso applies to a more general OA setting, based onthe full OA wave equation. Finally, we consider anOA signal resulting from an actual measurement onPVA hydrogel based tissue phantoms [15]. In thiscase we carefully estimated the apparative parameters p sys = (150000 cm / s , . 054 cm , . / s , − . µ a = 11 cm in the range z = 0 − . 095 cm,i.e. D ≈ . 73, in order to create a set of syntheticinput data by means of which an appropriate kernelgauge procedure can be carried out. The result ofthe procedure using p rec = (51 , . p E , we considered data withinthe interval cτ = [0 , . 15] cm, only. As evident from thefigure, the reconstructed stress profile p PL fits the signal p used in the gauge procedure remarkably well [17]. Conclusions. In the presented Letter we have intro-duced and discussed the kernel reconstruction problemin the paraxial approximation to the optoacoustic waveequation. We suggested a Fourier-expansion approachto approximate the Volterra kernel which takes a cen-tral role in the theoretical framework. The developedapproach proved useful as gauge procedure by meansof which the diffraction transformation experienced byOA signals can effectively be modeled, allowing to sub-sequently solve the source reconstruction problem in theunderlying apparative setting. From this numerical studywe found that the developed approach extends beyondthe framework of the paraxial approximation and alsoallows for the inversion of OA signals described by thefull OA wave equation. From a point of view of com-putational theoretical physics it would be tempting toexplore other kernel expansions in terms of generalizedFourier series as well as gauge procedures involving setsof measured pressure profiles only. Such investigationsare currently in progress with the aim to shed some morelight on this intriguing inverse problem in the field of op-toacoustics and to facilitate a complementary approachto conventional OA imaging. Acknowledgments. We thank A. Demircan for com-menting on an early draft of the manuscript and E. Blu-menr¨other for providing experimental data. This re-search work received funding from the VolkswagenS-tiftung within the “Nieders¨achsisches Vorab” programin the framework of the project “Hybrid Numerical Op-tics” (HYMNOS; Grant ZN 3061). Valuable discussionswithin the collaboration of projects MeDiOO and HYM-NOS at HOT are gratefully acknowledged. ∗ [email protected][1] G. J. Diebold, M. I. Khan, and S. M. Park, Science ,101 (1990); G. J. Diebold, T. Sun, and M. I. Khan, Phys.Rev. Lett. , 3384 (1991); I. G. Calasso, W. Craig, andG. J. Diebold, ibid . , 3550 (2001).[2] V. E. Gusev and A. A. Karabutov, Laser Optoacoustics (American Institute of Physics, 1993).[3] L. D. Landau and E. M. Lifshitz, Hydrodynamik (4th Ed.) (Akademie-Verlag (Berlin), 1981).[4] D. Colton and R. Kress, Inverse Acoustic and Electro-magnetic Scattering Theory (3rd Ed.) (Springer, 2013).[5] L. Wang, Photoacoustic Imaging and Spectroscopy , Opti-cal Science and Engineering (CRC Press, 2009).[6] P. Kuchment and L. Kunyansky, European Journal ofApplied Mathematics , 191 (2008).[7] M. Xu and L. V. Wang, Rev. Sci. Instr. , 041101(2006); L. V. Wang and S. Hu, Science , 1458 (2012); J.-M. Yang, C. Favazza, R. Chen, J. Yao, X. Cai,K. Maslov, Q. Zhou, K. K. Shung, and L. V. Wang, Na-ture medicine , 1297 (2012); L. Wang, J. Xia, J. Yao,K. I. Maslov, and L. V. Wang, Phys. Rev. Lett. ,204301 (2013); L. Wang, C. Zhang, and L. V. Wang, ibid . , 174301 (2014); I. Stoffels, S. Morscher, I. Hel-frich, U. Hillen, J. Leyh, N. C. Burton, T. C. P. Sardella,J. Claussen, T. D. Poeppel, H. S. Bachmann, A. Roesch,K. Griewank, D. Schadendorf, M. Gunzer, and J. Klode,Science Translational Medicine , 317ra199 (2015); ,319er8 (2015).[8] M. Agranovsky and P. Kuchment, Inverse Problems ,2089 (2007); X. L. De´an-Ben, A. Buehler, V. Ntziachris-tos, and D. Razansky, IEEE Transactions on MedicalImaging , 1922 (2012); Z. Belhachmi, T. Glatz, andO. Scherzer, Inverse Problems , 045005 (2016).[9] S. J. Norton and M. Linzer, IEEE Trans. Biomed. Eng., 202 (1981); M. Xu and L. V. Wang, Phys. Rev. E ,016706 (2005); P. Burgholzer, G. J. Matt, M. Haltmeier,and G. Paltauf, , 046706 (2007).[10] K. Chadan and P. Sabatier, Inverse Problems of Quan-tum Scattering Theory (Springer, 1989); B. Apagyi,G. Endr´edi, and P. Levay, Inverse and AlgebraicQuantum Scattering Theory , Lecture notes in physics (Springer, 1997); M. M¨unchow and W. Scheid, Phys.Rev. Lett. , 1299 (1980); O. Melchert, W. Scheid, andB. Apagyi, J. Phys. G: Nucl. Part. Phys. , 849 (2006).[11] A. C. Tam, Rev. Mod. Phys. , 381 (1986).[12] A. Karabutov, N. B. Podymova, and V. S. Letokhov,Appl. Phys. B , 545 (1996).[13] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in FORTRAN 77 (Cambridge Univer-sity Press, 1992).[14] J. Stritzel, O. Melchert, M. Wollweber, andB. Roth, “Direct and inverse solver for the 3D op-toacoustic Volterra equation,” (2016), (unpublished),arXiv:1606.04740.[15] E. Blumenr¨other, O. Melchert, M. Wollweber, andB. Roth, “Detection, numerical simulation and approxi-mate inversion of optoacoustic signals generated in multi-layered PVA hydrogel based tissue phantoms,” (2016),(unpublished), arXiv:1605.05657.[16] E. Hairer, S. P. Nørsett, and G. Wanner, Solving Or-dinary Differential Equations I (2nd rev. Ed.): NonstiffProblems (Springer, 1993).[17] A Python implementation of our code for the solu-tion of inverse problems (I.1) and (I.2) can be found at https://github.com/omelchert/INVERT.githttps://github.com/omelchert/INVERT.git