Implicit learning to determine variable sound speed and the reconstruction operator in photoacoustic tomography

Photoacoustic tomography (PAT) is an imaging method that uses non-ionized laser pulses and ultrasound to produce detailed images of the internal structure of biological tissue. The method is non-destructive, economical, and less harmful than other imaging options because it uses non-ionizing radiation [31]. For these reasons, it is used in a variety of biomedical applications, including skin melanoma detection, breast cancer detection, blood oxygenation mapping, tumor angiogenesis monitoring, functional brain imaging, and methemoglobin measurement [36].

In PAT, when the target object is irradiated with a non-ionizing laser pulse, the absorbed pulse generates a photoacoustic effect in which rapid thermal expansion leads to the production of acoustic waves, a phenomenon that was first discovered by Bell [7]. The resulting acoustic waves are measured using ultrasound detectors as the data and used to reconstruct images of the target object. This method is particularly useful for visualizing structures in optically opaque materials such as biological tissue, because it combines the high contrast of optical imaging with the high spatial resolution of ultrasound imaging (for more details, see [17, 21, 34]).

One of the goals of PAT is to reconstruct initial pressure $f$, which potentially contains important biological information such as the presence and location of cancer cells, from the acoustic waves measured by the detector. These acoustic waves, denoted by $p$, follow the wave equation:

$$\left\{\begin{array}[]{ll}\partial_{t}^{2}p({\mathbf{x}},t)=c({\mathbf{x}})\triangle_{\mathbf{x}}p({\mathbf{x}},t)&({\mathbf{x}},t)\in\mathbb{R}^{2}\times[0,\infty),\\ p({\mathbf{x}},t)|_{t=0}=f({\mathbf{x}})\text{ and }\partial_{t}p({\mathbf{x}},t)|_{t=0}=0&{\mathbf{x}}\in\mathbb{R}^{2}\end{array}\right.$$

(1)

Here, $c(\mathbf{x})$ represents the sound speed at location ${\mathbf{x}}$, which is assumed to be continuous and bounded by two constants $c_{M}$ and $c_{m}$ such that $c_{M}\geq c(\mathbf{x})\geq c_{m}\geq 0$. The measurements are made on the boundary of a region of interest $\Omega\subset\mathbb{R}^{2}$ under the reasonable assumption that $f$ has compact support within the bounded domain $\Omega$. We can then define wave forward operator $\mathcal{W}_{c}$, which maps initial pressure $f$ to solution $p$, i.e., $\mathcal{W}_{c}f=p$. The challenge then focuses on precisely reconstructing $f$ from the boundary measurements to derive an accurate internal image of the target object.

Most studies on the reconstruction of initial pressure $f$ from $\mathcal{W}_{c}{f}|_{\partial\Omega\times[0,\infty]}$ assume a known sound speed $c$. The use of a constant $c$ has been studied analytically in [12, 35, 37] (for further mathematical details, see [2, 20, 21] and references therein), while other analytical studies have employed a variable $c$. Agranovsky and Kuchment study reconstruction of the initial pressure from Dirichlet data with a known variable sound speed in three-dimensional space [1]. Moon et al. also study the singular value decomposition of the wave forward operator with a known radial sound speed in $n$-dimensional space [23]. Stefanov and Uhlmann address a more general problem with a Riemann metric instead of the sound speed [30]. Other studies detailing numerical approaches include [6, 15, 24]. Research has also been conducted on recovering $c$ when initial pressure $f$ is known [29], while the sufficient conditions for simultaneously reconstructing $f$ and $c$ are discussed in [22].

Recently, deep learning methods have been applied to PAT [32, 14] image reconstruction [3], handling limited data setups [4, 5, 13, 18, 26], and achieving a super-resolution [5]. Many of these studies are based on supervised learning, for which initial pressure $f$ is the target data. However, in practice, it cannot necessarily be assumed that the target data is known in PAT because the initial pressure represents the internal structure of the object. Therefore, methods for training the reconstruction operator without using the target data need to be considered. In line with this, in this paper, we propose an unsupervised learning method to implicitly estimate the variable speed of sound $c$ and the reconstruction operator. The proposed method only utilizes paired Dirichlet and Neumann data boundary values and contributes to PAT research by demonstrating the use of implicit learning techniques.

The rest of this paper is structured as follows. The next section presents the formulation of our problem. In Section 3, we describe our proposed framework and its loss function. The numerical results are presented in Section 4. Finally, Section 5 summarizes our research and its contributions to PAT, focusing on the variable speed of sound and initial pressure.

Building on (1), we consider the problem of estimating variable sound speed $c$ and the reconstruction operator that recovers initial pressure $f$ from Dirichlet and Neumann boundary data. We assume that $c\in C^{\infty}(\mathbb{R}^{2})$ satisfies $c_{M}\geq c(\mathbf{x})\geq c_{m}\geq 0$. Theorem 3.3 in [29] is a uniqueness theorem associated with this problem. According to this theorem, for a given paired dataset of initial values and Dirichlet data, sound speed $c$ and the reconstruction operator can be uniquely determined (see the Appendix for further details of the problem formulation, methods, and experimental results for this type of paired dataset is given). However, because it is an unreasonable assumption that the ground truth for the initial data can be used in PAT, obtaining paired data is impractical. Therefore, in this paper, we address this problem without using initial data.

For convenience, Dirichlet data $\mathcal{W}_{c}f|_{\partial B\times[0,\infty)}$ is denoted as $\mathcal{D}_{c}f$ and Neumann data $\partial_{\nu}\mathcal{W}_{c}f|_{\partial B\times[0,\infty)}$ is denoted as $\mathcal{N}_{c}f$, where $\nu$ is the unit outward normal vector at the detector surface $\partial B$ and $B$ is the unit ball in $\mathbb{R}^{2}$.

We can now discuss the uniqueness of $c$ in Problem 1. For this, we recall the Calderón problem on the conductivity equation [11, 33]: For the equation

$$\nabla\cdot(\gamma({\mathbf{x}})\nabla u({\mathbf{x}}))=0\qquad{\mathbf{x}}\in B,$$

(3)

with Dirichlet condition $u|_{\partial B}=g$, the problem is uniquely determining conductivity function $\gamma$ from the knowledge of the (bounded linear) Dirichlet to Neumann map $\Lambda_{\gamma}:H^{1/2}(\partial B)\to H^{-1/2}(\partial B)$ defined by

$$\Lambda_{\gamma}(g)=\partial_{\nu}u|_{\partial B},$$

where $H^{\alpha}(B),\alpha\in\mathbb{R}$ is the Sobolev space.

Note that conductivity equation (3) is equivalent to

$$(-\Delta+q({\mathbf{x}}))w({\mathbf{x}})=0\qquad{\mathbf{x}}\in B,$$

(4)

where $q({\mathbf{x}})=\dfrac{\nabla\gamma^{1/2}({\mathbf{x}})}{\gamma^{1/2}({\mathbf{x}})}$ and $w({\mathbf{x}})=\gamma^{1/2}({\mathbf{x}})u({\mathbf{x}})$.

Based on wave equation (1), we propose the following conjecture.

Conjecture: Assume that $c_{1}$ and $c_{2}$ are the known constant $c_{0}$ on $B^{c}$. If, for any $f\in L^{2}(B)$,

$$\mathcal{D}_{c_{1}}f=\mathcal{D}_{c_{2}}f\text{ and }\mathcal{N}_{c_{1}}f=\mathcal{N}_{c_{2}}f,$$

then $c_{1}=c_{2}$ on $\mathbb{R}^{2}$.

Because $c_{1}$ and $c_{2}$ are known constant $c_{0}$ on $B^{c}$, it suffices to show that $c_{1}=c_{2}$ on $B$. Let us define

$$\overline{p}({\mathbf{x}},t):=\left\{\begin{array}[]{ll}p({\mathbf{x}},t)|_{\overline{B}\times[0,\infty)}&t\geq 0,\\ p({\mathbf{x}},-t)|_{\overline{B}\times[0,\infty)}&t<0,\end{array}\right.$$

where $p$ is the solution for (1). Then, $\overline{p}$ satisfies

$$\left\{\begin{array}[]{ll}\partial_{t}^{2}\overline{p}({\mathbf{x}},t)=c({\mathbf{x}})\Delta_{\mathbf{x}}\overline{p}({\mathbf{x}},t)&({\mathbf{x}},t)\in B\times\mathbb{R},\\ \overline{p}({\mathbf{x}},0)=f({\mathbf{x}})\text{ and }\partial_{t}\overline{p}({\mathbf{x}},t)|_{t=0}=0&{\mathbf{x}}\in B.\end{array}\right.$$

(5)

Taking the Fourier transform of (5) with respect to the time variable $t$, we have

$$\left(-\Delta_{\mathbf{x}}+\frac{(2\pi\text{i}\omega)^{2}}{c({\mathbf{x}})}\right)\mathcal{F}_{t}(\overline{p}({\mathbf{x}},\cdot))(\omega)=0.$$

The Calderón problem appears to give us the uniqueness of $c$ in (5) given (2). However, there are gaps in this application. For $f\in L^{2}(B)$ we cannot be sure that $\mathcal{D}_{c}f\in H^{1/2}(\partial B)$. Moreover, it is not easy for set $\{\mathcal{D}_{c}f:f\in L^{2}(B)\}$ to be equal to $H^{1/2}(\partial B)$. Despite these gaps, our conjecture is validated experimentally in Section 4.

In practical terms, we focus on the following problem, emphasizing that only a finite collection of data can be used in real-world applications.

Our problem has the following difficulties:

1. 1.

   Target data $f$ in PAT is unavailable, and

2. 2.

   The explicit formula of the wave forward operator is unknown.

In this paper, we propose an implicit learning method for estimating $c$ and $\mathcal{D}_{c}^{-1}$. The proposed method uses a paired dataset of Dirichlet and Neumann data and an iterative method for the wave forward operator.

Our goal is to estimate $c$ and $\mathcal{D}_{c}^{-1}$ simultaneously from given data set $\mathscr{T}=\{(\mathcal{D}_{c}f_{i},\mathcal{N}_{c}f_{i})\}_{i=1}^{N}$. The proposed framework is shown in Figure 1 and consists of three components:

1. 1.

   Sound speed network $\tilde{c}$

2. 2.

   Reconstruction network $\mathcal{R}$

3. 3.

   Wave forward operator $\mathcal{W}_{\tilde{c}}$

Sound speed network $\tilde{c}$ is designed to estimate the unknown $c$ in (1). Simultaneously, reconstruction network $\mathcal{R}$ approximates ${\mathcal{D}_{c}^{-1}}$. This network takes Dirichlet data $\mathcal{D}_{c}f$ and outputs a reconstruction of initial data $\tilde{f}$. Wave forward operator $\mathcal{W}_{\tilde{c}}$, using an iterative scheme, solves wave equation (1) with the estimated $\tilde{c}$ and the reconstructed $\tilde{f}$. From solution $\tilde{p}$, Dirichlet $\mathcal{D}_{\tilde{c}}\tilde{f}$ and Neumann $\mathcal{N}_{\tilde{c}}\tilde{f}$ data are obtained.

If the proposed conjecture is valid, then

$$||\mathcal{D}_{\tilde{c}}\tilde{f}-\mathcal{D}_{c}f||_{2}=0\quad\text{and}\quad||\mathcal{N}_{\tilde{c}}\tilde{f}-\mathcal{N}_{c}f||_{2}=0\text{ for all }f\in L^{2}(B)$$

implying that $\tilde{c}$ and $c$ are the same. Theorem A guarantees that for $\tilde{c}=c$, if $||\mathcal{D}_{\tilde{c}}\tilde{f}-\mathcal{D}_{c}f||_{2}=0$ holds, then $\tilde{f}=f$. Thus, for given data set $\mathscr{T}=\{(\mathcal{D}_{c}f_{i},\mathcal{N}_{c}f_{i})\}_{i=1}^{N}$, we define the following loss function:

$$\mathcal{L}_{A}(\mathscr{T})=\sum\Big{[}\lambda_{\mathcal{D}}||\mathcal{D}_{\tilde{c}}\tilde{f}_{i}-\mathcal{D}_{c}f_{i}||_{2}^{2}+\lambda_{\mathcal{N}}||\mathcal{N}_{\tilde{c}}\tilde{f}_{i}-\mathcal{N}_{c}f_{i}||_{2}^{2}\Big{]}.$$

(7)

To effectively handle noisy data, we also incorporate a Total Variation (TV) regularization term for the reconstructed image [25]:

$$\mathcal{L}_{B}(\mathscr{T})=\lambda_{\text{TV}}\sum||\nabla\tilde{f}_{i}||_{1}\qquad\text{where }\tilde{f}_{i}=\mathcal{R}(\mathcal{D}_{c}f_{i}).$$

The overall loss function is then given by:

$$\mathcal{L}(\mathscr{T})=\mathcal{L}_{A}(\mathscr{T})+\mathcal{L}_{B}(\mathscr{T}).$$

In the following subsection, we describe each component of the architecture in detail.

We use Shepp-Logan phantoms for our simulation data. A Shepp-Logan phantom, introduced by Shepp and Logan in 1974 [27], is an artificial image representing a cross-section of the brain. It consists of several ellipses each defined by six parameters: the center coordinates of the ellipse, the major axis, the minor axis, the rotation angle, and the intensity. We create a set of phantoms $\left\{F_{i}\right\}_{i=1}^{4,096}$ by changing these parameters and obtain the Dirichlet data set

$$\left\{D_{c}^{r}F_{i}:\text{Dirichlet data for initial $F_{i}$ and sound speed $c$ on a ball of radius $r$ }\right\}_{i=1}^{4,096}$$

for a given sound speed $c$ by applying the iterative wave propagation method (8) to the phantoms. Noisy data are generated by adding Gaussian noise to 1% of the maximum value of the original data. We use 2,048 phantoms as a training set, 1,024 as a validation set, and the remaining 1,024 as a test set.

We set $c(x,y;h)=\exp(-h(x^{2}+y^{2}))$. Simulations are conducted using two different sound speed profiles:

	Type 1:	$\displaystyle c_{1}(x,y)=\frac{3}{4}\Big{(}c(x-0.5,y-0.6;2)+c(x-0.5,y+0.5;2)+c(x+0.6,y-0.45;2)$		(9)
		$\displaystyle\qquad\qquad\quad+c(x+0.55,y+0.6;2)\Big{)},$
	Type 2:	$\displaystyle c_{2}(x,y)=\frac{23}{35}\Big{(}c(x,y-0.675;2.5)+c(x-0.725,y-0.175;2.5)+c(x-0.4,y+0.625;2.5)$
		$\displaystyle\qquad\qquad\quad\quad+c(x+0.625,y-0.225;2.5)+c(x+0.45,y+0.575;2.5)\Big{)},$

which are characterized by a constant value outside of $B$.

In practice, for small a $h>0$, the Neumann data at the boundary of $B$ can be approximated by

$$\partial_{\boldsymbol{\nu}_{{\mathbf{x}}}}\mathcal{W}_{c}f({\mathbf{x}},t)\approx\frac{\mathcal{W}_{c}f({\mathbf{x}},t)-\mathcal{W}_{c}f({\mathbf{x}}-h\boldsymbol{\nu}_{{\mathbf{x}}},t)}{h},\quad{\mathbf{x}}\in\partial B$$

(10)

where $\boldsymbol{\nu}_{{\mathbf{x}}}$ is the unit outward normal vector at ${\mathbf{x}}\in\partial B$. Having Dirichlet and Neumann data pairs is equivalent to having Dirichlet data pairs ${(D_{c}^{1}F_{i},D_{c}^{1-h}F_{i})}$. Therefore, we modify the loss function described in (7) as follows:

$$\mathcal{L}=\lambda_{D}\sum\limits_{i}\left\|D_{c}^{1}F_{i}-D_{\tilde{c}}^{1}\tilde{F}_{i}\right\|_{2}^{2}+\lambda_{N}\sum\limits_{i}\left\|D_{c}^{1-h}F_{i}-D_{\tilde{c}}^{1-h}\tilde{F}_{i}\right\|_{2}^{2}+\lambda_{\text{TV}}\sum\limits_{i}||\tilde{F}_{i}||_{\text{TV}},$$

(11)

where $\tilde{F}=\mathcal{R}(D_{c}^{1}F)$.

We choose $h=0.05$, $\lambda_{D}=1$, $\lambda_{N}=1$, and $\lambda_{\text{TV}}=\sigma/5$, where $\sigma$ is the noise level. To minimize the loss function (11), we use the Adam optimizer [19] with a learning rate of $0.001$ and a batch size of $2$. The training takes $10^{5}$ iterations.

Table 1 shows the relative errors of the reconstructed image and the approximated sound speed. The results show that the proposed framework effectively reconstructs both the initial pressure and sound speed with high accuracy, maintaining its performance even with noisy data.

Figures 3 and 4 show the ground truth of the speed of sound along with their approximation output from the sound speed network. The relative error for $f$ with respect to ground truth $f_{\text{gt}}$ is defined by $\|f\|_{\text{rel}}=\dfrac{\|f_{\text{gt}}-f\|_{2}}{\|f_{\text{gt}}\|_{2}}$. The results demonstrate the ability of the network to accurately approximate the sound speed. This successful approximation over different sound speed scenarios ($c_{1},c_{2}$) illustrates the adaptability and robustness of the proposed network under varying conditions.

Figure 5 illustrates the ground truth and reconstruction results for the reconstruction network, for both noise-free and noisy data sets.

In this study, we introduce a framework for addressing the problem in PAT associated with the simultaneous estimation of sound speed $c$ and reconstruction operator $\mathcal{D}_{c}^{-1}$. Our approach employs implicit learning to approximate both $c$ and $\mathcal{D}_{c}^{-1}$ using only Dirichlet and Neumann boundary data. This greatly reduces the dependence on large labeled datasets that is typical of medical imaging.

In this Appendix, we address the problem of estimating the sound speed $c\in C^{\infty}(\mathbb{R}^{2})$ and the reconstruction operator $\mathcal{D}_{c}^{-1}$ when a paired dataset of initial pressure $f$ and Dirichlet data $\mathcal{D}_{c}f$ is given. With access to a sufficiently large dataset $\left\{(f_{i},\mathcal{D}_{c}f_{i})\right\}_{i=1}^{N}$, we can estimate the reconstruction operator using supervised learning. This involves minimizing the following loss function:

$$\frac{1}{N}\sum_{i}\|f_{i}-\mathcal{R}(\mathcal{D}_{c}f_{i})\|_{2}^{2}.$$

Here, $\mathcal{R}:\mathcal{D}_{c}f\mapsto f$ represents a neural network designed to approximate the reconstruction operator $\mathcal{D}_{c}^{-1}$.

From [29, Theorem 3.3], we have:

$$\text{If }\mathcal{D}_{\tilde{c}}f=\mathcal{D}_{c}f\text{ for all }f\in L^{2}(B),\text{ then }\tilde{c}=c\text{ on }B.$$

Based on this discussion, we define the loss function as follows:

$$\mathcal{L}=\lambda_{\mathcal{D}}||\mathcal{D}_{\tilde{c}}\tilde{f}-\mathcal{D}_{c}f||_{2}^{2}+\lambda_{I}||\tilde{f}-f||_{2}^{2}+\lambda_{\text{TV}}||\tilde{f}||_{\text{TV}},$$

(12)

where $\tilde{f}=\mathcal{R}(\mathcal{D}_{c}f)$.

Experiments are conducted using the same dataset and framework as described in the main text. For detailed results, please refer to Table 2, Figure 6, and Figure 8. The numerical results in the main text and Appendix show that the proposed implicit learning method, which does not use explicit target data, performs similarly to supervised learning. Additionally, the results in the Appendix are derived by simply modifying the loss function within the framework proposed in the main text, further confirming the applicability and robustness of our approach.

This work was supported by the National Research Foundation of Korea grant funded by the Korea government(MSIT) (NRF-2022R1C1C1003464, RS-2023-00217116 and RS-2024-00333393).