Photo-acoustic tomography in a rotating setting

TL;DR

This paper develops a direct inversion method for photo-acoustic tomography in a rotating setup, demonstrating unique recovery and optimal stability of optical parameters without relying on traditional two-step procedures.

Contribution

It introduces a novel direct inversion approach for rotating photo-acoustic tomography, modeling the process as a composition of pseudo-differential and Fourier integral operators, and proves injectivity and stability.

Findings

Unique reconstruction of optical coefficients in rotating setting

Establishment of optimal stability estimates

Extension of standard non-rotating results to rotating configurations

Abstract

Photo-acoustic tomography is a coupled-physics (hybrid) medical imaging modality that aims to reconstruct optical parameters in biological tissues from ultrasound measurements. As propagating light gets partially absorbed, the resulting thermal expansion generates minute ultrasonic signals (the photo-acoustic effect) that are measured at the boundary of a domain of interest. Standard inversion procedures first reconstruct the source of radiation by an inverse ultrasound (boundary) problem and second describe the optical parameters from internal information obtained in the first step. This paper considers the rotating experimental setting. Light emission and ultrasound measurements are fixed on a rotating gantry, resulting in a rotation-dependent source of ultrasound. The two-step procedure we just mentioned does not apply. Instead, we propose an inversion that directly aims to reconstruct the optical parameters quantitatively. The mapping from the unknown (absorption and diffusion) coefficients to the ultrasound measurement via the unknown ultrasound source is modeled as a composition of a pseudo-differential operator and a Fourier integral operator. We show that for appropriate choices of optical illuminations, the above composition is an elliptic Fourier integral operator. Under the assumption that the coefficients are unknown on a sufficiently small domain, we derive from this a (global) injectivity result (measurements uniquely characterize our coefficients) combined with an optimal stability estimate. The latter is the same as that obtained in the standard (non-rotating experimental) setting.

Equations

$$\begin{array}{r}-\nabla \cdot D\nabla u+\sigma u=0\text{in }\Omega ,\\ u=g\text{on }\partial \Omega .\end{array}$$

$$H(x)=\mu (x)\sigma (x)u(x),$$

$$\begin{array}{rl}& ({\partial }_t^2-c^2(x)\Delta )v(t,x)=0,\text{in }(0,\infty )\times {\mathbb{R}}^d,\\ & v(0,x)=H(x),x\in \Omega ,\\ & {\partial }_{t}v(0,x)=0,x\in {\mathbb{R}}^d,\end{array}$$

$$\begin{array}{rl}& ({\partial }_t^2-c^2(x)\Delta )v(t,x)=0,\text{in }(0,\infty )\times {\mathbb{R}}^d,\\ & v(0,x)=f(x),x\in \Omega ,\\ & {\partial }_{t}v(0,x)=0,x\in {\mathbb{R}}^d,\end{array}$$

$$\Lambda f=v(t,x)\big{|}_{[0,T]\times\partial\Omega_{M}}.$$

$$\begin{array}{rl}& \dot{X}(t)=\frac{c(X(t))}{2}\xi (t),\\ & \dot{\xi }(t)=-\frac{\nabla c(X(t))}{2}|\xi (t)|,\\ & X(0)=x_0,\xi (0)={\xi }_0.\end{array}$$

$$v(t,x)=\Lambda f(t,x)=(2\pi {)}^{-d}\sum _{\tau =\pm }\int {\mathrm{e}}^{i{\varphi }_{\tau }(t,x,\xi )-i\xi \cdot y}{a}_{\tau }(t,x,\xi )\hat{f}(\xi )\mathrm{d}\xi +Rf(t,x)$$

$$\mp {\partial }_t{\varphi }_{\pm }=c(x)|{\nabla }_x{\varphi }_{\pm }|,{\varphi }_{\pm }(0,x)=x\cdot \xi ,$$

$$\hat{f}(\xi )={\int }_{{\mathbb{R}}^d}{\mathrm{e}}^{-i\xi \cdot y}f(y)\mathrm{d}y.$$

$$\mathcal{G}=\{(t,x):x\in \Gamma ,0<t<s(x)\},$$

$${\tau }_{\pm }(x,\xi )=\max \{t\ge 0:{\gamma }_{x,\xi }(\pm t)\in {\Omega }_M\},$$

$$V_i={\chi }_i\Lambda f.$$

$$\begin{array}{r}-\nabla \cdot D\nabla u_{i,j}+\sigma u_{i,j}=0\text{in }\Omega ,\\ u_{i,j}=g_{i,j}\text{on }\partial \Omega ,j=1,2,...,N.\end{array}$$

$$H_{i,j}(x)=\sigma (x)u_{i,j}(x),$$

$$V_{i,j}={\chi }_i\Lambda H_{i,j},\text{for }i=1,2,...,M,j=1,2,...,N.$$

$$\delta V_{i,j}=V_{i,j}-{\overline{V}}_{i,j}={\chi }_i\Lambda \delta H_{i,j},\text{for }i=1,2,...,M,j=1,2,...,N.$$

$$\Vert \delta D{\Vert }_{L^2(\Omega )}+\Vert \delta \sigma {\Vert }_{H^1(\Omega )}\le C\left(\sum _{i=1}^M\sum _{j=1}^{2d}\Vert \delta V_{i,j}{\Vert }_{H^1([0,T]\times \partial {\Omega }_M)}\right).$$

$$\begin{array}{rl}-\nabla \cdot D\nabla \delta u_{i,j}+\sigma \delta u_{i,j}& =\nabla \cdot \delta D\nabla {\bar{u}}_{i,j}-\delta \sigma {\bar{u}}_{i,j}\text{in }\Omega ,\\ \delta u_{i,j}& =0\text{on }\partial \Omega .\end{array}$$

$$P_l(x,\partial )$$

$$P_r^{i,j}(x,\partial )$$

$$p_l(x,\xi )$$

$$p_r^{i,j}(x,\xi )$$

$$|D_x^{\alpha }{D}_{\xi }^{\beta }p(x,\xi )|\le C_{\alpha ,\beta ,\Omega }(1+|\xi |{)}^{m-|\beta |},\forall (x,\xi )\in \Omega \times {\mathbb{R}}^d.$$

$$\delta u_{i,j}=Q{P}_r^{i,j}\delta D-Q{\bar{u}}_{i,j}\delta \sigma .$$

$$r^{i,j}\sim \sum _{\alpha \in {\mathbb{N}}_0^d}\frac{1}{\alpha !}{\partial }_{\xi }^{\alpha }q{D}_x^{\alpha }{p}_r^{i,j},\text{where }{D}_x^{\alpha }=(-i{)}^{|\alpha |}{\partial }_x^{\alpha },$$

$$r_0^{i,j}(x,\xi )=\frac{-i\nabla {\bar{u}}_{i,j}(x)\cdot \xi }{{\xi }^2+1}\text{ mod }{S}^{-2}(\Omega ,{\mathbb{R}}^n).$$

$$\delta H_{i,j}=\sigma \delta u_{i,j}+\delta \sigma {\bar{u}}_{i,j},$$

$$\delta H_{i,j}=\sigma (Q{P}_r^{i,j}\delta D-Q{\bar{u}}_{i,j}\delta \sigma )+{\bar{u}}_{i,j}\delta \sigma ,$$

$${\mathcal{H}}_D^{i,j}=\sigma Q{P}_r^{i,j},{\mathcal{H}}_{\sigma }^{i,j}={\bar{u}}_{i,j}(I-Q\sigma ),\text{and}{\mathcal{H}}^{i,j}=[{\mathcal{H}}_D^{i,j},{\mathcal{H}}_{\sigma }^{i,j}].$$

$$\delta H_{i,j}={\mathcal{H}}^{i,j}\left[\begin{array}{c}\delta D\\ \delta \sigma \end{array}\right]={\mathcal{H}}_D^{i,j}\delta D+{\mathcal{H}}_{\sigma }^{i,j}\delta \sigma +S_{i,j}\delta p,$$

$$V_{i,j}={\chi }_i(2\pi {)}^{-d}\sum _{\tau =\pm }\int \int {\mathrm{e}}^{i{\varphi }_{\tau }(t,x,\xi )-i\xi \cdot y}{a}_{\tau }(t,x,\xi )H_{i,j}(y)\mathrm{d}y\mathrm{d}\xi +{\chi }_{i}R{H}_{i,j}.$$

Full text

This paper considers the rotating experimental setting. Light emission and ultrasound measurements are fixed on a rotating gantry, resulting in a rotation-dependent source of ultrasound. The two-step procedure we just mentioned does not apply. Instead, we propose an inversion that directly aims to reconstruct the optical parameters quantitatively. The mapping from the unknown (absorption and diffusion) coefficients to the ultrasound measurement via the unknown ultrasound source is modeled as a composition of a pseudo-differential operator and a Fourier integral operator. We show that for appropriate choices of optical illuminations, the above composition is an elliptic Fourier integral operator. Under the assumption that the coefficients are unknown on a sufficiently small domain, we derive from this a (global) injectivity result (measurements uniquely characterize our coefficients) combined with an optimal stability estimate. The latter is the same as that obtained in the standard (non-rotating experimental) setting.

††E-mail: guillaumebal@uchicago.edu and adrki@dtu.dk

1 Introduction

Photo-acoustic tomography (PAT) is a coupled-physics (also known as hybrid) imaging method that aims to reconstruct the optical parameters of biological tissues. The optical parameters are known to provide valuable, high-contrast information about for example cancerous tissue, and hence are of clinical interest [13, 34, 29, 12]. To image these parameters, PAT leverages the thermoelastic expansion generated by absorbed light. A domain of interest is illuminated by electromagnetic waves and the absorption of light within the sample causes a rapid expansion of the material proportional to the amount of absorbed photons; this conversion of light to ultrasound is what is known as the photo-acoustic effect. The expansion initiates a pressure wave, and this wave is recorded by acoustic detectors at the boundary of the sample. The objective of PAT is to reconstruct the optical parameters from such acoustic measurements.

A model approximating the propagation of light in the diffusive regime is the following:

[TABLE]

Here $\Omega\in\mathbb{R}^{d}$ (typically $d=2,3$) is the object domain, $u(x)$ is the light/photon intensity, $D(x)$ is the diffusion coefficient and $\sigma(x)$ is the absorption parameter. The illumination is modelled as the boundary condition $g$. The thermo-elastic expansion is proportional to the amount of absorbed light, and is given by

[TABLE]

where $\mu(x)$ is a proportionality coefficient known as the Grüneisen coefficient. It is known that the reconstruction of the three parameters $(\sigma,D,\mu)$ is not possible without prior assumptions or multi-color measurements [11, 8]. In the rest of the paper, we assume $\mu(x)$ known.

To model the acoustic wave resulting from the thermoelastic expansion, we consider the linear wave equation

[TABLE]

where $v(t,x)$ is the acoustic pressure and $c(x)$ is the sound speed, assumed to be known. For some domain $\Omega_{M}$ such that $\Omega\subset\Omega_{M}$, one measures $v(t,x)$ for $(t,x)\in[0,T]\times\partial\Omega_{M}$ for some sufficiently large duration $T>0$. We refer to this as the full data. The quantitative reconstruction of $D$ and $\sigma$ is then often considered in two successive steps:

The reconstruction of $H(x)$ from v(t,x)\big{|}_{[0,T]\times\partial\Omega_{M}}. In the case of full data, this is a well-posed problem and several inversion methods exist[20, 24, 31, 32]. 2. 2.

The reconstruction of $D(x)$ and $\sigma(x)$ from $H(x)$ (assuming $\mu(x)$ known). This is quantitative PAT, which is known to be reasonably well-posed [11, 12, 17, 24].

The fact that both steps are well-posed is what makes PAT an attractive medical imaging modality. However, the requirements for the decoupling into two separate steps can be hard to meet in an experimental situation. Having access to a full measurement of the acoustic wave requires the illuminated object to be fully surrounded by acoustic measurement devices, a situation which is hard to achieve in most situations of interest. Since Step 1 cannot be stably completed without access to (sufficiently) full data, one instead relies on doing several partial measurements, changing the position and the illumination patterns for each measurement, for instance by rotating both the light source and ultrasound detectors [13, 26]. In [10], such a situation, with $D$ constant, was considered, and it was shown that the absorption coefficient $\sigma$ can be stably determined in situations where the object (or measurement device) is rotating. In the current paper we investigate the situation were both $D$ and $\sigma$ are unknown, and the measurements and illuminations are rotating. We show that under certain assumptions on the parameters, this situation also allows for stable determination of both $D$ and $\sigma$.

Let us conclude this introductory section by mentioning that PAT is one of the many hybrid (coupled-physics) medical imaging modalities that have emerged in recent years. For a brief list of mathematical description of such methods, we refer the reader to, e.g., [4, 5, 6, 22, 28, 30, 33] and their multiple references. The rest of the paper is organized as follows. The setting and main uniqueness and stability result are presented in section 2. The details of the derivation are collected in section 3 while some concluding remarks on the imaging modality and the mathematical tools are given in section 4.

2 Background and main results

2.1 Sound propagation

The linear acoustic wave equation is a reasonable description of sound propagation [23, 31, 10]. The sound speed $c(x)$ is considered known in this paper. Taking $v$ to be the solution of

[TABLE]

we define the acoustic measurement operator $\Lambda$ as

[TABLE]

In the setting considered here, we do not have one, but rather a large number of ultrasound sources $f(x)=f_{i}(x)$ corresponding for instance to the rotating measurement setting, where a rotating set of light sources generates a different ultrasound source for each rotation. Since the ultrasound detectors also rotate in such a setting, $v(t,x)=\Lambda f$ is available only for $x$ in the support of the rotating detector, which is a small fraction of the boundary $\partial\Omega_{M}$. As a consequence, we cannot reconstruct the whole $f(x)$ stably from such measurements [4]. All that we can reconstruct is the singularities of $f$ that are visible from the available measurements. This mapping from the visible singularities of $f$ to the available measurements is described by a standard microlocal procedure, which we now recall [31].

Consider the Hamiltonian system:

[TABLE]

The curves $(X(t),\xi(t))$ are known as bi-characteristic curves, and each $X(t)$ is called a ray. We assume that $c(x)\in C^{\infty}(\mathbb{R}^{d})$ is non-trapping, that is $|X(t)|\to\infty$ as $t\to\infty$, for all rays [23]. An important feature of the rays is that they are curves along which the propagation of the singularities of a wave occurs, a fact that motivates the preceding assumptions on the measurement geometry, see [31, 10]. To that end, we use the unit speed geodesics $\gamma_{x_{0},\xi_{0}}(t)$, defined by the relation $\dot{\gamma}_{x_{0},\xi_{0}}(t)=\dot{X}(t)/|\dot{X}(t)|$, where $X(t)$ is the solution to (5) with initial conditions $(x_{0},\xi_{0})$.

We can then describe the measurement operator $\Lambda$ has the following Fourier integral operator(FIO) [15, 31]:

[TABLE]

for $(t,x)\in[0,T]\times\partial\Omega_{M}$, where $R$ is a linear operator with smooth Schwartz kernel. The phase functions $\phi_{\pm}$ are solutions to the eikonal equations

[TABLE]

and homogeneous of order one in $\xi$. As usual,

[TABLE]

Note that if $c=\text{const.}$, then $\phi_{\pm}=x\cdot\xi\pm c|\xi|t$, and the solution is exact. The function $a_{\pm}$ is called a classical amplitude of order 0, and satisfies a recursive transport equation; see ([15], page 128) for details.

We have assumed here that the FIO could be represented with a single (global) phase function $\phi_{\tau}$. This can always be done for sufficiently short times. For longer times, and to avoid the presence of caustics, the above operator should really be written as a composition of a finite number of such terms, or more generally as a globally defined FIO [15, 31]. To simplify notation, we represent our FIO, mapping all the singularities of $f$ (that is $f$ up to a smooth term) to the available measurements by (6).

2.2 Detector model

Assume a detector supported on $\Gamma_{d}$, a closed and bounded hypersurface in $\mathbb{R}^{d}$. The $i$’th measurement consists of acquiring $v(t,x)$ for $t\in[0,T]$ and $x\in\Gamma_{i}$, where $\Gamma_{i}$ is a translation and rotation of $\Gamma_{d}$.

To obtain enough information on the optical coefficients, we make the following assumption on the measurement setting. We require overlapping measurements, i.e., that $\mu(\Gamma_{i}\cap\left(\cup_{j\neq i}\Gamma_{j}\right))>0$, for all $i$, where $\mu$ is the Lebesgue (surface) measure on $\Gamma_{d}$. For a collection of $M$ measurements, we set the full measurement surface $\Gamma=\cup_{i=1}^{M}\Gamma_{i}$. We then assume we have a domain $\Omega_{M}$ such that $\Omega\subset\Omega_{M}$ and $\Gamma=\partial\Omega_{M}$. In the case of rotating measurements, we obtain that $\Gamma=\partial\Omega_{M}$ for some ball $\Omega_{M}=B_{r}$ of radius $r$ containing $\Omega$.

To describe our measurement setting, we follow the approach of [31] used to analyze the inversion of the wave equation with incomplete data. Our requirements on $\Gamma$ and $T$ are the same as those required to get stable inversions with incomplete data. Define

[TABLE]

where $s(x)$ is a continuous function determining the temporal measurement interval. We introduce the function $\tau_{\pm}(x,\xi)$, defined as

[TABLE]

where $\gamma_{x,\xi}(t)$ are the rays in (5). In order to make sure all of the singularities of the function $H(x)$ in (3) (and hence of the optical coefficients via (2)) are captured by the measurement, additional requirements on $\mathcal{G}$ (and hence $\Gamma$ and the measurement interval length $T$) are

• •

$\forall x\in\Omega,\quad\exists z\in\Gamma\text{ such that dist}(x,z)<s(z)$. Here dist$(x,z)$ is the length of the geodesic connecting $x$ and $z$, with respect to $\mathrm{d}s^{2}=c^{-2}(x)\mathrm{d}x^{2}$. The conditions says that all of the wave should reach $\Gamma$ during the measurement period.

• •

$\forall(x,\xi)\in T^{*}\Omega\backslash 0,\quad(\tau_{k}(x,\xi),\gamma_{x,\xi}(\tau_{k}(x,\xi))\in\mathcal{G},$ for $k=+$ or $k=-$.

Here $T^{*}\Omega=\Omega\times\mathbb{R}^{d}$ is the cotangent bundle of $\Omega$, where we consider $\xi\not=0$ only. The motivation for this condition is that every point in the wavefront set of $H$ should reach the measurement surface.

The above conditions are satisfied if for example $\Gamma$ is a sphere containing $\Omega$ and an interval such that $s(x)>\max_{x,z\in\bar{\Omega}_{M}}\text{dist}(x,z)/2$ when the wave speed is constant. These conditions are also satisfied for non-trapping speeds for sufficiently long measurement times [31].

Last, for the measurement at $\Gamma_{i}$, let $\{\varphi_{j}\}_{j=1}^{M}$ be a partition of unity subordinate to $\{\Gamma_{j}\}_{j=1}^{M}$. Note that there is a $K_{i}\subset\Gamma_{i}\setminus\left(\cup_{j\neq i}\Gamma_{j}\right)$ such that $\varphi_{i}|_{K_{i}}=1$, while $\varphi_{j}|_{K_{i}}=0$ for $j\neq i$. Next, we take $\psi\in C^{\infty}_{0}([0,2T])$ such that $\psi(t)=1$ for $t\in[\varepsilon,T]$, with $\varepsilon$ sufficiently small, and set $\chi_{i}(x,t)=\psi(t)\varphi_{i}(x)$. We now define our rotating measurement $V_{i}$ at $\Gamma_{i}$ by

[TABLE]

Note that $f\equiv H=\mu\sigma u$ above should really be read as $f=f_{i}=H_{i}=\mu\sigma u_{i}$, which is rotation dependent as the illumination $g=g_{i}$ in (1) rotates along with the detector $\Gamma_{i}$. There is therefore no hope to reconstruct all sources $f_{i}$ from measurements of the form (7) (with $f$ replaced by $f_{i}$) unless $\Gamma_{i}=\Gamma$. All the sources $f_{i}\equiv H_{i}$ have to be anchored to a rotation-independent object, namely the domain of interest modeled by the optical parameters $(D,\sigma)$. This is the objective of the next section.

2.3 Light propagation and boundary conditions

The second order elliptic PDE (1) serves as a reasonable model for propagation of light in in highly scattering media such as biological tissues.

We now aim to understand how the optical parameters $(D,\sigma)$ influence the ultrasound sources $H=\mu\sigma u$. Since $H\equiv H_{i}$ corresponds to rotating illuminations $g=g_{i}$ in (1), and rotating ultrasound measurements such as (7) do not allow full reconstructions of each $H_{i}$, we are also forced to understand such an influence locally, and in fact micro-locally. As in the derivation of (6), the resulting pseudo-differential calculus requires enough smoothness for all the required Taylor-type expansions to make sense. We therefore assume that $(D,\sigma)$ are smooth, as we did for the sound speed $c(x)$. A laborious, standard, tracking of all relevant calculations shows that only finitely many terms are necessary in each Taylor expansion. As a consequence, all results hold for $(D,\sigma,c)$ of class $C^{k}$ for $k$ sufficiently large. We will present all results assuming $k=\infty$ both to simplify and stress that in practice, the difference between large $k$ and yet larger $k$ is somewhat immaterial.

Suppose that the object of interest occupies an open domain $\Omega\in\mathbb{R}^{d},d=2,3$, where the boundary $\partial\Omega$ is $C^{\infty}$. Assume $D\in C^{\infty}(\mathbb{R}^{d})$, and such that there is a positive constant $C_{D}$ and $C_{D}^{-1}\leq D\leq C_{D}$, and that $\text{supp}(D-1)\subset\Omega$. For $0\leq\sigma\in C^{\infty}_{0}(\mathbb{R}^{d})$, assume there is some closed $\tilde{\Omega}$ such that $\Omega\subset\tilde{\Omega}$, and that $\text{supp}(\sigma)=\tilde{\Omega}$.

It is well-known in the non-rotating setting that stable reconstruction of both $D$ and $\sigma$ requires internal functionals $H(x)$ from multiple different illuminations, and that appropriate illuminations exist [12, 11]. Consider such a set of boundary conditions $\{g_{j}\}_{j=1}^{N}$. When the object is rotating and the detectors and illumination patterns are fixed, or equivalently, the detector and illumination patterns rotate, as for example in [14, 26], we get a different set of $N$ illuminations for each rotation $i=1,2,...M$, i.e., $\{g_{i,j}\}_{j=1}^{N}$. Now, for a fixed $i$, let $u_{i,j}$ be the solution to

[TABLE]

If the $g_{i,j}$ are smooth, we have that $u_{i,j}\in C^{\infty}(\Omega)$ ([16], Theorem 3, Chapter 6.3).

We already mentioned that, regardless of the number of illuminations, one cannot reconstruct both $D$,$\sigma$ and the Grüneisen coefficient $\mu$ [11]. To simplify notation, we assume $\mu$ known and set

[TABLE]

where $u_{i,j}$ is the solution to (8). Indexing the acoustic measurements accordingly, we have that

[TABLE]

This provides a full description of our measurement setting. The parameter $N$ provides the diversity in boundary illuminations that is necessary to stably reconstruct the optical parameters in the second step of standard quantitative PAT [12, 11]. The parameter $M$ indicates the number of rotations necessary to obtain measurements over all of $\Gamma$, and so roughly corresponds to the ratio between the size (volume of the hypersurface) of $\Gamma$ divided by that of the support of the rotating detector.

2.4 Main results

We can now state the main result of the paper. This is, under some restrictive assumptions, an injectivity and stability result. It does not provide a reconstruction algorithm, merely the reassurance that enough information has been collected to uniquely and stably characterize the unknown optical coefficients.

We therefore consider two pairs of admissible diffusion and absorption coefficients $(D,\sigma)$ and $(\bar{D},\bar{\sigma})$, satisfying the assumptions from Section 2.3. Since the PDO and FIO calculus we use here is not meant to handle boundaries, we have to assume that the values of the coefficients agree in a neighborhood of $\partial\Omega$. For fixed $i,j$, let $u_{i,j}$ and $\bar{u}_{i,j}$ be solutions to (8) with illumination $g_{i,j}$ and optical parameters $(D,\sigma)$, $(\bar{D},\bar{\sigma})$, respectively. Denote $\delta H_{i,j}=\sigma u_{i,j}-\bar{\sigma}\bar{u}_{i,j}$, $\delta D=D-\bar{D}$ and $\delta\sigma=\sigma-\bar{\sigma}$. We collect the measurements for each illumination and write their difference as

[TABLE]

We use here the fact that the wave propagation step is linear: the measurements are linear in the source terms $H_{i}$. The full inverse problem, which maps $(D,\sigma)$ to such measurements, is however nonlinear. In spite of this, we will show in the next section that the terms $\delta H_{i,j}$ can be written as a functional that is nonlinear in $(D,\sigma,\bar{D},\bar{\sigma},u_{i,j},{\bar{u}}_{i,j})$ but linear in $(\delta D,\delta\sigma)$, essentially as a generalization of the fact that for any polynomial $p(x)$, we can find another polynomial $q(x,y)$ such that $p(x)-p(y)=q(x,y)(x-y)$.

Moreover, we will show in the next section that there exist open sets of illuminations $\{g_{i,j}\}$ such that the mapping from $(\delta D,\delta\sigma)$ to $\{\delta H_{i,j}\}$ may be described as an elliptic pseudo-differential operator (PDO) with a symbol that depends on $(D,\sigma,\bar{D},\bar{\sigma},u_{i,j},{\bar{u}}_{i,j})$. The first result along these lines was obtained in the non-rotating setting (with $M\equiv 1$) in [24].

Combining the elliptic FIO in (6) with the above elliptic PDO shows that the mapping from $(\delta D,\delta\sigma)$ to $\delta V_{j}$ is itself an elliptic FIO with symbol that depends on the $(D,\sigma,\bar{D},\bar{\sigma},u_{i,j},{\bar{u}}_{i,j},c(x))$, which must all be sufficiently smooth for the calculus to apply.

Such micro-local results provide optimal stability estimates as well as an injectivity result provided that a smoothing compact operator does not have eigenvalue one. This assumption is notoriously difficult to verify, unless we have recourse to a smallness assumption somewhere. In [10], where the PAT problem with constant $D$ is considered, the smallness assumption was on the size of the absorption coefficient. Here, we make no assumption on the $O(1)$ size of the absorption and diffusion coefficients. Rather, we assume that the support of the domain where $\delta D=D-\bar{D}$ and $\delta\sigma=\sigma-\bar{\sigma}$ are unknown is itself sufficiently small. This will prove to be a sufficient assumption to obtain the following result.

Theorem 1**.**

If $\text{supp}(\sigma-\bar{\sigma})$ and $\text{supp}(D-\bar{D})$ are contained in sufficiently small ball $B_{\varepsilon}\subset\Omega$, there exist for each rotation $1\leq i\leq M$ an open set of $2d$ illuminations $\{g_{i,j}\}_{j=1}^{2d}$ such that the following estimate holds for the corresponding measurements.

[TABLE]

The constant $C$ depends on $D,\bar{D},\sigma,\bar{\sigma},c,\Omega$ and the illuminations.

Additional remarks on the set of necessary illuminations will be provided in the next section. We note that in dimension $d\geq 3$, no boundary conditions $g$ guarantee the necessary ellipticity assumptions to obtain the optimal stability estimates given in the above theorem independently of the coefficients $(\sigma,D)$ [3, 4].

3 Construction of elliptic operator

In lines with [24, 10], we show that the problem mapping the unknown coefficients to the (also unknown) ultrasound sources $\delta H_{i,j}$ is described by an elliptic pseudo-differential operator for an appropriate choice of illuminations and for smooth coefficients and solutions of (1).

3.1 Parametrix for $\delta H_{i,j}$

Let the situation be as described in Sections 2.3-2.4. Denote $\delta u_{i,j}=u_{i,j}-\bar{u}_{i,j}.$ Subtracting the equations (8) for $u_{i,j}$ and $\bar{u}_{i,j}$, we obtain an equation for $\delta u_{i,j}$ in terms of $\delta D$ and $\delta\sigma$:

[TABLE]

We define the operators

[TABLE]

with symbols

[TABLE]

Recall that a function $p(x,\xi)\in C^{\infty}(\Omega\times\mathbb{R}^{d})$ is said to be a symbol of class $S^{m}(\Omega\times\mathbb{R}^{d})$, $m\in\mathbb{R}$, if for any multi-indices $\alpha,\beta\in\mathbb{N}^{d}$, there is a constant $C_{\alpha,\beta,\Omega,}$ such that

[TABLE]

In particular, the symbols $p_{l}$ and $p_{r}^{i,j}$ are homogeneous in $\xi$, and their principal symbols are $p_{l,0}(x,\xi)=D(x)\xi^{2}$ and $p_{r,0}^{i,j}(x,\xi)=-i\nabla\bar{u}_{i,j}\cdot\xi$. Since $D(x)$ is positive on $\Omega$, $P_{l}$ is elliptic and $p_{l,0}(x,\xi)\sim\xi^{2}$. On the other hand, $p_{r,0}^{i,j}(x,\xi)$ is not elliptic at points $(x_{0},\xi_{0})$ such that $\nabla\bar{u}_{i,j}(x_{0})\cdot\xi_{0}=0$. Let $Q$ be the parametrix (an inverse modulo smooth terms) of $P_{l}$, with symbol $q$ and principal term $q_{0}$. Since $D(x),\sigma(x)>0$ in $\Omega$, we have that $q(x,\xi)=\frac{1}{\xi^{2}+1}\text{ mod }S^{-3}(\Omega,\mathbb{R}^{n})$, meaning that $q-\frac{1}{\xi^{2}+1}\in S^{-3}(\Omega,\mathbb{R}^{d}).$ Now we can solve (12) for $\delta u_{i,j}$, modulo smooth terms:

[TABLE]

The symbol $r^{i,j}$ of the composition $QP_{r}^{i,j}$ has the asymptotic expansion

[TABLE]

and hence the principal part is

[TABLE]

Next, consider $\delta H_{i,j}=H_{i,j}-\bar{H}_{i,j}=\sigma u_{i,j}-\bar{\sigma}\bar{u}_{i,j}$. Rearranging, we have

[TABLE]

and hence we can write $\delta H_{i,j}$ as the result of a pseudo-differential operator acting on $\delta D$ and $\delta\sigma$, i.e,

[TABLE]

modulo smooth terms. We define

[TABLE]

In terms of (13) we have

[TABLE]

where $S_{i,j}$ is a linear operator with smooth Schwartz kernel and $\delta p:=(\delta D,\delta\sigma)$. Let us finally note that the principal symbols of $\mathcal{H}^{i,j}_{D}$ and $\mathcal{H}^{i,j}_{\sigma}$ are $h^{i,j}_{D,0}(x,\xi)=r_{0}^{i,j}(x,\xi)$ and $h_{\sigma,0}^{i,j}(x,\xi)=\bar{u}_{i,j}(x)$, respectively.

3.2 Coupling with wave propagation

From Section 2.1 and (9), the measurement $V_{i,j}$ is given by

[TABLE]

Letting $V_{i,j}$ and $\bar{V}_{i,j}$ be solutions to (4) with initial conditions $H_{i,j}$ and $\bar{H}_{i,j}$ respectively, we set $\delta V_{i,j}=V_{i,j}-\bar{V}_{i,j}$. Due to the linearity of (4), we have

[TABLE]

The composition of an FIO with a PDO is well defined: from ([15], Theorem 4.2), the resulting operator is again a FIO with the same phase function and with a amplitude function with asymptotic expansion given by

[TABLE]

From [31], we know that $a_{\tau}$ is a zeroth order amplitude with $a_{\tau,0}=\frac{1}{2}$. Hence the principal term of (14) is $\frac{1}{2}(h_{D,0}^{i,j}(y,\xi)+h_{\sigma,0}^{i,j}(y,\xi))\text{ mod }S^{-2}(\Omega\times\mathbb{R}^{d}),\tau=\pm.$

The forward map now takes the form

[TABLE]

modulo a smooth term $Q_{i,j}\delta p$ involving a linear operator $Q_{i,j}$ with smooth Schwartz kernel.

3.3 Time-reversal

To bring us back from boundary measurements to objects (such as $\delta p$) that are defined on the spatial domain $\Omega$, we apply the time-reversal inversion to the composed operator in (15). The time reversal operator is an approximate inverse to $\Lambda$, which we denote by $A$ [31, 32, 20]. For given data $h\in H^{1}([0,T]\times\partial\Omega_{M})$, let $w$ be the solution of

[TABLE]

where $\phi$ solves $\Delta\phi=0,\phi|_{\partial\Omega_{M}}=h(T,\cdot)$. Then Ah=w\big{|}_{t=0,x\in\Omega}. For $f\in C^{\infty}_{0}(\Omega)$, we have that $f-A\Lambda f=Kf$, where $K$ is a compact operator with $\|K\|<1$. See [31] for details.

For a fixed measurement position $i$ and illuminations $g_{i,j}$, we now consider the system

[TABLE]

where we stacked the $\delta V_{i,j}$ and applied the time reversal operator to each side. We now define

[TABLE]

For a more compact notation in what follows, we write

[TABLE]

We will need the following lemma.

Lemma 1**.**

The operator $\kappa_{i}$ is a pseudo-differential operator with principal symbol $\kappa_{i}^{0}$ given by

[TABLE]

Proof.

It follows the linearity of $A$ and Theorem 3 in [31] that $A\left(\chi_{i}\Lambda\right)$ is a pseudo-differential operator of order zero with principal symbol $\tilde{\chi}_{i}(x,\xi)$. Since $A(\Lambda\mathcal{H}^{i,j})=(A\Lambda)\mathcal{H}^{i,j}$, we have that $A\left(\chi_{i}\Lambda\mathcal{H}^{i,j}\right)$ is again a pseudo-differential operator with symbol $\tilde{\chi}_{i}(x,\xi)(h^{i,j}_{D}(x,\xi)+h^{i,j}_{\sigma}(x,\xi))$. Considering the expressions for $h^{i,j}_{D}(x,\xi)$ and $h^{i,j}_{\sigma}(x,\xi)$ found in Section 3.1, equation (17) follows. ∎

By applying the adjoint of $\kappa_{i}$, $\kappa_{i}^{*}$ to (16) and summing over $i$ we get the system

[TABLE]

where

[TABLE]

Now $\delta\mathcal{V}$ contains all the measured data, and the linear operator $P$ has a smooth Schwartz kernel.

We thus have replaced our inversion problem by the analysis of equation (18).

3.4 Analysis of the system of pseudo-differential operators

We now investigate the system

[TABLE]

Such a system is called elliptic when the principal symbol $\sum_{i=1}^{M}(\kappa_{i}^{0})^{*}\kappa_{i}^{0}$ of the operator $\sum_{i=1}^{M}\kappa_{i}^{*}\kappa_{i}$ is full rank for all $(x,\xi)\in T^{*}\Omega$, [19, 2]. We show that there exist suitable sets of illuminations such that this is the case.

Since each of the terms in $\sum_{i=1}^{M}(\kappa_{i}^{0})^{*}\kappa_{i}^{0}$ is positive semi-definite, the sum will also be positive semi-definite. Hence it suffices to show that for each $(x,\xi)\in T^{*}\Omega$, at least one term is positive definite, as this will then guarantee the positive definiteness and full rank property of the sum. Since $(\kappa_{i}^{0})^{*}\kappa_{i}^{0}$ is positive when $\kappa_{i}^{0}$ is full rank, we aim to find illuminations such that for every $(x,\xi)\in T^{*}\Omega$, there is some $\kappa_{i}^{0}$ with full rank. Recall that

[TABLE]

From the assumptions on the non-trapping wave speed and sufficiently long measurement time, we know that for each $(x,\xi)\in T^{*}\Omega$, at least one $\tilde{\chi}_{i}(x,\xi)>0$. For $(x,\xi)\in T^{*}\Omega$ and such a rotation index $i=i(x,\xi)$, we need to show that $\kappa_{i}^{0}(x,\xi)$ is full rank.

We denote the determinant of $2\times 2$ matrix consisting of rows $m$ and $n$ of (21) by $q_{m,n}(x,\xi)$. Then

[TABLE]

For $\kappa_{i}^{0}(x,\xi)$ to be full-rank, we thus require that there is at least one $q_{m,n}(x,\xi)\neq 0$ for every $(x,\xi)\in\Omega\times S^{d-1}$. Here $S^{d-1}$ denotes the unit sphere in $\mathbb{R}^{d}$, and replaces $\mathbb{R}^{d}\setminus\{0\}$ since only the direction of $\xi$ is of importance.

Hence we seek to find illuminations $\{g_{i,j}\}$ such that for some indices $(m,n)$ it holds that

[TABLE]

Expressions similar to (23) appears in the literature on PAT, e.g., in [12, 11, 24, 4].

For (23) to hold, we need that the vector fields

[TABLE]

constitute a basis for $\mathbb{R}^{d}$ for each $x\in\Omega$. In that way, a direction $\xi_{0}\in S^{d-1}$ can never be orthogonal to all $v_{m,n}$. It is therefore clear that we must have at least $d+1$ different illuminations $\{g_{i,j}\}_{j=1}^{d+1}$. We will use complex geometric optics (CGO) solutions to show that there exists a set of $2d$ illuminations such that the vector fields $v_{m,n}$ does indeed form a basis for $\mathbb{R}^{d}$ for each $x\in\Omega$, and also show that under certain restrictions on the optical coefficients, the same boundary conditions will work for all rotations. We briefly introduce the CGO solutions first.

By the Liouville transformation $v=\sqrt{D}u$, equation (8) is written as

[TABLE]

where $q(x)=\Delta\sqrt{D}/\sqrt{D}+\sigma/D$. The CGO solutions are special solutions to (24) that are perturbations of complex plane waves $\mathrm{e}^{\rho\cdot x}$ of the form

[TABLE]

where $\rho\in\mathbb{C}^{d}$ and $\rho\cdot\rho=0$. For $|\rho|$ large enough, such solutions exist, and from ([12], Corollary 3.2), it is known that the perturbation term $\psi_{\rho}$ satisfies the bound

[TABLE]

where $s>d/2+k$ for $k\geq 1$. By choosing the parameter $\rho$ in a certain way, we can tailor these solutions to achieve solutions with the right properties. We present the result for the situation when $d=3$.

Theorem 2**.**

For any rotation $i$ there is a open set of illuminations $\{g_{i,j}\}_{j=1}^{6}$ such that for each $x\in\Omega$ the vector fields

[TABLE]

*form a basis for $\mathbb{R}^{3}$.

In addition, if $D$ is constant in a neighborhood of the boundary, there exists a $\delta$ such that if $\|q\|_{H^{d/2+k+\varepsilon}}\leq\delta$, then the same set of illuminations work for all rotations.*

By open set, we mean $\{\tilde{g}_{i,j}\}$ sufficiently close to $\{g_{i,j}\}$ in any topology of sufficiently smooth functions such as $C^{2}(\bar{\Omega})$. The proof combines ideas from [4] and [9].

Proof.

We choose CGO-parameters $\rho_{1}=t(\mathbf{e}_{2}+i\mathbf{e}_{1})$, $\rho_{2}=t(\mathbf{e}_{3}+i\mathbf{e}_{1})$ and $\rho_{3}=t(-\mathbf{e}_{2}+i\mathbf{e}_{1})$, and choose the corresponding (real-valued) solutions to be

[TABLE]

After some algebra, it follows that

[TABLE]

where $C>0$ is independent of $t$. We take the smallest $t$ so large that $|O(t^{-1})|\leq 1/2$ in all expressions above. Since the zero sets of the determinants are disjoint, it follows that there are complementary sets $\Omega_{1},\Omega_{2}$ such that $\Omega=\Omega_{1}\cup\Omega_{2}$ and that $\{\nabla u_{i,1},\nabla u_{i,2},\nabla u_{i,3}\}$ is a basis for $\mathbb{R}^{3}$ on $\Omega_{1}$, and $\{\nabla u_{i,1},\nabla u_{i,2},\nabla u_{i,4}\}$ is a basis for $\mathbb{R}^{3}$ on $\Omega_{2}$ (see [9] for details on the construction). Correspondingly, $u_{i,5}(x)\neq 0$ for $x\in\Omega_{1}$ and $u_{i,6}(x)\neq 0$ for $x\in\Omega_{2}$. As a consequence, the vector fields

[TABLE]

form a basis for $\mathbb{R}^{3}$ for every $x\in\Omega$. Hence we can choose illuminations such that $g_{i,j}=u_{i,j}|_{\partial\Omega}$. Further, by continuity of the mapping $g\in C^{2}(\overline{\Omega})\mapsto u\in C^{1}(\overline{\Omega})$, it suffices to take $g_{i,j}$ close to $u_{i,j}|_{\partial\Omega}$, i.e., $\|g_{i,j}-u_{i,j}\|_{C^{2}(\overline{\Omega})}\leq\varepsilon$, for some $\varepsilon>0$ sufficiently small [4].

Now, let $R\in\mathbb{R}^{3\times 3}$ be a rotation matrix, and set $\tilde{q}(x)=q(Rx)$. Then equation (24) with $q$ replaced by $\tilde{q}$ corresponds to a rotation $i$ of the object or the illumination patterns. Let

[TABLE]

where $\tilde{\psi}_{\rho_{i}}$ is the perturbation term with respect to $\tilde{q}$. We assume that $\|g_{i,1}-u_{i,j}|_{\partial\Omega}\|_{C^{2}(\overline{\Omega})}\leq\varepsilon/2$. Then

[TABLE]

where the last inequality follows by Sobolev embedding and the bound in (25), setting $s>3/2+2$, and the fact that $\|q\|_{H^{s}}=\|\tilde{q}\|_{H^{s}}$. The result now follows by requiring

[TABLE]

The same argument then works for all illuminations.

∎

The drawback with using CGO-solutions is that the results are not constructive, since they rely on the unknown parameters. Numerical simulations [11] show that most illuminations will in fact work, but we know that in dimension $d\geq 3$, we cannot find illuminations that provide ellipticity conditions independently of the coefficients [3].

3.5 Stability and uniqueness

Written out, the principal part of measurement object $\delta\mathcal{V}$ in (19) is of the form

[TABLE]

with components

[TABLE]

Note that for $\delta\mathcal{V}_{1}$ and $\delta\mathcal{V}_{2}$, we have by standard mapping properties of PDOs on Sobolev spaces that

[TABLE]

for $s\in\mathbb{R}$ and constants $C_{s},\tilde{C}_{s}>0$.

We now associate to the system (18) two sets of integers $s=(s_{1},s_{2})$ and $t=(t_{1},t_{2})$ such that for each entry of the $2\times 2$ symbol $\sum_{i=1}^{M}(\kappa_{i})^{*}\kappa_{i}$ we have that $\left(\sum_{i=1}^{M}(\kappa_{i})^{*}\kappa_{i}(\kappa^{0}_{i})^{*}\kappa^{0}_{i}\right)_{(i,j)}\in S^{s_{i}-t_{j}}(\Omega,\mathbb{R}^{n})$. Since for each $i$

[TABLE]

we can choose and $s=(-1,0)$ and $t=(1,0)$ so that the above entry $(k,l)$ is of order $s_{k}-t_{l}$.

The operator $(\kappa^{0}_{i})^{*}\kappa^{0}_{i}$ is said to be elliptic of type $(s,t)$ in the Douglis-Nirenberg sense if it is invertible for every $(x,\xi)\in T^{*}\Omega$, see [19, 1]. When $(\kappa^{0}_{i})^{*}\kappa^{0}_{i}$ is elliptic, it has a parametrix $Q(x,\xi)$, an inverse modulo smooth terms, and applying $Q$ to (20) yields the estimate (since $QP$ is smooth as well)

[TABLE]

for any $s^{\prime}\in\mathbb{R}$ ([19], Lemma 1.0.2’, or [24]). The constant $C$ depends on $D,\bar{D},\sigma,\bar{\sigma},\Omega$, on $s^{\prime}$ and the illuminations. We can now complete the proof of Theorem 1.

Proof.

That boundary conditions such that $(\kappa^{0}_{i})^{*}\kappa^{0}_{i}$ is elliptic of the $(s,t)$ type exist follows from Theorem 3.4, and hence the estimate (28) holds. Let $s^{\prime}=-1$ in (28). For $u\in L^{2}(\Omega)\subset H^{-1}(\Omega)$, $\text{supp }u\subset B_{\varepsilon}$, we have

[TABLE]

where we use the Poincaré inequality to get the $\varepsilon$-dependence; see [16]. A priori, $u\in H_{0}^{1}(B_{\varepsilon})$, and we have that

[TABLE]

Using the estimates in (27) we get that that for some constant $C_{v}>0$

[TABLE]

Last, the time-reversal operator $A$ is bounded ([25], Theorem 2.4 with $\theta=1$), and we have $\|A\delta V_{i,j}\|_{H^{1}(\Omega)}\leq C_{A}\|\delta V_{i,j}\|_{H^{1}([0,T]\times\partial\Omega_{M})}.$ Under the assumption that $\text{supp }\delta D,\text{supp }\delta\sigma\subset B_{\varepsilon}$, we make use of the inequalities (29) and rearrange (28) to get

[TABLE]

We can then choose $\varepsilon$ such that $(1-C\tilde{C}\varepsilon)>0$. Recalling Theorem 3.4, we have $N=2d$ and we get (11). ∎

4 Remarks and Conclusions

In a sufficiently idealized setting, where the ultrasound propagation may be modeled by a wave equation with reasonably well known sound speed (and nagging attenuation effects may be neglected [21]), quantitative PAT displays favorable mathematical properties, as a composition of two reasonably well-posed inverse problems. This holds in the setting where (sufficiently) full ultrasound measurements are available, which is not always practical [14, 27].

We show in this paper that PAT enjoys favorable stability properties as well when a similar amount of measurement is collected for optical illuminations and ultrasound detectors that are allowed to rotate during the acquisition procedure so that different parts of the ultrasound measurements are generated by different optical illuminations.

The absence of complete ultrasound measurement for each optical illumination renders the inversion of a single wave equation ill-posed. Only when measurements from all optical sources are accounted for can one expect stable reconstructions. In such settings, where only local pieces of information are available for each illumination, it is difficult to envision a direct global inversion procedure. Rather, it is the ideal setting to apply micro-local methods, whose roles are precisely to propagate (phase-space-) local information through wave or elliptic equations.

Writing the difference of sufficiently smooth nonlinear functionals ${\cal F}(u)-{\cal F}(v)$ as a general functional ${\cal G}(u,v,u-v)$ that is linear in its last component, a procedure that applies in a great variety of contexts, we can write the difference of measurements as a (standard) elliptic Fourier integral operator applied to $u-v$ with a symbol that depends on $(u,v)$. This imposes that $(u,v)$, for us here the sound speed and the optical coefficients as well as the solutions to the diffusion equation, be sufficiently smooth. This procedure has been applied for the second step of many hybrid inverse problems in [24]. The inversion procedure then necessarily provides conditional stability, that is to say stable reconstructions up to the possibly non-trivial kernel of a compact operator. In this paper, we chose to make a smallness assumption on the support of the coefficients of interest to show that the kernel of said operator was necessarily trivial, a trick that is certainly not new [18].

These are our two main somewhat unnecessary condition: very large smoothness of coefficients and smallness of support.

Such strong smoothness assumptions are not necessary in standard PAT [11] (although no known results hold for arbitrary coefficients; for instance $D$ and $\sigma$ arbitrary measurable functions that are bounded between $1$ and $2$, say). They can also be avoided for more general hybrid inverse problems by writing the inverse problem as a coupled system of nonlinear partial differential equations for $(D,\sigma,u_{j})$, where the PDO calculus of the second step of PAT is replaced by a better behaving potential theory [7]. In that setting, one can also sometimes apply a unique continuation principle that allows one to obtain an injectivity result independently of the size of the support of the unknown coefficients.

However, in our setting of rotating measurements, it seems unclear how the micro-local pursuit of propagation of singularities can be totally avoided. Within that context, unnecessary smoothness assumptions seem necessary (realizing the inherent contradiction). As we mentioned in the introduction, all the smoothness we need is for $s^{\prime}$ in (28) to be strictly negative. This can be achieved by Taylor expansions involving finitely many terms in the definition of PDO and FIO operators as well as in their composition. How many (a sufficiently large finite number) dictates how much smoothness our coefficients need to verify.

These technical constraints notwithstanding, we expect standard quantitative PAT and QPAT in a rotating setting to display very similar resolution capabilities.

Acknowledgment

The work of GB was partially funded by the US National Science Foundation and the US Office of Naval Research.