Quantitative reconstructions in multi-modal photoacoustic and optical coherence tomography imaging

In this paper we perform quantitative reconstruction of the electric susceptibility and the Gr\"uneisen parameter of a non-magnetic linear dielectric medium using measurement of a multi-modal photoacoustic and optical coherence tomography system. We consider the mathematical model presented in [11], where a Fredholm integral equation of the first kind for the Gr\"uneisen parameter was derived. For the numerical solution of the integral equation we consider a Galerkin type method.


Introduction
Tomographic imaging techniques visualize the inner structure of probes. Particularly relevant for this work are Optical Coherence Tomography (OCT) and Photoacoustic (PAT). In OCT a sample is placed in an interferometer and is illuminated by light pulses. Then, the backscattered light is measured far from the medium, see for instance [7,9,14]. PAT visualizes the capability of a medium to transform optical (infrared) waves into ultrasound waves to be measured on the surface of the medium [17,28,30]. PAT is called coupled physics imaging technique since it combines two kind of waves [1]. As stand alone imaging techniques PAT and OCT are not capable of recovering all diagnostically relevant physical parameters, but only some combinations of them, see [3] for PAT and [12] for OCT.
Recently setups which combine different imaging modalities, have been investigated mathematically with the objective to reconstruct more diagnostically relevant physical parameters from the measurements. Particular applications are coupled physics imaging systems and elastography [2,20,29], to name but a few. We refer to these techniques as hybrid imaging or multi-modal imaging systems. Note that in the mathematical literature the name hybrid imaging is also used for coupled physics imaging.
In this work we consider the multi-modal PAT/OCT system, developed for imaging biological tissues, see [10,21,22,23,31]. We show that with such a system, in contrast to the single modality setups, we obtain sufficient measurements which allow us to extract quantitative information on the electric susceptibility and the Grüneisen parameter of the sample. In the multi-modal PAT/OCT system, two different excitation laser systems, both operating in the same wavelength range, are used. The PAT and OCT scans are performed sequentially and vary a lot in acquisition times (around 5 minutes in PAT and less than 30 seconds in OCT). The obtained PAT and OCT images are co-registered afterwards.
In Section 2, we describe mathematically the multi-modal PAT/OCT setup. We use the model, from [11], based on Maxwell's equations for the electric permittivity. In Section 3, we present the equivalence of the inverse problem of recovering both optical parameters with the solution of a Fredholm integral equation of the first kind for the Grüneisen parameter. Here the kernel of the integral operator depends on the PAT measurements.
We propose a numerical reconstruction method based on a Galerkin method using a series expansion of the unknown functions with respect to Hermite functions, see Section 4. The discretization of the continuous integral operator results in a system of linear algebraic equations. We solve the matrix equation using Tikhonov regularization. Numerical results which justify the feasibility of the proposed method are presented in Section 5.

The multi-modal PAT/OCT system
We consider the two modalities independently. Full field illumination is used in PAT and focused in OCT. The medium in OCT is illuminated by a Gaussian light. However, we can assume that the plane wave illumination is still valid [14].

Light propagation.
We consider macroscopic Maxwell's equations in order to model the interaction of the incoming light with the sample. These equations describe the time evolution of the electric and magnetic fields E and B for given charge density ρ and electric current J: where D ≡ E + 4πP is the electric displacement and H ≡ B − 4πM denotes the effective magnetic field, related to the electric and magnetic polarization fields P and M, respectively. We specify the material properties of the medium. Definition 2.1.
The medium is called non-magnetic if M = 0, and perfect linear dielectric and isotropic if there exist a scalar function χ ∈ C ∞ c (R × R 3 ; R) the electric susceptibility, with χ(t, x) = 0 for all t < 0, x ∈ R 3 , (this property is referenced as causality), such that The electric susceptibility describes the optical properties of the medium and is the parameter to be determined. In addition, we assume that the medium has no free charges, meaning ρ = 0 in (1a).

PAT measurements.
Let the medium be defined as in Definition 2.1. Then, we estimate the averaged change in energy density around a point x, for every ν, by In order to derive the above formula we have to consider the interaction of the medium with the incoming electromagnetic wave locally. For a derivation, using microscopic Maxwell's equations, see for instance [11,Section 4].
The laser pulse is absorbed by the medium and part of it is transformed into heat. This generates a pressure wave which is then measured on the object surface. Since the laser pulse is typically very short, the propagation of the acoustic wave during thermal absorption can be neglected. Then, we consider as PAT measurements the initial pressure density p which is proportional to the absorbed energy The proportionality factor Γ is the Grüneisen parameter, a parameter which, together with the susceptibility χ, describes the optical properties of our medium.

OCT measurements.
In the frequency domain, the equation (3) and the condition (5) result to an integral equation of Lippmann-Schwinger type [8,12].
Let the medium be defined as in Definition 2.1 and E (0) ν as in Definition 2.2. If E ν is a solution of (3) with initial values (5), then its Fourier transform solves the Lippmann-Schwinger integral equationÊ Due to the limiting penetration depth of OCT (1 to 2 millimeters), the medium can be considered as weakly scattering, since only single scattering events will be measured. In addition, in OCT the measurements are performed in a distance much larger compared to the size of the medium.
The Born approximation allows us to obtain an explicit form forÊ ν from the Lippmann-Schwinger equation (10). In the limiting caseχ → 0, we take the first order approximation of the electric field by replacingÊ ν withÊ (0) ν in the integrand of (10).
We write x in spherical coordinates x = ρϑ, ρ > 0, ϑ ∈ S 2 . Under the far-field approximation, we consider the asymptotic behavior of the expression (10) for ρ → ∞, uniformly in θ Then, we definê as the electric field considering both approximations.
The approximated backscattered lightÊ is combined with a known back-reflected field and its correlation is measured at each point on the detector surface. Under some assumptions on the incident illumination we state that what we actually measure in OCT is the backscattered light at a detector placed far from the medium [12,Proposition 8].
Then, we formulate the direct problem as: Given a medium as in Definition 2.1 with susceptibility χ and Grüneisen parameter Γ, and incident illumination E (0) ν of the form (6), the direct problem is to find the PAT measurements p ν (x), x ∈ Ω, ν > 0, given by (9), and the OCT measurements

The Inverse Problem
In the following the assumptions on the medium (Definition 2.1) hold and especially the causality of χ. We denote byχ the three-dimensional Fourier transform ofχ with respect to spacẽ (11) and simple calculations, see [12,Proposition 9], provide us with the dataχ However, in practice, these data are incomplete because of the band-limited source and size of the detector. Thus, we get the spatial and temporal Fourier transform of χ only in a subset of Then, the inverse problem we address here reads: Given a medium as in Definition 2.1 and incident fields E (0) ν of the form (6) for all ν > 0, the inverse problem is to recover the parametersχ and Γ given the internal PAT measurements p ν (x), for x ∈ Ω, and all ν > 0, given by (9), and the external OCT dataχ(ω, ω c (ϑ + e 3 )), ω ∈ R \ {0}, ϑ ∈ S 2 + , given by (12).
Similar inverse problems have been considered in [4,6] where the far-field measurements from OCT are replaced by boundary measurements and in [5] for the diffusion approximation of the radiative transfer equation.
To present an equivalent formulation of the inverse problem, we assume that in both imaging techniques, we illuminate with multiple laser pulses with small spectrum centered around different frequencies. This setup describes swept-source OCT and multi-frequency PAT measurements.
First we describe the PAT measurements for multiple laser pulses. We combine (8) and (9) to get where E ν is the electric field generated by the laser pulse E (0) ν . Using the Fourier transform of (2a) we derive

Remark 3.2:
In the case of nonlinear medium, the polarization field P is usually expressed as a power series of the electric field E. Then, the third order term contributes to the so-called two-photon absorbed energy [15]. We refer to [27] for reconstructions in two-photon PAT.
As in OCT, we replace E ν by the initial pulse E (0) ν and we approximate the PAT data by The support off ν is localized around the frequency ν, see (7). Thus, we get in the limit ε → 0 (for constant norm f ν 2 ) that . Then, we get asymptotically We assume measurements for all frequencies ν > 0. Recall that χ is a causal real valued function. Then the real part ofχ can be completely determined from the imaginary part via the Kramers- Here, H denotes the Hilbert transform with respect to frequency We use the above two equations in order to describe the OCT data (12). Then we end up with the Fredholm integral equation for the Grüneisen parameter Γ. Once (14) is solved, we can easily recover the imaginary part ofχ from equation (13).

Remark 3.3:
If the medium is a perturbation of a single material then the above equation is transformed to a Fredholm integral equation of the second kind for a new function depending on 1 Γ [11]. At least in this simplified setting, we find that using the multi-modal model PAT/OCT we can (uniquely) determine the Grüneisen parameter and the susceptibility χ describing the absorption and scattering properties of the medium.
In the following section we present a Galerkin type method for the numerical solution of equation (14) considering two types of media. There exist also other projection methods for the numerical solution of integral equations, the collocation method and the method of moments and quadrature methods, like the Nyström method [16,19,25].

Numerical Implementation
Without loss of generality we set c = 1 and we specify Ω = [−l, l] 3 . For the numerical examples we have to introduce the parameterΓ, related to the physical parameter Γ, which satisfies , for x ∈ Ω, and suppΓ ⊂ Ω L , where Ω L = [−L, L] 3 , for L > l. This is possible since Γ ≥ Γ 0 > 0, with Γ 0 ∼ 1 in biological tissues. In addition, we do not consider the restrictions on the frequency in the following analysis.

Medium with depth-dependent coefficients.
In the first example, we assume that both parameters Γ andχ are only depth-dependent, meaning that they vary only in the incident direction. Then,Γ andχ admit the forms respectively. Here 1 denotes the characteristic function. This case represents media which have a multilayer structure with depth-dependent properties, like the human skin. If the illumination is focused to a small region inside the object and this region is small enough such that the functions can be assumed constant in both directions e 1 and e 2 , we get the above forms.
Thus the problem reduces to the problem of recovering a one-dimensional function γ. We do not consider the two-dimensional detector array but only the measurements at the single point detector located at the position (0, 0, d), meaning we set ϑ = e 3 . Then, the equation (14) takes the simplified form where m(ω) := ( π l ) 2χ (ω, 2ωe 3 ). Let p ∈ (L 2 (R)) 2 and γ ∈ L 2 (R). Since the kernel of the integral operator and the right-hand side in (15) have specific structures containing Hilbert and Fourier transforms we consider as orthonormal basis of L 2 (R) the Hermite functions h k , k ∈ N 0 . In addition the multi-dimensional Hermite functions can be written as sum of products of the usual Hermite functions. Their properties are given in Section 6. Other choices are also possible, especially when we treat the three-dimensional problem with real data, for instant using wavelets as basis functions.
Let x ∈ R. The Hermite polynomials are defined by the formula The normalized Hermite functions are given by where α k = (2 k k! √ π) − 1 2 . The functions h k satisfy the orthonormality condition Proposition 4.1.
Let x ∈ R. We consider the expansion with coefficients γ k ∈ R, k ∈ N 0 and with coefficients p k,l ∈ R, k, l ∈ N 0 . Then, if γ satisfies the integral equation (15), the coefficients γ k , k ∈ N 0 solve the equation where , and q := j + l − 2n − 2r.
Before proving this Proposition, we state the following lemma. Its proof is presented in Section 6.
Proof (Proposition 4.1): Using the expansion (18) and considering (44), we get The coefficientsp k,l , using (45), are given bỹ We substitute the above expansions and (17) in (15) and we obtain We rewrite the product of the two Hermite functions in the integrand using the formula (41) and we change variables to get Then, equation (21) using (20) takes the form wherem(ω) = 2e ω 2 4 m(ω).
Using (16) and (43), we get We substitute this expansion in (22) where for simplicity we set q := j + l − 2n − 2r. Again the last product using (41) admits the form We expand also the data using the same basis functions and in order to obtain a linear equation for γ j we have to enlarge the index of the last sum. We set s := k + q − 2u and for k+q−s 2 ∈ N 0 we reformulate (24) using the above formulas as Equating the coefficients in the above equation yields (19).
The final step, for the Galerkin method, is to consider a finite dimensional subset of L 2 (R), meaning restrict ourselves to a finite number of coefficients. Let j, k, l = 0, ..., N − 1. Then the definitions used in the above analysis gives q = 0, ..., 2(N − 1) and s = 0, ..., 3(N − 1). Finally, the discrete linear system of (19) reads where

Medium with coefficients constant in one direction.
In this case, we assume that both parameters are constant only in one direction, let us say in e 1 . Then,Γ andχ take the forms respectively. This assumption results to a two-dimensional function γ, a case more involved compared to Section 4.1 that approximates better the unconditional general problem. Here, we need two-dimensional data, thus we have to consider measurements for all frequencies in a one-dimensional array, modeling measurement points on a line.
The equation (14), for c = 1, now takes the form
We use (46), for k = (k, l), and we expand γ as where the coefficients γ k,l are defined by and we assume the expansion where Then, if γ solves the integral equation (26), its coefficients γ k,l , k, l ∈ N 0 satisfy the equation for (p a,n,u + ip a,n,u ) min(k,n) r=0 β k,n,r ζ k+n−2r ϑ k+n−2r Proof: The equation (26) using the expansions (27) and (28)  We apply twice the formula (41) for the product of two Hermite functions, to obtain ∞ k,l=0 γ k,l ∞ a,n,u=0 (p a,n,u + ip a,n,u )h a (ω) The last two integrals can be again simplified using lemma 4.2. We get We expand again the product h a (ω)h s−2j (ω) using (41) as We set µ := a + s − 2j − 2t and for a+s−2j−µ 2 ∈ N 0 the above sum can be rewritten as We expand the right-hand side of (30) using the same basis functions Then, the equation (30) using the above formulas and equating the coefficients results in equation (29).

Numerical Results
Both linear systems derived in the previous section admit the general form In the case of depth-dependent coefficients, see Section 4.1, we have and in the case of constant in one direction coefficients, see Section 4.2, we get We approximate the solution of (33) by minimizing the Tikhonov functional where I is the identity matrix with dimensions depending on each problem. We consider also noisy data for both measurements data, the pressure p and the OCT data m, with respect to the L 2 norm for given noise levels δ p , δ m and v = v 1 + iv 2 , w = w 1 + iw 2 , for v 1 , v 2 , w 1 and w 2 normally identically distributed, independent random variables.
We present reconstructions for different functions γ (related to 1/Γ) and ψ (related toχ) for both cases of media. As OCT data we consider the functionχ (using the Fourier transform and the Kramers-Kronig relation) and to construct the simulated PAT data we have to assume that both functions have similar behavior such that ratioχ/γ (see (13)) is still integrable. In all figures we plot the spatial domain Ω L . and We set Ω = [−3.5, 3.5] and Ω L = [ −4, 4] such that supp ψ(ω, ·) ⊂ Ω, and supp γ ⊂ Ω L , and we restrict ourselves to ω ∈ W := [ −4, 4]. We consider data with δ p = δ m = 3% noise. The results are presented in Figure 1 for regularization parameter λ = 10 −4 and different values of N. Here, we see the improvement in the reconstructions as N increases.
In the second example, we use the same function ψ and we consider as γ the function We keep all the parameters the same as in the first example. In Figure 2, we see the reconstructions for different number of coefficients. In the third example, we consider such that again supp ψ(ω, ·) ⊂ [−3.5, 3.5], see the left picture in Figure 3. We present the reconstructions of m(ψ) using the form (34) for γ, while keeping all the other parameters the same. We set N = 15 coefficients. We present the results for m(ψ(ω, 1)), ω ∈ W, (center picture) and m(ψ (3, x)), x ∈ Ω, (right picture) in Figure 3.

Examples with coefficients constant in one direction (see Section 4.2).
Here, the measurements are given at points on a line. We consider the minimum amount of measurement points in order to have an exactly determined system (32) in our examples. In the following examples we keep the same noise levels δ p = δ m = 3% and we obtain the regularization parameter using the L-curve criterion [18].
Let x, y ∈ R. In the fourth example, we consider  Figure 4. The results for the cross-section of the imaginary part of ψ, given by equation (38), at frequency ω = 0 are presented in Figure 5.
In the last example the unknown function is given by The size of the medium is kept the same as in the previous example and we set W = [−2, 2]. Here, we want to test the performance of our numerical scheme with respect to the number of the detection directions ϑ (k) . In the case of three measurement directions, see (32), we consider ϑ (1) = 0, cos( 5π 12 ), sin( 5π 12 ) , ϑ (2) = (0, 0, 1) , ϑ (3) = 0, − cos( 5π 12 ), sin( 5π 12 ) .  The reconstructions for N = 5 coefficients are presented in Figure 6, where we set to zero the negative values. We set the imaginary part of ψ to be The reconstruction for W = [−3, 3] are given in Figure 7, where we see the improvement of the results with respect to the Fourier coefficients. In the first case we set N = 5 and we consider one detection direction. In the second case, we use N = 10 coefficients and two measurement points in the directions: ϑ (1) = 0, cos( 7π 16 ), sin( 7π 16 ) , ϑ (2) = 0, − cos( 7π 16 ), sin( 7π 16 ) .

Conclusions
In this work we considered the inverse problem to reconstruct quantitatively the electric susceptibility and the Grüneisen parameter of a non-magnetic linear dielectric medium from measurements with the multi-modal tomographic system of Photoacoustic and Optical Coherence Tomography. Our scheme is based on the numerical solution of a Fredholm integral equation of the first kind for the Grüneisen parameter using a Galerkin type method. We presented numerical results for different kinds of media.
where the coefficients f k are defined by The Hilbert transform of f admits the expansion wheref k are given by [26] For k ∈ N d 0 and x ∈ R d , we define the kth Hermite polynomial as Now we present the proof of Lemma 4.2.
Proof (Lemma 4.2): We consider the convolution theorem for the inverse Fourier transform and the above properties. To compute the last integral we apply the formula (42) The last two equations result to (20).