Photoacoustic tomography with direction dependent data: an exact series reconstruction approach

Photoacoustic image reconstruction often assumes that the restriction of the acoustic pressure on the detection surface is given. However, commonly used detectors often have a certain directivity and frequency dependence, in which case the measured data are more accurately described as a linear combination of the acoustic pressure and its normal derivative on the detection surface. In this paper, we consider the inverse source problem for data that are a combination of the acoustic pressure and its normal derivative. For the special case of a spherical detection geometry we derive exact frequency domain reconstruction formulas. We present numerical results showing the robustness and validity of the derived formulas. Moreover, we compare several different combinations of the pressure and its normal derivative showing that used measurement model significantly affects the recovered initial pressure.


Introduction
Photoacoustic tomography (PAT) is a hybrid imaging technique that combines high optical contrast with good ultrasonic resolution. It is based on the generation of an ultrasound wave inside an object of interest by pulses optical illumination. The initial pressure distribution of the induced sound wave encodes the optical absorption properties of the object, which are of great interest in medical diagnostics. Studies have shown that PAT allows high resolution imaging of small animals, and medical applications including early cancer diagnostics [2,29].
In a typical PAT setup, the induced acoustic waves are recorded by several point-like detectors outside the support of the investigated sample. Typically, the detectors are located on a surface S that fully or partially encloses the volume Ω in which the object of interest is contained (see figure 1). In most reconstruction approaches, the measured data are identified with the restriction of the acoustic pressure wave to the surface ∂Ω or sampled values of it, possible convolved in time with the detector impulse response function. Typical piezoelectric sensors, which are commonly used in PAT, however, have a directional dependence and are most sensitive in the direction orthogonal to the sensor surface. Moreover, at the resonance frequency they are more sensitive than at lower frequencies. As noted in [30], such measurement data are rather modeled by the combination of the pressure field and its normal derivative, than the pressure alone. This is also suggested by visual comparison with real data [26,27]. Systematic theoretical and experimental studies on PAT sensor modeling are interesting lines of future research.
Let f be the smooth compactly supported initial pressure distribuition, and p : R n × (0, ∞) → R the induced pressure wave that satisfies (1.1) In order to incorporate directional dependence, we model the data of a detector located at x ∈ ∂Ω by c 1 p(x, t) + c 2 n(x) · ∇p(x, t) for (x, t) ∈ ∂Ω × [0, T], (1.2) where n(x) · ∇p(x, t) is the normal derivative of the pressure, n(x) the outwards pointing normal of the surface at x, T the final measurement time, and c 1 , c 2 ∈ R are constants. The goal is to recover the initial pressure distribution f in (1.1) from data in (1.2). To the best of our knowledge, this paper is the first work investigating PAT with direction dependent data of the form (1.2) allowing arbitrary values of c 1 , c 2 . The standard image reconstruction problem in PAT corresponds to the special case c 1 = 0 and c 2 = 0 in our data model (1.2). Many methods for reconstructing the initial pressure distribution have been derived in the recent years in various situations. This includes, for example, different detection geometries with variable or constant speed of sound, or limited view situations. Theoretical questions concerning uniqueness and stability of the inverse source problem have also been investigated [3,14,16,18,28]. A practically important case assumes a constant speed of sound. In this case, several exact reconstruction formulas have been developed [1, 8, 9, 11-13, 15, 17, 19, 20, 23-25, 31]. Among these formulas so-called series expansion formulas provide very fast and accurate reconstructions. They give an expansion of the initial pressure in terms of eigenfunctions of the Laplacian [1,11,15,17,31].
The case where measurements are modeled by the normal derivative of the pressure (c 1 = 0 and c 2 = 0 in (1.2)) is studied much less. It has been considered in [6,7], where an explicit inversion formula of the backprojection type is derived for the case that the detection surface is a sphere in three spatial dimensions. Besides that, we are not aware of any results for the case c 2 = 0. We are not aware of any results when c 1 , c 2 are both non-vanishing.
In this paper, we provide series expansion reconstruction formulas for the direction dependent data model (1.2) allowing arbitrary values of c 1 , c 2 ∈ R, for the case that the measurement surface is a sphere. Our approach is based on expansions in spherical harmonics and an explicit formula relating the spherical harmonics coefficients of the direction dependent data, as a function of time, and the spherical harmonics coefficients of the initial pressure distribution f , as a function of distance to the origin (see section 2). By using Fourier Bessel series in the radial variable, in section 3 we derive inversion formulas that allows for a stable implementation. Numerical implementation and results are presented in section 4. The paper ends with a conclusion and outlook presented in section 5.

Preliminaries
In PAT, sound propagation is commonly described by an acoustic pressure wave p : R n × [0, ∞) → R that satisfies the initial value problem (1.1). In this work, we assume that the initial pressure distribution satisfies f ∈ C ∞ 0 (Ω). In PAT with direction dependent data we model measurement data by (1.2), where p is the solution of (1.1). The aim is to recover the initial pressure distribution f from such data. To the best of our knowledge, this is a new inverse problem for the wave equation that has not been considered so far. For the practical application, the cases n = 2 and n = 3 are of relevance. In particular, the case n = 2 appears from direction dependent measurements with integrating line detectors [4,15].
In particular, we study the case where Ω = B R (0) is the ball of radius R centered at the origin. In this situation, we can write the measurement data in the form (2.1) We write M c1,c2 for the operator that takes the initial data f in (1.1) to the corresponding boundary data. The constants c 1 , c 2 ∈ R are weights describing the contribution of the pressure and its normal derivative, respectively, to the measures data. The operator M 1,0 is the standard PAT forward operator, and the operator M 0,1 corresponds to the case where the data only consists of the normal derivative. Clearly, we have M c1,c2 = c 1 M 1,0 + c 2 M 0,1 .
For the derivation of the inversion formulas we use T = ∞. Because of the strict Huygens principle, in the case of odd spatial dimension, we have p(x, t) = 0 for t > 2R, which yields exact inversion for any measurement time T 2R. For even spatial dimension, this is not the case and, strictly taken, our inversion formulas are exact only for T = ∞. Deriving exact formulas using data over a finite time interval [0, T] in even dimension is an interesting open issue. For the standard PAT operator M 1,0 such a formula has been derived in [7, theorem 1.4].

Notation
Our approach is based on Fourier methods that lead to a relation between measurement data and the initial pressure distribution in the frequency domain. We denote the involved transforms bŷ which are the Fourier, cosine and sine transforms, respectively. Moreover, we denote by Y ,k (θ) the spherical harmonics [22], which form a complete orthonormal system in L 2 (S n−1 ). In particular, where N(n, ) = (2 + n − 2)(n + − 3)!/( !(n − 2)!) for ∈ N and N(n, 0) := 1. We write J ν for the ν th order Bessel function of the first kind and denote by w j, the j th positive root of J +(n−2)/2 .

Relations between spherical harmonics coefficients
For the following we use the spherical harmonics expansions of the measurement data g = M c1,c2 f in (2.1) and the Fourier transform f of the initial pressure, The following lemma is the key to our results.
be the initial pressure distribution in (1.1) and let g = M c1,c2 f be the corresponding measurement data given by (2.1). Then, (2.4) Proof. Let p be the solution of (1.1).
Changing to spherical coordinates, λω ← ξ, and using the spherical harmonics expansion (2.2) we can write the direction dependent measurements as where we used the abbreviations According to the Fuck-Hecke theorem [22] A ,k (θ) = (2π) The last two displayed equations imply (2.6) Comparing expansions (2.2) and (2.6) we conclude that This becomes (2.4) after applying the inversion formula g = 2 π CCg for the cosine transform; see [10, equation (7.28)]. □ As a first application of lemma 1, we obtain a range condition for the classical PAT forward operator M 1,0 that maps the initial data f to the solution of (1.1) on the boundary. More precisely, we have the following result.

Corollary 1 (Range condition). Let
As another application of lemma 1, we derive the following inversion formula.
Together with (2.3), this gives the desired result. □ Corollary 2 gives an exact inversion formula for reconstructing any smooth compactly supported function f ∈ C ∞ 0 (R n ) from data M c1,c2 f . In order to actually implement the formula one requires a proper discretization of the integral. In real situations, only noisy data g δ ∈ L 2 (S R × (0, ∞)) are available and the cosine transform C{g δ ,k }(λ) will not vanish at the roots of the denominator (c 1 + c 2 R −1 )J + n−2 2 (Rλ) − c 2 λJ +n/2 (Rλ) that appears in formula (2.7). This means that corollary 2 behaves unstable close to the roots and cannot be directly used to reconstruct the initial pressure density f . In the following section, we present a strategy to overcome this issue.

Stable series inversion formulas
In practice, reconstructing the initial pressure by the inversion formula stated in corollary 2 is unstable. For its stabilization, in this section we derive a series inversion formula using a Fourier Bessel series, similar to [15,31]. We derive different inversion formulas for the cases c 2 = 0 and c 2 = 0, respectively. For the following recall that w j, denotes the j th positive zero of J +(n−2)/2 .
The following theorem is the main result of this paper.
Proof. Because f ,k is of class C ∞ and is compactly supported, we can expand f ,k (ρ)ρ n−2 2 into a Fourier Bessel series [10], where for the second line, we used (2.8). From lemma 1, we have We can therefore evaluate (3.4) at w = w j, and solve for f ,k w R . Together with (3.3), we obtain .

(3.5)
For the case c 2 = 0, equation (3.4) giveŝ For the zeros w = w j, of J + n−2 2 (w) this is an indeterminate form, which can be evaluated with L'Hospital's rule. We have Using this and the identity ∂ w C{g ,k }(wR −1 ) = −R −1 S{tg ,k (t)}(wR −1 ), and applying L'Hospital's rule give . Therefore equation (3.3) yields By combing (3.5) and (3.6) with the spherical harmonics expansion (2.3), we obtain the desired inversion formulas (3.1) and (3.2). □ Note that the inversion formula (3.1) for M c1,c2 might be slightly surprising because it holds for any c 2 = 0 regardless of the value of c 1 . However, this independence is an immediate consequence of the decomposition M c1,c2 f = c 1 M 1,0 + c 2 M 0,1 , the range condition in corollary 1 for the operator M 1,0 , and the particular form of the right hand side in (3.1). Equation

Numerical experiments
In this section, we present reconstruction results with the inversion formulas in theorem 1 using data g = M c1,c2 f for different combinations of c 1 , c 2 . We consider the case of n = 2 spatial dimensions and use R = 1. The 2D case arises in applications where the acoustic pressure is measured by integrating line detectors [4,26].

Discretization and data simulation
In order to compute the pressure field and its gradient restricted to the unit circle, we use a discrete wave propagation model on a 2D equidistant grid. To numerically compute the solution at the discretization points we use the k-space method [5,21] implemented as described in [14]. Using the k-space method, we obtain the values of the acoustic pressure field and compute its gradient by symmetric differences. The pressure and its normal derivative on the circle are obtained by linear interpolation. For the presented numerical results, the initial pressure (figure 2) is given on 280 × 280 Cartesian grid points in [−1,1] 2 and the data M c1,c2 f are simulated for N θ = 300 equidistant sensor locations on unit circle. At each sensor, we use N t = 1600 equidistant time samples in the measurement interval [0,6].
The top left image in figure 2 shows the initial pressure distribution f used our simulations. The top right image shows the classical PAT data M 1,0 f + ξ 1,0 , and the bottom right image shows the normal derivative M 0,1 f + ξ 0,1 , which has been re-scaled such that M 1,0 f and M 0,1 f have the same 2 -norm. The bottom left image shows the mixed data M 1,1 f + ξ 1,1 . In all cases, we have added Gaussian white noise ξ c1,c2 with a standard deviation of 50% of the 2 -norm M c1,c2 f 2 := ( m,n |M c1,c2 f [m, n]| 2 /(N θ N t )) 1/2 , resulting in a relative 2 -data error ξ c1,c2 2 / M c1,c2 f 2 of 50% in all cases. Note that we simulated the data until T = 6 but only show data until t = 2.2 in figure 2, because after T = 2 the data M c1,c2 f are smooth and monotonically decreasing.

Implementation of the inversion formulas
In order to reconstruct the initial pressure distribution, we implement discrete versions of the reconstruction formulas (3.1) and (3.2), which for n = 2 and R = 1 are given by   e ikϕ (for c 2 = 1),   e ikϕ (for c 1 = 1).
Below we briefly describe the discretization of inversion formula (4.1). The numerical reconstruction uses discrete samples g where θ m = 2π (m − 1) /N θ and t n = T (n − 1) /N t . First, the Fourier coefficients g k in the angular variable are computed by the FFT algorithm. Next, the cosine transform is approximated by C{g k }[ j] := T Nt Nt n=1 g k [n] cos(ω j,k t n ). This is implemented by matrix-vector multiplication where the entries cos(ω j,k t n ) are pre-computed and stored. Finally, we evaluate (4.1) by truncating both sums. Formula (4.2) is discretized in an analogous manner. This results in the discrete versions of (4.1) and (4.2)   e ikϕ (for c 1 = 1).
(4.4) In both formulas, the inner sums are evaluated by matrix-vector multiplications where the required matrix elements are pre-computed and stored. The outer sums are evaluated with the inverse FFT algorithms. Application of (4.3) (or (4.4)) with a standard Matlab implementation on a desktop PC with 16 GB RAM and 3.40 MHz eight-core processor with N θ = 300, N r = 180 and N t = 1200 takes about 1.2 s. Note that ω j,k are the positive roots of J k and therefore J k+1 (ω j,k ) = 0 which implies that neither (4.3) nor (4.4) suffer from a division by zero problem.  figure 4, which again support exactness of inversion formula and further shows their stability with respect to Gaussian noise. Figure 5 shows a quantitative error analysis in which the relative reconstruction error is shown depending on the noise level of the data.

Discussion
First, note that the theoretically exact inversion formulas (3.1) and (3.2) use data for all times whereas the discrete counterparts only use data up to finite time. Nevertheless, the discrete formula (4.3) applied to M 1,0 f and (4.4) applied to M 0,1 f give visually satisfactory results. This is consistent with our previous observations [4,8,15]. However, the relative reconstruction error is quite large, even for data without noise. We mainly address this due to the required truncation of the data. In future work we will therefore address this issue and aim deriving exact inversion formula which only use data until a finite time T.
Second, up to discretization and truncation errors, due to the range condition in corollary 1, the inversion formula (4.3) applied to classical PAT data M 1,0 f should yield the zero image. From figures 3 (top right) and 4 (top right) we see that this is clearly not the case. Again, this  2) (exact for c 2 = 0) and the right column to reconstructions with formula (3.1) (exact for c 2 = 0). Top: standard PAT data M 1,0 f + ξ 1,0 . Middle: mixed data M 1,1 f + ξ 1,1 . Bottom: normal derivative M 0,1 f + ξ 0,1 . In all cases, the added Gaussian white noise ξ c1,c2 has a standard deviation of 50% of the 2 -norm M c1,c2 f 2 . Only for the top left, middle right and bottom right images the applied inversion formula matches the data.
is mainly due to the data truncation, which we have verified (results not shown) by varying T. The introduced violation of the range condition is also the reason that the inversion formula applied to M 1,1 f yields worse results than the result for M 0,1 f . Finally, applying the inversion formula (3.2) to data M 1,1 f results in amplified boundary structures. Visually, the amplification is quite appealing which may be the reason that commonly data are modeled by M 1,0 f instead of the more general model M c1,c2 f . Quantitative analyzing the effect of M 0,1 f component to (3.2) is an interesting open issue.

Conclusion
We investigated PAT with the direction dependent data model (2.1), which uses linear combinations M c1,c2 f = c 1 M 1,0 f + c 2 M 0,1 f of the acoustic pressure and its normal derivative. We developed an exact and stable reconstruction formula for the special case of spherical detection geometry. Numerical results show the validity of the proposed approach. Investigating such data for more general detection geometries in PAT is an interesting line of future research. Moreover, deriving an explicit inversion formula that only used data for t T < ∞ is an important future aspect. Finally, in future work we will test our formulas on experimental data and investigate which values of c 1 , c 2 actually are the best to accurately model experimental PAT data, for example using small piezoelectric sensors or piezoelectric line sensors. Figure 5. Quantitative error analysis: the plots show the relative reconstruction error f − f rec 2 / f 2 (vertical axis) versus the relative data error g − g δ 2 / g 2 (horizontal axis). Note that the axis are not in percent, such that value 1 on the axis corresponds to 100%. Left: results for formula (3.2). Right: results for formula (3.1).