# An optimal control approach to photoacoustic tomography

Maîtine Bergounioux, Thomas Haberkorn, Yannick Privat

► To cite this version:

Maîtine Bergounioux, Thomas Haberkorn, Yannick Privat. An optimal control approach to photoacoustic tomography. 55th IEEE Conference on Decision and Control, Dec 2016, Las Vegas, France. pp.1620 - 1624, 10.1109/CDC.2016.7798497. hal-01582986

HAL Id: hal-01582986

https://hal.science/hal-01582986

Submitted on 6 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An optimal control approach to Photoacoustic Tomography

Maïtine Bergounioux, Thomas Haberkorn and Yannick Privat

**Abstract**—Quantitative photoacoustic tomography is a hybrid imaging technique for soft tissues. It consists in exciting the body to reconstruct with a laser pulse and measuring the induced acoustic waves due to the inhomogenous heating and then expansion of the tissues. We present a simplified mathematical model of this phenomenon, which writes as a system of two coupled equations, namely a wave equation and a diffusion approximation of a radiative transfer equation. Following [5], the inverse problem of reconstructing the optical absorption coefficient is written as an optimal control problem where the control is the parameter we seek to reconstruct. Necessary conditions for optimality of a control and numerical results for this approach are given in the case of a small number of sensors.

## I. INTRODUCTION

Photoacoustic imaging is a biomedical imaging technology that combines the optical absorption contrast of different media with the high resolution of acoustic waves. It is aimed at imaging soft tissues and, since it is non-ionizing and non-invasive, it has great potential as a safe way to localize heterogeneities in a body.

The general idea of the photoacoustic effect is simple. The tissue to be investigated is irradiated with a short-pulse laser at a given wavelength, usually near infrared where the optical absorption contrast is the highest. This energy is then converted to heat which expands inhomogeneously the medium. This expansion generates a pressure wave that is captured by transducers outside the body.

This phenomenon is used to reconstruct the heterogeneities in the body from measurements of the acoustic waves outside the body. This gives rise to an inverse problem on which various methods have been successfully applied. To cite a few, we can mention the popular filtered backprojection approach (see for example [7], [8], [12]), that assumes a closed observation surface and are primarily used to approximate the initial condition of the acoustic wave phenomenon. The time reversal method (see for example [9], [10]) also uses a closed observation surface to estimate the initial pressure distribution for the wave phenomenon. Extended literature exists on those two approaches. We will use the approach presented in [5] that, contrary to the two mentioned previously do not restrict the observation domain to be a closed surface and is therefore more flexible.

This work was supported by ANR Aventures ANR-12-BS01-0001
M. Bergounioux and T. Haberkorn are with Université d'Orléans, FDP-MAPMO, France
maïtine.bergounioux@univ-orleans.fr,
thomas.haberkorn@univ-orleans.fr

Y. Privat is with Université Pierre et Marie Curie, LJLL, Paris, France
yannick.privat@upmc.fr

This paper is structured as follows. In section II, we introduce the equations of the model, a diffusion equation for the illumination of the body and a wave equation for the propagation of the acoustic wave. In section III, the inverse problem we propose to solve is presented, namely the reconstruction of the optical absorption coefficient from acoustic wave measurements. In section IV, we present the optimal control formulation to this inverse problem. We also give the optimality system the optimal control has to satisfy. The originality of this approach is that we do not decompose the inverse problem in two separate ones but tackle it by directly using the coupled system of equations modeling the phenomenon, a similar approach can be found in [11]. Section V gives numerical results of the solving of our optimal control problem, when using only a small number of measurements. Finally we conclude and give some ideas of future works.

## II. PHOTOACOUSTIC MODEL

The photoacoustic effect consists in 3 main processes: the light transport, the thermal expansion and the pressure wave generation and propagation.

The light transport is the optical effect in which the source, a nanosecond laser pulse, propagates. As it propagates, the light is subjected to absorption and diffusion, and is governed by the Boltzmann equation [1]. However, we use a simplified version of this equation and model the optical effect with the diffusion equation

$$ \left\{ \begin{array}{l} \mu_a \Phi - \operatorname{div}(D \nabla \Phi) = S \quad , \text{ in } \Omega \\ \Phi(\partial \Omega) = 0 \end{array} \right. \qquad (1) $$

Where $\Phi(x)$ is the fluence, $\mu_a(x)$ the unknown absorption coefficient, $D(x)$ the known diffusion coefficient and $S(x)$ the known source. The domain $\Omega$ is the illuminated body.

The second effect is the thermal expansion that, due to its time-scale, is neglected in our model. Finally, the third effect, which is in fact coupled with the thermal expansion, is the generation and propagation of the acoustic wave. The pressure jump is generated by the energy deposition $H(x) = \Gamma(x)\mu_a(x)\Phi(x)$ that is used as the source of a wave equation

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial^2 p}{\partial t} - v_s^2 \Delta p = 0, & \text{in } (0,T) \times \mathcal{B}, \\ p([0,T], \partial \mathcal{B}) = 0 & \\ p(0, \mathcal{B}) = 1_{\Omega} \Gamma \mu_a \Phi & \\ \displaystyle \frac{\partial p}{\partial t} (0, \mathcal{B}) = 0 & \end{array} \right. \qquad (2) $$

Where $p(t,x)$ is the cumulative pressure, $\Gamma(x)$ is the known Grüneisen coefficient, $v_s(x)$ is the known speed of sound and
$1_Ω$ is the characteristic function of $Ω$ (= 1 on $Ω$, 0 outside). The set $\mathcal{B}$ contains $Ω$ and is large enough so that the acoustic wave initially in $Ω$ does not reach $\partial \mathcal{B}$.

### III. INVERSE PROBLEM

In this work, we want to reconstruct the absorption coefficient $\mu_a$ in $Ω$, from acoustic wave measurements outside $Ω$ (in $\mathcal{B} \setminus Ω$) and the knowledge of the light source $S$ that induced the wave. The setting we will consider in our numerical experiments is described in Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Setting of the 2 domains and distribution of $\mu_a$ and $D$.</figcaption></figure>

The domain in which $\mu_a$ is unknown is $Ω$, a 2 by 2 cm square. The domain $\mathcal{B}$ simply has to be large enough to allow the acoustic wave to travel without rebounding on $\partial \mathcal{B}$, this will be numerically taken care of with appropriate boundary conditions on $\partial \mathcal{B}$.

In $\mathcal{B} \supset Ω$, $\mu_a = \mu_a^0 = 0.1 \text{ cm}^{-1}$ except in an ellipse on the left side and a disk on the right side of $Ω$ where $\mu_a = \mu_a^0 + 10\%$, as pictured in Figure 1. The diffusion coefficient $D$ is assumed to be equal to $D^0 = 0.03 \text{ cm}^2$ except on a disk centered in $Ω$ where it is equal to $D^0 + 10\%$.

The measurements are taken on $\omega = \bigcup_{i \in \{N,E,S,W\}} \omega_i$ or $\omega = \bigcup_{i \in \{E,W\}} \omega_i$, see Figure 1. The $\omega_i$, $i \in \{N,E,S,W\}$, are small disks that in our numerical experiments will be represented by one pixel of the finite difference grid. So there will only be a small amount of measurements available to reconstruct $\mu_a$ with an observation domain that is not closed. We denote by $p^{\text{obs}}$ the acoustic wave measurements in $[0, T] \times \omega$.

Finally, we will consider 4 light sources $S_k$, $k = 1, \dots, 4$, of the diffusion equation (1). Those sources are located inside $Ω$ near it's boundary so that the diffusion equation approximation is valid. Depending on the numerical experiments, the 4 light sources will either all be placed on the left side of $Ω$ or evenly distributed near $\partial Ω$, with the same angular positions as the $\omega_i$, $i \in \{N,E,S,W\}$. Each light source will give rise to a set of measurements $p_k^{\text{obs}}$ on $[0, T] \times \omega$. This idea of considering multiple sources has already been explored in [3] for instance.

### IV. OPTIMAL CONTROL PROBLEM

As said in the introduction, the main difference between our approach and the classical ones is that we don't consider 2 distinct inverse problems: one to reconstruct the initial condition $H$ of the wave equation (2), and another one to reconstruct the $\mu_a$ map from this initial condition. We only consider one inverse problem to reconstruct the $\mu_a$ map directly from the pressure measurements and the 2 coupled equations (1) and (2). The problem we consider is then the following optimal control problem

$$ \min_{\mu_a \in \mathcal{W}_{\text{ad}}} \frac{1}{2} \sum_{k=1}^{s} \int_{[0,T] \times \omega_k} (p_k[\mu_a] - p_k^{\text{obs}})^2 + \frac{\alpha}{2} \int_{\Omega} (\mu_a - \mu_a^0)^2, \quad (3) $$

where $p_k[\mu_a](t,x)$, $k = 1, \dots, s$ is the solution of the coupled system for the given light source $S_k$ and the control $\mu_a$. The integer $s$ is the number of performed experiments, that is the number of different light sources $S_k$ to consider (4 in our case). And the admissible control set is

$$ \mathcal{W}_{\text{Ad}} = \{\mu_a \in L^{\infty}(\Omega) | \mu_a^{\min} \le \mu_a(x) \le \mu_a^{\max}, \text{ a.e. in } \Omega\} \quad (4) $$

The reason for the square on the data fitting term $(p_k[\mu_a] - p_k^{\text{obs}})^2$ is physical because it is then akin to an energy term, as well as mathematical since it adds regularity. The penalization term is added to ensure the existence of a solution to the problem. Note however that it would be interesting to use a less smooth penalization term such as an $L^1$-norm or the total variation.

To see problem (3) in the framework of optimal control problem, one has to see $\mu_a$ as the control and equations (1) and (2) as state equations. In [5], we considered a related optimal control problem but with a non-stationnary Fluence equation and with $\mu = (\mu_a, D)$ as the control. From Theorem 3.1 in [5], and with some modifications because we only consider $\mu_a$ as control and not $(\mu_a, D)$, we can deduce the following result.

**Theorem 1:** Assume that $\alpha > 0$. Then, problem (3) has at least one solution $\bar{\mu}_a$.

Additionnaly, in [6], the authors prove the following uniqueness result in a very similar case.

**Theorem 2:** Assume that $\alpha$ is large enough, then problem (3) has a unique solution $\bar{\mu}_a$.

Note that in [6], the Fluence equation is a non stationary diffusion equation and the large enough statement is explicitly quantified.

In [5], we derived first order necessary conditions for a $\bar{\mu}$ to be optimal. In this simplified framework, we can still derive first order necessary conditions and they lead to 2 adjoint equations, one for the wave equation (2), the other for the diffusion equation (1). We denote by $q^p$ in $[0, T] \times \mathcal{B}$, and $q^\Phi$ in $Ω$ the adjoint state to the pressure $p$ and to the Fluence $\Phi$. The adjoint state to the pressure $q^p$, satisifes the
following equation

$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial^2 q^p}{\partial t} - v_s^2 \Delta q^p = \mathbb{1}_{\omega}(p - p^{\text{obs}}), \quad \text{in } (0, T) \times \mathcal{B}, \\[1.5ex]
\displaystyle q^p([0, T], \partial \mathcal{B}) = 0 \\[1.5ex]
\displaystyle q^p(T, \mathcal{B}) = \frac{\partial q^p}{\partial t}(T, \mathcal{B}) = 0
\end{array}
\right.
\qquad (5)
$$

It is important to note that the initial condition of this
equation is given at time *T* and not time 0.

The adjoint of the Fluence, $q^\Phi$, satisfies the following
equation

$$
\left\{
\begin{array}{l}
\mu_a q^\Phi - \operatorname{div}(D \nabla q^\Phi) = - \frac{\partial q^p}{\partial t} (0, \cdot) \Gamma \mu_a \quad , \text{ in } \Omega \\
q^\Phi (\partial \Omega) = 0
\end{array}
\right.
\qquad (6)
$$

As mentioned before, we use the subscript *k* to refer to the
*k*-th experiment, with source $S_k$. That is we denote by $q_k^p$ the
adjoint state to $p_k$ and $q_k^\Phi$ the adjoint state to $\Phi_k$.

Again, from [5], we can deduce the following necessary
conditions for a control $\bar{\mu}_a$ to be a solution of (3).

Theorem 3: Assume $\bar{\mu}_a$ is an optimal solution to problem (3). Then, there exists $q_k^p, q_k^\Phi, k = 1, \dots, s$ such that

* The 2s equations (2) for the pressure and (1) for the Fluence are satisfied (with the s sources $S_k$, $k = 1, \dots, s$)

* The 2s adjoint state equations (5)-(6) are satisfied by $q_k^p$ and $q_k^\Phi$ respectively, for $k = 1, \dots, s$.

* For every $\mu_a \in \mathcal{U}_{\text{ad}}$, the following inequality holds

$$
\begin{equation}
\begin{aligned}
\langle -\sum_{k=1}^{s} (\partial_t q_k^p(0, \cdot)\Gamma + q_k^\Phi)\Phi_k & \\
& \quad + \alpha(\bar{\mu}_a - \mu_a^0) \rangle, \mu_a - \bar{\mu}_a \rangle_{L^2(\Omega)} \ge 0.
\end{aligned}
\tag{7}
\end{equation}
$$

V. NUMERICAL RESULTS

To illustrate this approach, we will use Theorem 3 to
reconstruct the $\mu_a$ map given in Figure 1. Note that $\mu_a$
is known in $\mathcal{B} \setminus \Omega$ and is equal to $\mu_a^0$. Also note that all
the other parameters of equations (1) and (2) are assumed
to be known, in particular the diffusion coefficient $D$, the
Grüneisen coefficient $\Gamma$ and the speed of sound $v_s$. We use
$D$ as in Figure 1, $\Gamma = 0.225$ and $v_s = 1500$ m/s. The distance
unit is normalized so that $\Omega$ is the unit square (so the distance
unit is 2 cm). The time unit is also normalized so that the
speed of sound is one. The duration $T$ of an experiment is
set so that the acoustic wave emanating from $\Omega$ has enough
time to completely exit $\Omega$ and encounter the sensors $\omega$
($T = \sqrt{2}$ after normalization). The sources $S_k$ are taken as
small (0.05 normalized distance of radius) disks inside $\Omega$ of
unit magnitude after normalization.

The numerical simulations are performed using the *Scilab*
software. We use finite difference discretization in time and
space to compute the solutions of state equations (1)-(2)
and adjoint state equations (5)-(6). Equations (2) and (5)
use a staggered scheme coupled with a perfectly matched
layer absorbing boundary for the spatial domain, see [4],
to avoid having to consider a large domain $\mathcal{B}$. Indeed, the
only requirement for domain $\mathcal{B}$ is that the acoustic wave
emanating from $\Omega$ does not rebound on $\partial \mathcal{B}$. $\Omega$ is discretized
with a 64 by 64 grid.

The following steps are used for the simulations :

1) Solve equations (1)-(2) for each source $S_k (k = 1, \dots, s)$ with the real $\mu_a$ map to obtain the measurements $p_k^{\text{obs}}$ in $[0, T] \times \omega$. Note that $\omega$ is the reunion of 4 (first example) or 2 (second example) pixels of the spatial discretization.

2) Use $\mu_a = \mu_a^0$ as initial guess for the solution of problem (3) and $\alpha = 0.1$ (for example) for the initial regularization weight.

3) Use a nonlinear conjugate gradient to cancel, up to a given precision,

$$
- \sum_{k=1}^{s} (\partial_t q_k^p(0, \cdot)\Gamma + q_k^\Phi)\Phi_k + \alpha(\bar{\mu}_a - \mu_a^0), \quad (8)
$$

with respect to $\mu_a$. Let us call $\bar{\mu}_a^\alpha$ the obtained solution.

4) If $\alpha$ is small enough, then $\bar{\mu}_a^\alpha$ is taken as our solution. Otherwise, take $\alpha = \alpha/10$ and return to step 3 with $\mu_a = \bar{\mu}_a^\alpha$ as initial guess.

Note that we do not use the projection of $\mu_a$ into $\mathcal{W}_{ad}$ so that (7) is equivalent to the cancellation of (8). Also note that, in order to compute (8) we need to first solve equation (1) and use the result $\Phi_k$ as the initial condition of equation (2) that is solved from time 0 to T and gives $p_k$. Then, equation (5) is solved from time T to 0, using the source $\mathbb{1}_{\omega}(p_k^{\text{obs}} - p_k)$ and gives $q_k^p$. Finally, $q_k^\Phi$ is computed from equation (6) and $\partial_t q_k^p(0, \cdot)$.

We will use two numerical examples. The first, in section V-A considers the case of 4 simulated experiments in which the 4 sources are all in the left side of $\Omega$. This choice will lead to the right side of $\Omega$ not being well enough illuminated. The second example, in section V-B considers the case of 4 simulated experiments in which the 4 sources are evenly distributed in $\Omega$, thus illuminating it more homogeneously.

A. 4 LEFT SIDED SOURCES

In this simulation, we consider 4 sources, located on
the left side of Ω. Figure 2 gives the sum, in log-scale
of the simulated initial conditions to equation (2), that is
$\sum_{k=1}^{4} \Gamma \mu_a \Phi^k$ on Ω. From this Figure, one can see the position
of the 4 sources, they correspond to the highest values.

The measurements are performed thanks to 4 sensors
located at $\omega_{N,E,S,W}$ from Figure 1. The computation of
the solution of problem (3) is done following the steps
introduced above, from $\alpha = 0.1$ to $\alpha = 10^{-6}$. The solution
$\bar{\mu}_a$ found is pictured on Figure 3. The regions outlined in
white represent the location of the real $\mu_a$ changes, not
of the reconstructed one. In this reconstruction of $\mu_a$, the
left ellipse is appropriately found, but the right disk does
not even start to show. This absence is not surprising since
the light sources barely illuminate the place where this
disk is. The reconstruction of the left ellipse, though not
perfect, is encouraging considering the small number of
measurements used. Note that we cannot expect to have a
$\mu_a$ with sharp changes because of the choice of the $L^2$ norm
as regularization. In order to have a $\mu_a$ solution of problem
(3) with sharp changes we could for instance consider a total
variation regularization, see [2].
<figure><img src="image_3.png" /><figcaption>Fig. 2. Sum of $\log_{10}$ of initial condition of Wave equation (2) for 4 sources distributed on the left ($\sum_{k=1}^{s}(\log_{10}\Gamma\mu_a\Phi^k + 10^{-6})$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. $\mu_a$ reconstruction for 4 sources distributed on the left and 4 sensors.</figcaption></figure>

Moreover, the maximum value of $\mu_a$ in the ellipse is preserved. Note that the shape of the reconstructed $\mu_a$ hints at a wave originating from the sensor $\omega_W$.

## B. EVEN DISTRIBUTION OF 4 SOURCES

In this simulation, we consider 4 evenly distributed sources. Figure 4 gives the sum, in log-scale, of the simulated initial conditions of equation (2), that is $\sum_{k=1}^{4}\Gamma\mu_a\Phi^k$ on $\Omega$. Again, the 4 highest values in Figure 4 give the locations of the 4 sources. We can see that the body is illuminated more homogeneously than in the previous example. In particular, the left ellipse of $\mu_a$ as well as its right disk are well illuminated.

The measurements are performed thanks to only 2 sensors, located at $\omega_{E,W}$ from Figure 1. The computation of the solution is done using the same steps as in the previous example. This solution is pictured on Figure 5. Again, the real $\mu_a$ changes are outlined in white. This reconstruction of $\mu_a$ shows both the left ellipse and the right disk that are accurately located. Again, the $\mu_a$ hints at the location of the

<figure><img src="image_5.png" /><figcaption>Fig. 4. Sum of $\log_{10}$ of initial condition of Wave equation (2) for 4 sources distributed evenly ($\sum_{k=1}^{s}(\log_{10}\Gamma\mu_a\Phi^k + 10^{-6})$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. $\mu_a$ reconstruction for 4 evenly distributed sources and 2 sensors.</figcaption></figure>

2 sensors since we can see waves in it.

Note that if in the numerical solving, we do not use the known diffusion coefficient $D$ but a uniform one equal to $D^0$ everywhere, the reconstruction is barely different. By this, we mean that for the generation of the measurements we use the non uniform $D$ (with a $\delta D$ disk at the center of $\Omega$) but for the reconstruction algorithm we use a uniform $D$.

# VI. CONCLUSIONS AND FUTURE WORKS

## A. Conclusions

In this work, we gave new numerical results stemming from the optimal control formulation of a Photoacoustic tomography problem introduced in [5]. We showed that we are able to properly reconstruct the absorption coefficient $\mu_a$ from a small set of measurements. It is expected that using more measurements will lead to a sharper reconstruction. This approach as many advantages since it is flexible in term of the observation domain. It can also be used with variable
speed of sound or other parameters like the Grüneisen one or the Diffusion.

## B. *Future Works*

Many directions can be taken to build on those results. One could consider a more complete inverse problem in which $\mu_a$ and $D$ are unknown. Indeed, this was the framework in [5] but the first numerical simulations showed that the reconstruction of $D$ does not yet yield good results. We can however expect that by multiplying the measurements and the number of different sources (and thus experiments), we should be able to improve this reconstruction of $D$. Another direction would be to change the regularization function on $\mu_a$ to allow for sharp changes of the solution. A good candidate for such a regularization term would be the total variation of $\mu_a$, see [2].

## VII. ACKNOWLEDGMENTS

The authors gratefully acknowledge the contribution of Agence Nationale pour la Recherche, through ANR project AVENTRURES, ANR-12-BS01-0001.

## REFERENCES

[1] S.R. Arridge, Optical tomography in medical imaging, *Inverse Problems*, 15, no. 2, R41, 1999

[2] H. Attouch, G. Buttazzo and G. Michaille, Variational analysis in Sobolev and BV spaces: applications to PDEs and optimization, *MPS-SIAM series on optimization*, 2006.

[3] G. Bal and K. Ren, Multi-source quantitative photoacoustic tomography in a diffusive regime, *Inverse Problems*, 27, 2011.

[4] J. P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, *Journal of Computational Physics*, Vol. 114, Is. 2, October 1994.

[5] M. Bergounioux, X. Bonnefond, T. Haberkorn and Y. Privat, An optimal control problem in Photoacoustic Tomography, *M3AS: Mathematical Models and Methods in Applied Sciences*, Vol. 24, No. 12, November 2014.

[6] M. Bergounioux and E. L. Schwindt. On the uniqueness and stability of an inverse problem in photo-acoustic tomography, *Journal of Mathematical Analysis and Applications*, 431, 2015.

[7] D. Finch, M. Haltmeier and R. Rakesh, Inversion of spherical means and the wave equation in even dimensions, *SIAM J. Appl. Math.*, Vol. 68, No. 2, 2007.

[8] M. Haltmeier, T. Schuster and O. Scherzer, Filtered backprojection for thermoacoustic computed tomography in spherical geometry, *Math. Meth. Appl. Sci.*, Vol. 28, No. 16, 2005.

[9] Y. Hristova, P. Kuchment and L.Nguyen, Reconstruction and time reversal in thermoacoustic tomography in acoustically homogeneous and inhomogeneous media, *Inverse Problems*, Vol. 24, No. 5, 2008.

[10] Y. Hristova, Time reversal in thermoacoustic tomography: an error estimate, *Inverse Problems*, Vol. 25, No. 5, 2009.

[11] N. Song, C. Deumié and A. Da Silva, Considering sources and detectors distributions for quantitative photoacoustic tomography, *Biomedical Optics Express*, Vol. 5, No. 11, Novembre 2014.

[12] M. Xu and L.V. Wang, Universal back-projection algorithm for photoacoustic computed tomography, *Physical Review E*, 71, 2005.