Quantitative Photoacoustic Imaging in the Acoustic Regime using SPIM

TL;DR

This paper introduces a method for quantitative photoacoustic imaging that simultaneously reconstructs absorption, sound speed, and density in a medium using a generalized wave equation and SPIM techniques.

Contribution

It presents analytical formulas for reconstructing multiple parameters in photoacoustic imaging based on a linearized model with SPIM, extending standard methods.

Findings

Reconstruction formulas for absorption, sound speed, and density.

Simultaneous imaging of multiple medium parameters.

Extension of photoacoustic imaging capabilities.

Abstract

While in standard photoacoustic imaging the propagation of sound waves is modeled by the standard wave equation, our approach is based on a generalized wave equation with variable sound speed and material density, respectively. In this paper we present an approach for photoacoustic imaging, which in addition to recovering of the absorption density parameter, the imaging parameter of standard photoacoustics, also allows to reconstruct the spatially varying sound speed and density, respectively, of the medium. We provide analytical reconstruction formulas for all three parameters based in a linearized model based on single plane illumination microscopy (SPIM) techniques.

Equations

$$E_{r,\theta }=\{x\in \text{mathds}{R}^3\mid x\cdot \theta =r\}$$

$$P_{r,\theta }^{(0)}(x)=\gamma (x)\mu (x){\Phi }_{r,\theta }(x).$$

$${\sigma }_{r,\theta ;ij}(t,x)=K(x){\mathrm{div}}_x{u}_{r,\theta }(t,x){\delta }_{ij},$$

$$\rho (x){\partial }_{tt}{u}_{r,\theta ;i}(t,x)=\sum _{j=1}^3{\partial }_{x_j}{\sigma }_{r,\theta ;ij}(t,x)={\partial }_{x_i}(K(x){\mathrm{div}}_x{u}_{r,\theta }(t,x)).$$

$$P_{r,\theta }(t,x)=-K(x){\mathrm{div}}_x{u}_{r,\theta }(t,x).$$

$$\rho (x){\partial }_{tt}{\mathrm{div}}_x{u}_{r,\theta }(t,x)+\nabla \rho (x)\cdot {\partial }_{tt}{u}_{r,\theta }(t,x)={\Delta }_x(K(x){\mathrm{div}}_x{u}_{r,\theta }(t,x)).$$

$$\frac{\rho (x)}{K(x)}{\partial }_{tt}{P}_{r,\theta }(t,x)+\frac{\nabla \rho (x)}{\rho (x)}\cdot {\nabla }_x{P}_{r,\theta }(t,x)={\Delta }_x{P}_{r,\theta }(t,x).$$

$${\partial }_{tt}{P}_{r,\theta }(t,x)=\rho (x)\alpha (x){\mathrm{div}}_x\left(\frac{1}{\rho (x)}{\nabla }_x{P}_{r,\theta }(t,x)\right).$$

$$\begin{array}{rl}{\partial }_t{P}_{r,\theta }(0,x)& =0,\\ P_{r,\theta }(0,x)& =P_{r,\theta }^{(0)}(x).\end{array}$$

$$m_{r,\theta }(t,\xi )=P_{r,\theta }(t,\xi )$$

$$\rho (x;\epsilon )$$

$$\alpha (x;\epsilon )$$

$$P_{r,\theta }(t,x;\epsilon )=P_{r,\theta ,0}(t,x)+\epsilon P_{r,\theta ,1}(t,x)+{\epsilon }^2{P}_{r,\theta ,2}(t,x;\epsilon ),$$

$${\partial }_{tt}{P}_{r,\theta ,0}(t,x)$$

$${\partial }_t{P}_{r,\theta ,0}(0,x)$$

$$P_{r,\theta ,0}(0,x)$$

$${\partial }_{tt}{P}_{r,\theta ,1}(t,x)$$

$${\partial }_t{P}_{r,\theta ,1}(0,x)$$

$$P_{r,\theta ,1}(0,x)$$

$${\partial }_{tt}{P}_{r,\theta ,2}(t,x;\epsilon )$$

$${\partial }_t{P}_{r,\theta ,2}(0,x;\epsilon )$$

$$P_{r,\theta ,2}(0,x;\epsilon )$$

$$F(t,x;\epsilon )=\frac{1}{\rho (x;\epsilon )}\frac{{\rho }_1(x){\alpha }_0-{\rho }_0{\alpha }_1(x)}{{\rho }_0}\nabla {\rho }_1(x)\cdot \nabla P_{r,\theta ,0}(t,x)+{\alpha }_1(x)\Delta P_{r,\theta ,1}(t,x)-\frac{\alpha (x;\epsilon )}{\rho (x;\epsilon )}\nabla {\rho }_1(x)\cdot \nabla P_{r,\theta ,1}(t,x)$$

$$P_{r,\theta }(t,x;\epsilon )=\frac{{\partial }^4}{\partial t^4}\left({\int }_{B_t(x)}K(t,x,y;\epsilon )P_{r,\theta }^{(0)}(y)\mathrm{d}y\right)+{\epsilon }^2{P}_{r,\theta ,2}(t,x;\epsilon ),$$

$$\begin{array}{rl}K(t,x,y;\epsilon )& =\frac{(t-|y-x|{)}^2}{8\pi |y-x|}(1-\epsilon {\alpha }_1(y))+{\int }_{{\mathcal{E}}_t(x,y)}\frac{1}{16{\pi }^2|z-x||z-y|}\epsilon {\alpha }_1(z)\mathrm{d}z\\ & +{\int }_{{\mathcal{E}}_t(x,y)}\frac{t-|z-x|-|z-y|}{16{\pi }^2|z-x||z-y|}\epsilon \nabla {\rho }_1(z)\cdot \frac{z-y}{|z-y|}\mathrm{d}z\\ & +{\int }_{{\mathcal{E}}_t(x,y)}\frac{(t-|z-x|-|z-y|{)}^2}{32{\pi }^2|z-x||z-y{|}^2}\epsilon \nabla {\rho }_1(z)\cdot \frac{z-y}{|z-y|}\mathrm{d}z,\end{array}$$

$${\mathcal{E}}_t(x,y)=\{z\in \text{mathds}{R}^3\mid |z-x|+|z-y|\le t\}$$

$$\underset{0\le t\le T}{\sup }\left(\Vert P_{r,\theta ,2}(t,\cdot ;\epsilon ){\Vert }_{H^{k+1}({\mathbb{R}}^3)}+\Vert {\partial }_t{P}_{r,\theta ,2}(t,\cdot ;\epsilon ){\Vert }_{H^k({\mathbb{R}}^3)}\right)\le C_{k,T;\epsilon }{\int }_0^T\Vert F(t,\cdot ;\epsilon ){\Vert }_{H^k({\mathbb{R}}^3)}\mathrm{d}t$$

$$P_{r,\theta ,0}(t,x)=\frac{\partial }{\partial t}\left(\frac{1}{4\pi t}{\int }_{\partial B_t(x)}{P}_{r,\theta }^{(0)}(y)\mathrm{d}s(y)\right)=\frac{{\partial }^2}{\partial t^2}{\int }_{B_t(x)}\frac{P_{r,\theta }^{(0)}(y)}{4\pi |y-x|}\mathrm{d}y,$$

$$P_{r,\theta ,1}(t,x)$$

$$={\int }_{B_t(x)}\frac{1}{4\pi |y-x|}\left({\alpha }_1(y)\Delta P_{r,\theta ,0}(t-|y-x|,y)-\nabla {\rho }_1(y)\cdot \nabla P_{r,\theta ,0}(t-|y-x|,y)\right)\mathrm{d}y.$$

$${\int }_{B_t(x)}$$

$$=\frac{\partial }{\partial t}\left({\int }_{B_t(x)}\frac{{\alpha }_1(y){\partial }_t{P}_{r,\theta ,0}(t-|y-x|,y)}{4\pi |y-x|}\mathrm{d}y\right)$$

$$=\frac{{\partial }^2}{\partial t^2}\left({\int }_{B_t(x)}\frac{{\alpha }_1(y)P_{r,\theta ,0}(t-|y-x|,y)}{4\pi |y-x|}\mathrm{d}y\right)-\frac{\partial }{\partial t}\left({\int }_{\partial B_t(x)}\frac{{\alpha }_1(y)P_{r,\theta }^{(0)}(y)}{4\pi |y-x|}\mathrm{d}s(y)\right).$$

$$P_{r,\theta ,0}={\partial }_{tt}{\tilde{P}}_{r,\theta ,0}\text{with}{\tilde{P}}_{r,\theta ,0}(t,x)={\int }_0^t{\int }_0^{\tau }{P}_{r,\theta ,0}({\tau }_1,x)\mathrm{d}{\tau }_1\mathrm{d}\tau ={\int }_{B_t(x)}\frac{P_{r,\theta }^{(0)}(y)}{4\pi |y-x|}\mathrm{d}y.$$

$${\int }_{B_t(x)}\frac{{\alpha }_1(y)P_{r,\theta ,0}(t-|y-x|,y)}{4\pi |y-x|}\mathrm{d}y=\frac{{\partial }^2}{\partial t^2}\left({\int }_{B_t(x)}\frac{{\alpha }_1(y){\tilde{P}}_{r,\theta ,0}(t-|y-x|,y)}{4\pi |y-x|}\mathrm{d}y\right).$$

$$\begin{gathered} {\int }_{B_t(x)}\frac{\nabla {\rho }_1(y)\cdot \nabla P_{r,\theta ,0}(t-|y-x|,y)}{4\pi |y-x|}\mathrm{d}y \\ =\frac{{\partial }^2}{\partial t^2}\left({\int }_{B_t(x)}\frac{1}{16{\pi }^2|y-x|}{\int }_{B_{t-|y-x|}(y)}\nabla {\rho }_1(y)\cdot \frac{\nabla P_{r,\theta }^{(0)}(z)}{|z-y|}\mathrm{d}z\mathrm{d}y\right). \end{gathered}$$

Full text

[TABLE]

[TABLE]

1. Introduction

Photoacoustic imaging (PAI) is a novel imaging technique which uses pulsed laser excitation in the visible or infrared frequency regime to illuminate a specimen and measures the acoustic response of the medium (see [26] and some mathematical survey references [16, 25, 17]). Mathematical models of standard PAI do not take into account variable sound speed and elastic material parameters of the medium, such as compressibility and mass density [26]. More sophisticated models take into account spatially varying sound speed, but assume the sound speed to be known from measurements of other modalities (like for instance elastography and ultrasound experiments). PAI which takes into account given spatially varying sound speed and elasticity parameters has been considered in [1, 14, 13, 22, 24, 20, 18, 9, 10]. Little is known, how and when it is (at least) theoretically possible, to recover in addition to the photoacoustic imaging parameter, the absorption density, also sound speed and elastic parameters. [15] provided a method for parallel estimation of the sound speed and the photoacoustic imaging parameter via focusing to planes, which is called single plane illumination microscopy for a Born-approximation of the wave equation. This paper is the starting point for the present one, which shows that in addition to the sound speed, also elastic parameters can be recovered. Parallel sound speed and photoacoustic imaging parameter estimation has been considered as a nonlinear inverse problem (without making the assumption of a Born-approximation) in [21, 23, 19].

The main goal of this paper is to perform quantitative imaging in the acoustic regime: The photoacoustic imaging process consists of three processes, the acoustic wave propagation, the optical illumination and the visco-elastic part to transform optical energy into acoustic waves. Thus our topic might not be confused with quantitative imaging in the optical illumination part, where optical parameters of the medium are identified (see e.g. [6, 5, 7, 3, 4, 2, 8]).

In the following we describe the proposed experimental setup. We consider the specimen to be photoacoustically imaged to be supported in a bounded domain $\Omega_{0}\subset\mathds{R}^{3}$ embedded into a known environment. To analyse the specimen’s interior structure via photoacoustic imaging, it is illuminated with a laser pulse. This produces a local light fluence $\Phi_{r,\theta}:\mathds{R}^{3}\to[0,\infty)$ (integrated over the short time interval of the pulse), where we want to assume that the laser can be tuned in such a way that the function $\Phi_{r,\theta}$ is mainly supported in a small vicinity of the plane

[TABLE]

for some parameters $r\in[0,\infty)$ and $\theta\in\mathbf{S}^{2}$. Ideally, we are illuminating only single planes.

According to the photoacoustic effect, the object will absorb at each point parts of the light depending on the spatially varying absorption coefficient $\mu:\mathds{R}^{3}\to[0,\infty)$ and transform the absorbed energy first into heat, which is then transformed proportional to the Grüneisen parameter $\gamma:\mathds{R}^{3}\to[0,\infty)$ into a locally varying pressure distribution

[TABLE]

To avoid getting signals from everywhere, we want to assume that the absorption is confined to some bounded domain $\Omega\supset\overline{\Omega_{0}}$, that is, $\operatorname{supp}\mu\subset\Omega$.

This initial pressure distribution $P^{(0)}$ will then propagate as an elastic wave through the medium. To describe this wave, we model the object as an elastic body with smooth bulk modulus $K:\mathds{R}^{3}\rightarrow[0,\infty)$ and with a vanishing shear modulus. Then we know from linear elasticity theory that the stress tensor $\sigma_{r,\theta}\in C^{2}([0,\infty)\times\mathds{R}^{3};\mathds{R}^{3\times 3})$ as a function of time and space has the simple form

[TABLE]

where $u\in C^{3}([0,\infty)\times\mathds{R}^{3};\mathds{R}^{3})$ denotes the displacement vector of the material as a function of time and space. Introducing the mass density $\rho\in C^{1}(\mathds{R}^{3})$ of the material, this leads to the equation of motion:

[TABLE]

Moreover, the relation between the displacement vector $u_{r,\theta}$ and the pressure distribution $P_{r,\theta}\in C^{2}([0,\infty)\times\mathds{R}^{3})$ is given by

[TABLE]

To rewrite the equation of motion in terms of the pressure, we take the divergence of (1.1) and obtain

[TABLE]

Using again the equation (1.1), we can rewrite the term $\partial_{tt}u_{r,\theta}$ on the left hand side as an expression containing $u_{r,\theta}$ only in the form of $\operatorname{div}_{x}u_{r,\theta}$ and find with the relation (1.2) the equation of motion for the pressure $P_{r,\theta}$:

[TABLE]

Introducing the parameter $\alpha(x)=\frac{K(x)}{\rho(x)}$ corresponding to the square of the speed of sound, this can be written in the form

[TABLE]

Our aim is to reconstruct the product $\gamma\mu$ of absorption coefficient and Grüneisen parameter, the speed of sound $\sqrt{\alpha}$, and the mass density $\rho$ from measurements of this acoustic wave $P_{r,\theta}$. To accomplish this, we assume that we use the measurements

[TABLE]

for all points $\xi\in\partial\Omega$, all times $t\in[0,\infty)$, and all choices of illumination planes, that is, for all $r\in[0,\infty)$ and $\theta\in\mathbf{S}^{2}$.

In section 2 we formulate the forward problem, taking into account small perturbations of the acoustic parameters. We consider an asymptotic expansion of the solution of the acoustic wave equation induced by the deviation of sound speed and density around known constant functions. Proposition 2 combines the zeroth and first order term of the expansion as an integral operator acting on the initial acoustic pressure wave. The kernel of this operator is of utmost importance as it contains information about the seeked acoustic parameters. In section 3 we exploit asymptotic behaviors of this kernel which will be the main tools for the final reconstruction of the parameters and the source. In the final two section 4 and section 5 we consider the idealized forward problem where we are omitting the higher order terms of the expansion derived in section 2, which is reasonable when the distortions of the acoustic parameters become sufficiently small. The SPIM method then allows us to extract the product of the kernel and the source by applying the inverse Radon transform. By combining various measurement data, we are able to reconstruct the source and the distortions of the sound speed and density.

2. First Order Approximation

To simplify the problem, we will only consider the equation of motion (1.3a) in the case where $\rho$ and $\alpha$ are small perturbations of constant functions. To formulate this more precisely, let us introduce an artificial parameter $\varepsilon>0$ and assume that we have a family of media with parameters

[TABLE]

with $0<\rho_{\varepsilon,\text{min}}\leq\rho(x;\varepsilon)\leq\rho_{\varepsilon,\text{max}}$ and $0<\alpha_{\varepsilon,\text{min}}\leq\alpha(x;\varepsilon)\leq\alpha_{\varepsilon,\text{max}}$. The perturbations will be considered sufficiently smooth and supported in the domain of absorption, that is, $\operatorname{supp}(\rho_{0}-\rho)\subset\Omega$ and $\operatorname{supp}(\alpha_{0}-\alpha)\subset\Omega$.

Then, we can also expand the solution $P_{r,\theta}(\cdot,\cdot;\varepsilon)$ of the initial value problem (1.3) in $\varepsilon$:

[TABLE]

where the zeroth and first order terms fulfil

[TABLE]

and

[TABLE]

The higher order term satisfies the equation

[TABLE]

where

[TABLE]

sheds light on the coupling with the zeroth and first order expansion. The parameters $\rho_{0}$ and $\alpha_{0}$ are assumed to be known and we normalise them to $\rho_{0}=1$ and $\alpha_{0}=1$ henceforth. The initial value problems (2.3) and (2.4) can be explicitly solved, see for example [11, Section 2.4, Theorem 2], which gives us the following representation for $P_{r,\theta}$.

Proposition \theproposition.

Let $\alpha_{1}$, $\rho_{1}$ and $P_{r,\theta}^{(0)}$ be smooth functions. Then, the solution $P_{r,\theta}(\cdot,\cdot;\varepsilon)$ of the initial value problem (1.3) with the parameters $\rho$ and $\alpha$ given by (2.1) and (2.2) has the form

[TABLE]

where the kernel $K$ is given by

[TABLE]

where

[TABLE]

is the prolate spheroid with focal points $x$ and $y$ and larger semi-axis $\frac{t}{2}$, and $P_{r,\theta,2}$ remains bounded,

[TABLE]

with $C_{k,T;\varepsilon}=\mathcal{O}(1)$ as $\varepsilon\rightarrow 0$.

Proof*:*

According to [11, Section 2.4, Theorem 2], the solution $P_{r,\theta,0}$ of the initial value problem (2.3) is given by

[TABLE]

where we used the coarea formula to write the surface integral over the sphere as the derivative of the integral over the ball.

Similarly, the solution $P_{r,\theta,1}$ of (2.4) is explicitly known, see for example [11, Section 2.4, Theorem 4], and we have

[TABLE]

To write this as an integral operator on $P_{r,\theta}^{(0)}$, we integrate by parts to pull all derivatives acting on $P_{r,\theta,0}$ to the other terms.

• •

For the first term, we use that $P_{r,\theta,0}$ is a solution of (2.3) and get

[TABLE]

To reduce $P_{r,\theta,0}$ in the first summand to a term involving only the initial data $P_{r,\theta}^{(0)}$, we write, according to (2.9),

[TABLE]

Then, we pull the two time derivatives outside the integral and find because of $\partial_{t}\tilde{P}_{r,\theta,0}(0,x)=\tilde{P}_{r,\theta,0}(0,x)=0$ for the first integral in (2.11) that

[TABLE]

• •

For the second term in (2.10), we write $P_{r,\theta,0}$ in the form (2.12) and pull the time derivatives out of the integral to get

[TABLE]

Applying the divergence theorem to the inner integral, we obtain that

[TABLE]

We now insert (2.12) into (2.13) and this into (2.11). Moreover, we combine (2.15) with (2.14). The results, we then put into (2.10) and obtain the expression

[TABLE]

Combining the formulas (2.9) and (2.16) and replacing the surface integrals again with the coarea formula with time derivatives of integrals over the corresponding balls, we find that

[TABLE]

Remarking that we can explicitly calculate the primitive functions with respect to time of these integrals, for example:

[TABLE]

we integrate the result four times with respect to time to get rid of the time derivatives and then interchange the order of integration. This yields the expression

[TABLE]

with the kernel $K$ given by (2.6).

For the estimate of the higher order term we refer to [12, Section 23.2, Lemma 23.2.1]. Note that the coefficients of the homogeneous part of equation (2.5) depend continuously on $\varepsilon$. In the limit case $\varepsilon=0$, equation (2.5) results in the non-homogeneous wave equation with constant sound speed. $\square$

3. Properties of the Kernel

We want to study in this section the kernel $K$ defined in (2.6). In particular, we are interested in the behaviour of $K(t,x,y;\varepsilon)$ when $t$ gets much larger than the support of the involved functions $\alpha_{1}$ and $\rho_{1}$ and in the limit when $t$ approaches the distance $|y-x|$ from above (the kernel vanishes for $t<|y-x|$).

For large values of $t$, the change in the domain $\mathscr{E}_{t}(x,y)$ as $t$ varies does not change the value of the integral and $K$ behaves as a quadratic polynomial in $t$. Explicitly, we get with this argument the following expansion in $t$.

Lemma 3.1**.**

Let $\alpha_{1}\in C_{\mathrm{c}}(\mathds{R}^{3})$, $\rho_{1}\in C_{\mathrm{c}}^{1}(\mathds{R}^{3})$, and $K$ be defined by (2.6). For arbitrary values $x,y\in\mathds{R}^{3}$, let $T_{x,y}>0$ be chosen so that $\mathscr{E}_{T_{x,y}}(x,y)\supset\operatorname{supp}\rho_{1}\cup\operatorname{supp}\alpha_{1}$.

Then, we have for all $t\geq T_{x,y}$ the relation

[TABLE]

Proof*:*

We replace in (2.6) the integration over $\mathscr{E}_{t}(x,y)$ with the integral over $\mathds{R}^{3}$, which does not change the value of the integral, since we have $\operatorname{supp}\rho_{1}\cup\operatorname{supp}\alpha_{1}\subset\mathscr{E}_{T_{x,y}}(x,y)\subset\mathscr{E}_{t}(x,y)$, and reorder the terms as coefficients of a polynomial in $t$. This yields

[TABLE]

To simplify the coefficient of the linear term, we use the divergence theorem and that $\frac{1}{4\pi|x|}$ is the fundamental solution for the negative Laplace opertor $-\Delta$ and find

[TABLE]

Plugging this into (3.2), we arrive at (3.1). $\square$

The other limit, where the integral over the spheroid $\mathscr{E}_{t}(x,y)$ simplifies is $t\downarrow|y-x|$. In this case, the spheroid shrinks to the line from $x$ to $y$.

Lemma 3.2**.**

Let $x,y\in\mathds{R}^{3}$ be two arbitrary, different points and $\psi\in C^{2}(\mathds{R}^{3})$. Then, the function $F:[|y-x|,\infty)\to\mathds{R}$ defined by

[TABLE]

with $\mathscr{E}_{t}(x,y)=\{z\in\mathds{R}^{3}\mid|z-x|+|z-y|\leq t\}$ as before, is two times differentiable and can be expanded around $t=|y-x|$ as

[TABLE]

Proof*:*

We introduce the new coordinates $r_{1}\geq 0$, $r_{2}\geq 0$, and $\varphi\in[0,2\pi]$ of a point $\phi(r_{1},r_{2},\varphi)\in\mathds{R}^{3}$ such that

[TABLE]

To this end, let $e_{1}=\frac{y-x}{|y-x|}$, $e_{2}$, $e_{3}$ be an orthonormal basis of $\mathds{R}^{3}$ and set

[TABLE]

with the functions $\xi$ and $\eta$ explicitly given by

[TABLE]

see Figure 1.

The volume element by switching to the new coordinates $r_{1}$, $r_{2}$, $\varphi$ is given by

[TABLE]

Using the relations

[TABLE]

following directly from (3.4), we find that

[TABLE]

Thus, switching to these new coordinates (and remembering that a point with distance $r_{1}$ to $x$ and $r_{2}$ to $y$ only exists if the triangle inequalities $||y-x|-r_{1}|\leq r_{2}\leq|y-x|+r_{1}$ are fulfilled), the integral in (3.3) becomes

[TABLE]

To expand $F$ now in a Taylor polynomial around $t=|y-x|$, we first realise that $F(|y-x|)=0$ and then calculate the derivatives at $t=|y-x|$. For the first derivative, we find (the inequality $||y-x|-r_{1}|\leq t-r_{1}\leq|x-y|+r_{1}$ is equivalent to $\frac{t-|y-x|}{2}\leq r_{1}\leq\frac{t+|y-x|}{2}$) that

[TABLE]

Since $\phi(r_{1},|y-x|-r_{1},\varphi)=x+r_{1}e_{1}$ for all $r_{1}\in[0,|y-x|]$ and all $\varphi\in[0,2\pi]$, we get

[TABLE]

where $L_{x,y}$ denotes the straight line between the points $x$ and $y$.

For the second derivative, we obtain

[TABLE]

To calculate the limit $t\downarrow|y-x|$, we remark that

[TABLE]

where we introduced the manifold $D(t,r_{1})$ with boundary given as the disc normal to $e_{1}$ with center in $\frac{1}{2}(x+y)+\xi(r_{1},t-r_{1})e_{1}$ and radius $\eta(r_{1},t-r_{1})$, whose boundary $\partial D(t,r_{1})$ is parametrised by $\phi(r_{1},t-r_{1},\cdot)$ and its unit normal vector field $\nu$ is given by $\nu(\phi(r_{1},t-r_{1},\varphi))=\cos\varphi e_{2}+\sin\varphi e_{3}$. Now, using Stokes’ theorem, we can rewrite the boundary integral to an integral over the two-dimensional disc and find

[TABLE]

Plugging all this in our formula for $F^{\prime\prime}$, we can take the limit $t\downarrow|y-x|$ and get with

[TABLE]

that

[TABLE]

Since an integration by parts gives us

[TABLE]

all but the first term in the integral cancel each other, and the expression simplifies to

[TABLE]

$\square$

Corollary \thecorollary.

Let $\alpha_{1}\in C^{2}(\mathds{R}^{3})$, $\rho_{1}\in C^{1}(\mathds{R}^{3})$ and $K(t,x,y;\varepsilon)$ be given by (2.6). Then

[TABLE]

and

[TABLE]

Proof*:*

The first identity follows immediately by applying the previous lemma on the kernel. For the second identity, note that

[TABLE]

$\square$

4. Forward Problem

We will now consider the idealised problem where the light illumination $\Phi_{r,\theta}$ of the laser pulse is perfectly focused on the plane $E_{r,\theta}$ and where the perturbations of the parameters $\rho$ and $\alpha$ are so small that we can ignore the error term $P_{r,\theta,2}$ for the measurements. Henceforth we shall assume that we are given the measurements

[TABLE]

with $f(y)=\Phi^{(0)}\gamma(y)\mu(y)$ and $\mathcal{R}_{3}$ being the three-dimensional Radon transform. Note that the kernel vanishes for $|y-x|>t$.

For the following discussion, we shall assume that we are given the measurement data $M_{r,\theta}(t,x)$ for all $r\geq 0$, $\theta\in{\mathbf{S}}^{2}$, and $t>0$ at distinct points $x\in\Sigma$, where $\Sigma$ is a detector surface enclosing the compact support of $f$. We assume the seeked functions to be sufficiently smooth and that the perturbation of the sound speed and density are supported in a strict subset of $\operatorname{int}(\operatorname{supp}f)$.

5. Reconstruction Algorithm

Under the assumptions of the previous section, we are going to present a reconstrucion method for the functions $f$, $\alpha_{1}$ and $\rho_{1}$. Given our measurement data, we define the functions

[TABLE]

where $x\in\Sigma$ and $y\in\mathbb{R}^{3}$. Then, by virtue of Lemma 3.1 and Lemma 3.2 together with integration by parts, we obtain

[TABLE]

In the first step, we are going to reconstruct $\alpha_{1}$. Equation (5.4) gives us the values of $f$ on $\operatorname{supp}f\setminus(\operatorname{supp}\alpha_{1}\cup\operatorname{supp}\rho_{1})$. By knowledge of $f$ in a vicinity close to the boundary of $\operatorname{supp}f$ together with equation (5.1), we infer to know all line integrals of $\alpha_{1}$.

Proposition \theproposition.

*Assume $\operatorname{supp}\alpha_{1}\cup\operatorname{supp}\rho_{1}\subsetneq\operatorname{int}(\operatorname{supp}f)$ and suppose that the data $M_{0,x}$ and $N_{\infty,x}$ is given for all $x\in\Sigma$. Then the three-dimensional $X$-ray transform $X_{3}[\varepsilon\alpha_{1}](y,v)$ is known for all $y\in\mathbb{R}^{3}$ and $v\in\mathbf{S}^{2}$. *

Proof*:*

Let $[0,\infty)\ni\lambda\mapsto\gamma_{v}(\lambda)=x+\lambda v$, $v=(y-x)/|y-x|$, be the parametrisation of the half-ray connecting $x\in\Sigma$ with a point $y\in\operatorname{supp}f$. Thus, $y=\gamma_{v}(\mu)$ for some $\mu>0$ and the line integral in (5.1) becomes

[TABLE]

Choose $x^{\prime}\in\Sigma$ such that $y\in L_{x,x^{\prime}}$. Then

[TABLE]

Let $y_{0}\in L_{x,x^{\prime}}\cap\left(\operatorname{supp}f\setminus\left(\operatorname{supp}\alpha_{1}\cup\operatorname{supp}\rho_{1}\right)\right)$ for which $f(y_{0})\neq 0$ is known. Division through $f(y_{0})$ gives us

[TABLE]

$\square$

By inverting the $X$-ray transform, we obtain $\alpha_{1}$.

Proposition \theproposition.

For $\alpha_{1}$ sufficiently small in the $C^{2}$-norm we get

[TABLE]

Proof*:*

Combining (5.2) with (5.4), we get

[TABLE]

and the first identity follows. The distortion of the density can be expressed by virtue of (5.4). $\square$

Acknowledgement

The work of AB and OS has been supported by the Austrian Science Fund (FWF), Project P26687-N25 (Interdisciplinary Coupled Physics Imaging).