On multi-spectral quantitative photoacoustic tomography in diffusive regime

In photoacoustic tomography (PAT), near infrared (NIR) light propagates into a medium of interest. As a fraction of the incoming light energy is absorbed, the medium heats up. This results in mechanical expansion and the generation of compressive (acoustic) waves. The acoustic waves propagate to the boundary of the medium where they are measured. From knowledge of these (time-dependent) acoustic measurements, one attempts to reconstruct the diffusion, absorption and Grüneisen thermodynamic properties of the medium of interest [2, 7, 19, 31, 33, 36, 39, 47, 55–57].

The reconstruction problem in PAT is a two-step process. In the first step, one reconstructs the absorbed energy map (multiplied by a coefficient), H defined in (2) below, from the measured acoustic signals on the surface of the medium. This is a relatively well-known inverse source problem for the wave equation that has been extensively studied in the past [1, 16, 21, 22, 25–28, 32, 34, 37, 39, 40, 50–53]. In the second step, often called quantitative PAT (qPAT), the objective is to reconstruct the diffusion coefficient D in (1), the absorption coefficient σ in (1) and the Grüneisen coefficient Γ in (2), from the reconstructed energy data H. This step has recently attracted significant attention from both mathematical [3, 4, 8–10] and computational [11, 17, 18, 20, 24, 23, 35, 45, 46, 59, 61] perspectives. Mathematically, the reconstruction problem in this step is very similar to the diffuse optical tomography (DOT) problem except that in DOT we have boundary current data, while here in qPAT we have interior absorbed energy data. The use of interior data improves the stability, and thus the resolution, of the reconstruction significantly. This is the main advantage of PAT over DOT.

Most of the studies in the past have assumed that the Grüneisen coefficient is a constant and is known so that only the absorption and diffusion coefficients are to be reconstructed. As recently pointed out in [17, 19, 48], it is more realistic to consider Γ as a function of space. In this work, we treat Γ also as an unknown. We thus aim at both Γ and the optical coefficients. Note that unlike the reconstruction of the absorption coefficient σ and the diffusion coefficient D, the reconstruction of Γ is a linear inverse problem since the data H we have are linearly proportional to Γ. This means that if both σ and D are known, then we can simply solve the diffusion equation (1) and compute the energy E = σu. Then Γ is reconstructed as Γ = H/E.

The work in [9] has shown that with data at only one wavelength, the reconstruction of the three coefficients simultaneously is not unique. In fact, only two of the three coefficients can be reconstructed without other a priori information. The objective of this work is exactly to reconstruct all three coefficients by using measured data of different wavelengths and a priori information on the form of the coefficients. We show that the reconstruction can be done when the coefficients take specific forms. One specific case of practical importance is that if the dependence of the coefficient on the spatial variable and on the wavelength variable can be separated, then they can be reconstructed with spectral data stably.

The paper is structured as follows. In section 2, we reformulate the inverse problem of quantitative PAT in the setting of multiple wavelength illuminations and present uniqueness and stability theory regarding the inverse problem in simplified settings. We then present in section 3 numerical procedures to reconstruct the unknown coefficients. We present some numerical simulation results with synthetic data in section 4. Conclusions and further remarks are offered in section 5.

Let us denote by X a bounded domain in $\mathbb R^d$ (d = 2, 3) with the smooth boundary ∂X, and $\Lambda \subset \mathbb R_+$ the set of wavelengths at which the interior data are constructed. In the diffusive regime, the density of photons at the wavelength λ, $u(\mathbf x,\lambda )$ solves the following diffusion equation:

$$\begin{eqnarray} &&\begin{array}{l} -\nabla \cdot D(\mathbf x,\lambda )\nabla u(\mathbf x,\lambda ) + \sigma (\mathbf x,\lambda ) u = 0,\quad X\times \Lambda \\ u(\mathbf x,\lambda ) = g(\mathbf x,\lambda ),\quad \partial X\times \Lambda, \end{array} \end{eqnarray} \tag{ 1 }$$

where $D(\mathbf x,\lambda )>0$ and $\sigma (\mathbf x,\lambda )>0$ are the wavelength-dependent diffusion and absorption coefficients, respectively, and $g(\mathbf x,\lambda )$ is the illumination pattern in photoacoustic experiments. The wavelength-dependent interior data constructed from the inversion of the acoustic problem are given by

$$\begin{equation} H(\mathbf x,\lambda ) = \Gamma (\mathbf x,\lambda )\sigma (\mathbf x,\lambda )u(\mathbf x,\lambda ), \quad (\mathbf x,\lambda )\in X\times \Lambda , \end{equation} \tag{ 2 }$$

where $\Gamma (\mathbf x,\lambda )>0$ is the non-dimensional Grüneisen coefficient that relates the absorbed optical energy σu to the increase in pressure H. The thermodynamic coefficient Γ measures the photoacoustic efficiency of the medium.

The problem of multi-spectral qPAT is to reconstruct the coefficients $D(\mathbf x,\lambda )$, $\sigma (\mathbf x,\lambda )$ and $\Gamma (\mathbf x,\lambda )$ from interior data of the form (2). The results in [9] state that without further a priori information, we cannot uniquely reconstruct all three coefficients. In fact, only two quantities related to the three coefficients, i.e. the functionals ${\chi }(\mathbf x,\lambda )$ and $q(\mathbf x,\lambda ),$ defined in (5) below can be reconstructed. The objective of this work is to show that when data of multiple wavelengths are available, we can recover the uniqueness in the reconstructions with relatively little additional a priori information. To obtain results that are as general as possible, we allow the diffusion and absorption coefficients to be arbitrary functions of the wavelength and assume that the Grüneisen coefficient can be written as the product of a function of space and a function of wavelength. More precisely, the coefficients take the following form:

$$\begin{equation} D=D(\mathbf x,\lambda ), \quad \sigma =\sigma (\mathbf x,\lambda ),\quad \Gamma (\mathbf x,\lambda )=\gamma (\lambda )\Gamma (\mathbf x), \end{equation} \tag{ 3 }$$

where γ(λ) is a function assumed to be known. This assumption is not sufficient to guarantee unique reconstructions. We need to assume a little more on the coefficients. Our main assumption is as follows.

• (A1)  
  There exist two known wavelengths λ₁, λ₂ ∈ Λ so that $D(\mathbf x,\lambda _1)=\rho D(\mathbf x,\lambda _2)$ for some known positive constant ρ, while $\frac{\sigma }{\sqrt{D}}(\mathbf x,\lambda _1)\ne \frac{\sigma }{\sqrt{D}}(\mathbf x,\lambda _2)$ for all $\mathbf x\in X$.This assumption essentially requires that the dependence of the diffusion and absorption coefficients is different for at least two wavelengths. Otherwise, we would simply be able to eliminate the spectral dependence from equation (1). We will see later on that this assumption is satisfied by most coefficient models employed by researchers in the community. We also need some assumptions on the regularity and boundary values of the coefficients. We assume, denoting by $\mathcal C^p(X)$ the space of p-times differentiable functions in X, that
• (A2)  
  the function $0<D(\cdot , \lambda )\in \mathcal C^2(\bar{X})$, and its boundary value D_|∂X is known. The functions $0<\sigma (\cdot , \lambda ), 0<\Gamma (\cdot , \lambda ) \in \mathcal C^1(\bar{X})$.

We now present the main uniqueness result on reconstructions with multi-spectral data.

Theorem 2.1. Let (D, σ, Γ) and $(\tilde{D}, \tilde{\sigma }, \tilde{\Gamma })$ be two sets of coefficients given in (3), with γ(λ) known, satisfying assumptions (A1) and (A2). Then there exists an open set of illuminations $\big (g_1(\mathbf x,\lambda ),\ g_2(\mathbf x,\lambda )\big )$ such that the equality of the data

$$\begin{equation*} \lbrace H_1(\mathbf x,\lambda ),\ H_2(\mathbf x,\lambda )\rbrace = \lbrace \tilde{H}_1(\mathbf x,\lambda ),\ \tilde{H}_2(\mathbf x,\lambda )\rbrace \end{equation*}$$

in X × Λ implies that

$$\begin{equation*} \lbrace D(\mathbf x,\lambda ),\ \sigma (\mathbf x,\lambda ),\ \Gamma (\mathbf x)\rbrace = \lbrace \tilde{D}(\mathbf x,\lambda ),\ \tilde{\sigma }(\mathbf x,\lambda ),\ \tilde{\Gamma }(\mathbf x)\rbrace \end{equation*}$$

in X × Λ provided that the values of the diffusion coefficient agree on the boundary: $D_{|\partial X}=\tilde{D}_{|\partial X}$ for all λ ∈ Λ.

Proof. With the regularity assumptions on the diffusion coefficient, we can recast equation (1), using the Liouville transform $v=\sqrt{D} u$, as

$$\begin{equation} \begin{array}{@{}ll}\Delta v(\mathbf x,\lambda ) + q(\mathbf x,\lambda ) v(\mathbf x,\lambda ) = 0, & X\times \Lambda \\ v(\mathbf x,\lambda ) = \tilde{g}(\mathbf x,\lambda ):=\sqrt{D}(\mathbf x,\lambda ) g(\mathbf x,\lambda ), & \partial X\times \Lambda , \end{array} \end{equation} \tag{ 4 }$$

with the interior data H(x, λ) = v(x, λ)/χ(x, λ), where χ and q are defined respectively as

$$\begin{equation} {\chi }=\frac{\sqrt{D}}{\Gamma \sigma }, \qquad -q=\frac{\Delta \sqrt{D}}{\sqrt{D}}+\frac{\sigma }{D} . \end{equation} \tag{ 5 }$$

Let us denote by $v_1(\mathbf x,\lambda )$ and $v_2(\mathbf x,\lambda )$ the solutions of (4) with the illuminations $\tilde{g}_1(\mathbf x,\lambda )$ and $\tilde{g}_2(\mathbf x,\lambda ),$ respectively. Then it is straightforward to check that [9]

$$\begin{equation} \begin{array}{@{}ll}-\nabla \cdot v_1^2 \nabla \dsty\frac{v_2}{v_1} = 0, & X\times \Lambda \\\ms v_1^2 (\mathbf x,\lambda ) = {\tilde{g}}_1^2, & \partial X\times \Lambda . \end{array} \end{equation} \tag{ 6 }$$

Since $\frac{v_2}{v_1}=\frac{H_2}{H_1}$ is known, this is a transport equation for $v_1^2(\mathbf x,\lambda )$ with the vector field $\nabla \frac{H_2}{H_1}$. By the results of [9], there exist the well-chosen illuminations g₁ and g₂ (and thus $\tilde{g}_1$ and $\tilde{g}_2$) such that the vector field is smooth enough with the positive modulus $|\nabla \frac{H_2}{H_1}|>0$. This ensures that the transport equation is uniquely solvable for v²₁. Once v²₁ is reconstructed, we can reconstruct χ = v₁/H₁ and then q = −Δv₁/v₁.

Now let λ₁, λ₂ ∈ Λ be two wavelengths as in assumption (A1). We first reconstruct the two pairs of functionals (${\chi }(\mathbf x,\lambda _1),\ q(\mathbf x,\lambda _1)$) and (${\chi }(\mathbf x,\lambda _2),\ q(\mathbf x,\lambda _2)$). With the assumption that $D(\mathbf x,\lambda _1)=\rho D(\mathbf x,\lambda _2)$, we obtain after some algebra using the above equations that

$$\begin{equation} \begin{array}{@{}ll}\Delta \sqrt{D}(\mathbf x,\lambda _1) + Q(\mathbf x)\sqrt{D}(\mathbf x,\lambda _1) =0, & X\\\ms \sqrt{D}(\mathbf x,\lambda _1) = \sqrt{D}(\mathbf x,\lambda _1)_{|\partial X}, & \partial X, \end{array} \end{equation} \tag{ 7 }$$

where the function $Q(\mathbf x)$ is defined as

$$\begin{equation} Q(\mathbf x) = \frac{\gamma (\lambda _1){\chi }(\mathbf x,\lambda _1)q(\mathbf x,\lambda _1)-\sqrt{\rho }\gamma (\lambda _2){\chi }(\mathbf x,\lambda _2)q(\mathbf x,\lambda _2)}{\gamma (\lambda _1){\chi }(\mathbf x,\lambda _1)-\sqrt{\rho }\gamma (\lambda _2){\chi }(\mathbf x,\lambda _2)}. \end{equation} \tag{ 8 }$$

Using the conditions in assumptions (A1) and (A2), the denominator in Q does not vanish and Q is bounded. The elliptic equation (7) can then be solved uniquely as shown in [9] to reconstruct $\sqrt{D}(\mathbf x,\lambda _1)$. We then find

$$\begin{eqnarray} \sigma (\mathbf x,\lambda _1)&=&-(Dq)(\mathbf x,\lambda _1) - (\sqrt{D} \Delta \sqrt{D}) (\mathbf x,\lambda _1), \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \Gamma (\mathbf x) &=& \frac{\sqrt{D}(\mathbf x,\lambda _1)}{{\chi }(\mathbf x,\lambda _1)\sigma (\mathbf x,\lambda _1)} . \end{eqnarray} \tag{ 10 }$$

Once $\Gamma (\mathbf x)$ is reconstructed, then so is $\Gamma (\mathbf x,\lambda )=\Gamma (\mathbf x)\gamma (\lambda )$. The results in [9, corollary 2.3] then allow us to reconstruct $D(\mathbf x,\lambda )$ and $\sigma (\mathbf x,\lambda )$ for all λ ∈ Λ.

It has been shown in [9] that the reconstruction of v₁ (and thus that of χ and q as well) is Lipschitz stable for each fixed wavelength. Because the solution of (7) is stable with respect to Q, we deduce that the reconstruction of the $D(\mathbf x,\lambda _1)$ (and thus $\Gamma (\mathbf x)$) is stable. Precise stability estimates similar to those in [9, theorem 2.4] can be derived.

The coefficient model (3) is quite general and covers many of the models used in the literature. For instance, we may consider the following standard model [13–15, 17, 19, 29, 35, 41, 42, 48, 49, 54, 58, 60]:

$$\begin{equation} \begin{array}{@{}l}D(\mathbf x,\lambda )=\alpha (\lambda ) D(\mathbf x), \quad \Gamma (\mathbf x,\lambda )=\gamma (\lambda )\Gamma (\mathbf x)\\ \sigma (\mathbf x,\lambda )=\displaystyle \sum _{k=1}^K\beta _k(\lambda )\sigma _k(\mathbf x), \end{array} \end{equation} \tag{ 11 }$$

where the functions α(λ), {β_k(λ)}^K_{k = 1} and γ(λ) are assumed to be known. In other words, all three-coefficient functions can be written as products of functions of different variables. Moreover, the absorption coefficient contains multiple components. This is the parameter model proposed in [17, 19, 29, 35, 42, 49] to reconstruct chromophore/fluorochrome concentrations from optical measurements. The following result regarding model (11) is a natural extension of theorem 2.1.

Corollary 2.2. Let $D(\mathbf x,\lambda )$, $\sigma (\mathbf x,\lambda )$ and $\Gamma (\mathbf x,\lambda )$ be as in (11) and satisfying assumptions (A1) and (A2). Assume further that there exist K wavelengths such that the matrix B with elements B_kj = β_k(λ_j) (1 ⩽ k, j ⩽ K) is non-singular. Then the coefficients

$$\begin{equation*} \lbrace D(\mathbf x),\ \lbrace \sigma _k(\mathbf x)\rbrace _{k=1}^K,\ \Gamma (\mathbf x)\rbrace \end{equation*}$$

can be uniquely reconstructed from the well-chosen data $(H_1(\mathbf x,\lambda ),H_2(\mathbf x,\lambda ))$, $(\mathbf x,\lambda )\in X\times \Lambda$.

Proof. Theorem 2.1 allows the unique reconstruction of $D(\mathbf x)$, $\Gamma (\mathbf x)$ and $\lbrace \sigma (\mathbf x,\lambda _k)\rbrace _{k=1}^K$. Since B is non-singular, we can invert the linear system $\sigma (\mathbf x,\lambda _j)= \sum _{k=1}^K\beta _k(\lambda _j)\sigma _k(\mathbf x)$, j = 1, ..., K, to recover $\lbrace \sigma _k(\mathbf x)\rbrace _{k=1}^K$.

We now present two numerical reconstruction methods for the multi-spectral qPAT problem.

3.1. The vector field method

3.2. The nonlinear least-squares method

The second approach to solve numerically the inverse problem is to reformulate the reconstruction problem as a data-driven minimization problem. More precisely, we look for the coefficients in the form of (3) that minimize the following mismatch functional:

$$\begin{equation} \Phi (D,\sigma ,\Gamma )= \frac{1}{2}\displaystyle \sum _{i=1}^N \int _{X}\int _\Lambda \big (H_i(\mathbf x,\lambda )-H_i^*(\mathbf x,\lambda )\big )^2\, {\rm d}\mathbf x\,{\rm d}\lambda , \end{equation} \tag{ 14 }$$

where N ⩾ 2 is the number of different spatial source patterns used for a fixed wavelength. In the theoretical results we presented in section 2, N = 2. The data predicted by the model are denoted as $H_i(\mathbf x,\lambda )=\Gamma (\mathbf x,\lambda )\sigma (\mathbf x,\lambda )u_i(\mathbf x,\lambda )$, with $u_i(\mathbf x,\lambda )$ being the solution of the diffusion equation (1) with the ith source pattern $g_i(\mathbf x,\lambda )$. The real interior data, that is, the data constructed from acoustic measurement, are denoted by $H_i^*(\mathbf x,\lambda )$.

We solve the minimization problem by a quasi-Newton method with the BFGS updating rule on the Hessian operator [12, 30, 38, 44]. This method requires the Fréchet derivatives of the functional Φ with respect to the unknowns D, σ and Γ. These derivatives can be computed in the following way. Let us denote by $w_i(\mathbf x,\lambda )$ the solution of the adjoint problem

$$\begin{equation} \begin{array}{l}-\nabla \cdot D(\mathbf x,\lambda )\nabla w_i(\mathbf x,\lambda ) + \sigma (\mathbf x,\lambda ) w_i = \Gamma \sigma Z_i,\quad X\times \Lambda \\ w_i(\mathbf x,\lambda ) = 0,\quad \partial X\times \Lambda , \end{array} \end{equation} \tag{ 15 }$$

with $Z_{i}\equiv H_i(\mathbf x,\lambda )-H_i^*(\mathbf x,\lambda )$. Then it is standard to show that

Theorem 3.1. The nonlinear least-squares functional $\Phi : [L^2(X\times \Lambda )]^3 \mapsto \mathbb R$ is Fréchet differentiable with respect to D, σ and Γ. The Fréchet derivatives are given respectively by

$$\begin{eqnarray} &&\left\langle \frac{\partial \Phi }{\partial D}, \hat{D} \right\rangle = \displaystyle \sum _{i=1}^N\langle \nabla u_i \cdot \nabla w_i, \hat{D}\rangle _{L^2(X\times \Lambda )}, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&\left\langle \frac{\partial \Phi }{\partial \Gamma }, \hat{\Gamma }\right\rangle = \displaystyle \sum _{i=1}^N \langle Z_i \sigma u_i,\hat{\Gamma }\rangle _{L^2(X\times \Lambda )}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} &&\left\langle \frac{\partial \Phi }{\partial \sigma }, \hat{\sigma }\right\rangle = \displaystyle \sum _{i=1}^N \langle \Gamma Z_iu_i-w_iu_i,\hat{\sigma }\rangle _{L^2(X\times \Lambda ),} \end{eqnarray} \tag{ 18 }$$

where $\langle \cdot \rangle _{L^2(X\times \Lambda )}$ denotes the usual inner product in L²(X × Λ).

When model (11) is employed, the derivatives with respect to σ have to be replaced by the derivatives with respect to each σ_k (1 ⩽ k ⩽ K). These derivatives can be simply computed by using (18) and the chain rule. More precisely,

$$\begin{equation} \left\langle \frac{\partial \Phi }{\partial \sigma _k}, \hat{\sigma }_k \right\rangle = \displaystyle \sum _{i=1}^N \langle \Gamma Z_iu_i-w_iu_i,\beta _k(\lambda )\hat{\sigma }_k\rangle _{L^2(X\times \Lambda )}, \quad 1\le k\le K. \end{equation} \tag{ 19 }$$

Compared with the vector field method in the previous section, the nonlinear least-squares method is a very general method. However, it is computationally more expensive. In each quasi-Newton iteration, we need to solve several forward and adjoint problems to evaluate the value of the objective function and the Fréchet derivatives. The method is thus slow. However, since the diffusion equation (1) and its adjoint (15) are elliptic problems that can be solved efficiently with finite-element methods, the overall reconstruction cost is reasonably low.

We now present some numerical simulations using synthetic data to support the theory that is developed in the previous section. To simplify the presentation, we consider only two-dimensional problems, even though the uniqueness results and the reconstruction methods hold in three-dimensional domains as well. The domain we consider is the square: X = (0, 2 cm) × (0, 2 cm). Material properties of the medium will be provided in specific cases.

In the first set of numerical experiments, we show reconstruction results for smooth coefficients. The coefficients take the forms given in (11) with two components in the absorption coefficient, i.e. K = 2. The spectral components of the coefficients are given as follows:

$$\begin{equation} \alpha (\lambda )=\left(\frac{\lambda }{\lambda _0}\right)^{3/2}, \quad \gamma (\lambda ) = 1,\quad \beta _1(\lambda )=\frac{\lambda }{\lambda _0},\quad \beta _2(\lambda )=\frac{\lambda _0}{\lambda }, \end{equation} \tag{ 20 }$$

where the normalization wavelength λ₀ is included to control the amplitude of coefficients. These weight functions are by no means what they should exactly be in practical applications. However, the specific forms do not have impacts on the results of the reconstruction. The spatial components of the coefficients are given as

$$\begin{equation} \fl \begin{array}{l}\Gamma (\mathbf x)=0.8+0.4\tanh (4x-4), \qquad D(\mathbf x)=0.03+0.02\sin (\pi x)\sin (\pi y)\\ \sigma _1(\mathbf x)=0.2+0.1 {\rm e}^{-(x-1)^2-(y-1)^2}, \qquad \sigma _2(\mathbf x)=0.2-0.1 {\rm e}^{-(x-1)^2-(y-1)^2}. \end{array} \end{equation} \tag{ 21 }$$

Four illuminations are used, two at each wavelength. The reconstruction results with clean data, i.e. synthetic data without being polluted by additional random noise, are shown in the second row of figure 1. Note that the data used in this case indeed contain noise that comes from the discretization of the continuous PDE model, which is very small since we used a very fine discretization. The reconstructions are almost perfect, with overall quality very similar to those presented in mono-spectral data [9], as can be seen from the plots and the cross-sections (red dashed) in figure 2. The third row of figure 1 presents one typical realization of reconstruction results with data containing 5% extra additive random noise. The cross-sections of the reconstruction along the line y = 1 are shown in figure 2 (blue dashed). The results are still very accurate, with overall quality comparable to those results obtained in [9].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the second set of numerical experiments, we show reconstruction results for piecewise constant coefficients with the numerical minimization algorithm that we introduced in the previous section. The true coefficients are shown in the first row of figure 3. We again use four illuminations of two different wavelengths. The reconstructions with clean and noisy data are presented on the second and third rows of figure 3, respectively. In the case of noisy data, we regularized the nonlinear least-squares problem with the Tikhonov regularization term $\frac{\rho }{2}(\Vert \nabla D\Vert _{L^2}^2+\Vert \nabla \sigma \Vert _{L^2}^2+ \Vert \nabla \Gamma \Vert _{L^2}^2)$. The strength ρ of the regularization is taken after a few tests with different values. The effect of the regularization is mainly to smooth out the additive noise in the data. It also smooths out part of the jumps across the interfaces of coefficient discontinuities, as can be seen from the plots in the third row of figure 3. To better visualize the quality of the reconstructions, we again plot the cross-sections of the reconstructed coefficients in figure 4. We observe that discontinuous coefficients can be reconstructed as accurately as smooth coefficients, even though the discontinuity is smeared out in the former case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We presented in this paper some theoretical results for quantitative photoacoustic tomography with multi-spectral interior data to reconstruct the Grüneisen, absorption and diffusion coefficients simultaneously. We showed the uniqueness of the solution to the inverse problem under minor assumptions on the form of the coefficients. One specific case of practical importance is that, if the dependence of the coefficient on the spatial variable and on the wavelength variable can be separated, the spatial components of all three coefficients can be reconstructed with spectral data uniquely and stably.

We presented a vector field-based algorithm and a general nonlinear least-squares-based method for the numerical reconstruction of smooth and rough coefficients. We showed by numerical simulations that the reconstructions are very accurate for both types of coefficients, assuming that the interior data H constructed from measured acoustic data are accurate enough. The reconstructions are stable with Lipschitz type stability estimates that are identical to those in [9].

The uniqueness results that we obtained in this paper depend on the illuminations that are selected. In other words, uniqueness only holds when appropriate illuminations are used. Fortunately, the characterization of these 'well-selected' illuminations is identical to that presented in [9]. In fact, as pointed out in [9], many illumination pairs that we have used work, in the sense that they generate the data H₁ and H₂ such that the condition $|\nabla \frac{H_2}{H_1}|>0$ holds, and thus the transport equation (6) is uniquely and stably solvable.

The results of this paper and [9] are all based on diffusion theory for light propagation in tissues. It is well known that the phase-space radiative transport model [5, 6, 43] is a more accurate light propagation model for this application. Uniqueness and stability results for the transport model have been developed in [8] in the situation where measurements from all possible illuminations are available. Numerical simulations have also been reported [20]. It would be interesting to extend the study in [8] to include an unknown variable Grüneisen coefficient.

We would like to thank Professor Simon Arridge (University College London) for useful discussion on the simplified coefficient model (11) in optical and photoacoustic tomography. KR was supported in part by NSF grant DMS-0914825. GB was supported in part by NSF grants DMS-0554097 and DMS-0804696.