An optimal control approach to photoacoustic tomography

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 noninvasive, 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 phenomemon 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 phenomemon.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 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 Where Φ(x) is the fluence, µ a (x) the unknown absorption coefficient, D(x) the known diffusion coefficient and S(x) the known source.The domain Ω is the illuminated body.
The second effect is the thermal expansion that, due to it's 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) = Γ(x)µ a (x)Φ(x) that is used as the source of a wave equation Where p(t, x) is the cumulative pressure, Γ(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 B contains Ω and is large enough so that the acoustic wave initially in Ω does not reach ∂ B.

III. INVERSE PROBLEM
In this work, we want to reconstruct the absorption coefficient µ a in Ω, from acoustic wave measurements outside Ω (in B \ Ω) 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.The domain in which µ a is unknown is Ω, a 2 by 2 cm square.The domain B simply has to be large enough to allow the acoustic wave to travel without rebounding on ∂ B, this will be numerically taken care of with appropriate boundary conditions on ∂ B.
In B ⊃ Ω, µ a = µ 0 a = 0.1 cm −1 except in an ellipse on the left side and a disk on the right side of Ω where µ a = µ 0 a + 10%, as pictured in Figure 1.The diffusion coefficient D is assumed to be equal to D 0 = 0.03 cm except on a disk centered in Ω where it is equal to D 0 + 10%.The measurements are taken on 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 µ a with an observation domain that is not closed.We denote by p obs the acoustic wave measurements in [0, T ] × ω.
Finally, we will consider 4 light sources S k , k = 1, • • • , 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 ∂ Ω, with the same angular positions as the ω i , i ∈ {N, E, S,W }.Each light source will give rise to a set of measurements p obs k on [0, T ] × ω.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 µ a map from this initial condition.We only consider one inverse problem to reconstruct the µ 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 where is the solution of the coupled system for the given light source S k and the control µ 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 The reason for the square on the data fitting term 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 µ 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 µ = (µ a , D) as the control.From Theorem 3.1 in [5], and with some modifications because we only consider µ a as control and not (µ a , D), we can deduce the following result.
Theorem 1: Assume that α > 0.Then, problem (3) has at least one solution μa .Additionnaly, in [6], the authors prove the following uniqueness result in a very similar case.
Theorem 2: Assume that α is large enough, then problem (3) has a unique solution μa .Note that in [6], the Fluence equation is a non stationnary diffusion equation and the large enough statement is explicitely quantified.
In [5], we derived first order necessary conditions for a μ 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 ] × B, and q Φ in Ω the adjoint state to the pressure p and to the Fluence Φ.The adjoint state to the pressure q p , satisifes the following equation 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 Φ , satisfies the following equation 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 p k the adjoint state to p k and q Φ k the adjoint state to Φ k .Again, from [5], we can deduce the following necessary conditions for a control μa to be a solution of (3).
Theorem 3: Assume μa is an optimal solution to problem (3).Then, there exists q p k , q Φ k , k = 1, • • • , 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, • • • , s) • The 2s adjoint state equations ( 5)-( 6) are satisfied by q p k and q Φ k respectively, for k = 1, • • • , s. • For every µ a ∈ U ad , the following inequality holds To illustrate this approach, we will use Theorem 3 to reconstruct the µ a map given in Figure 1.Note that µ a is known in B \ Ω and is equal to µ 0 a .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 Γ and the speed of sound v s .We use D as in Figure 1, Γ = 0.225 and v s = 1500 m/s.The distance unit is normalized so that Ω 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 Ω has enough time to completely exit Ω and encounter the sensors ω (T = √ 2 after normalization).The sources S k are taken as small (0.05 normalized distance of radius) disks inside Ω 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 B. Indeed, the only requirement for domain B is that the acoustic wave emanating from Ω does not rebound on ∂ B. Ω is discretized with a 64 by 64 grid.
The following steps are used for the simulations : 1) Solve equations ( 1)-( 2) for each source with the real µ a map to obtain the measurements p obs k in [0, T ] × ω.Note that ω is the reunion of 4 (first example) or 2 (second example) pixels of the spatial discretization.2) Use µ a = µ 0 a as initial guess for the solution of problem (3) and α = 0.1 (for example) for the initial regularization weight.
3) Use a nonlinear conjugate gradient to cancel, up to a given precision, with respect to µ a .Let us call μα a the obtained solution.4) If α is small enough, then μα a is taken as our solution.Otherwise, take α = α/10 and return to step 3 with µ a = μα a as initial guess.Note that we do not use the projection of µ a into U 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 Φ 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 1 ω (p obs k − p k ) and gives q p k .Finally, q Φ k is computed from equation ( 6) and ∂ t q p k (0, •).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 Ω.This choice will lead to the right side of Ω 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 Ω, 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 ∑ 4 k=1 Γµ a Φ 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 ω N,E,S,W from Figure 1.The computation of the solution of problem (3) is done following the steps introduced above, from α = 0.1 to α = 10 −6 .The solution μa found is pictured on Figure 3.The regions outlined in white represent the location of the real µ a changes, not of the reconstructed one.In this reconstruction of µ 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 µ a with sharp changes because of the choice of the L 2 norm as regularization.In order to have a µ a solution of problem (3) with sharp changes we could for instance consider a total variation regularization, see [2].Moreover, the maximum value of µ a in the ellipse is preserved.Note that the shape of the reconstructed µ a hints at a wave originating from the sensor ω 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 ∑ 4 k=1 Γµ a Φ k on Ω.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 µ a as well as its right disk are well illuminated.
The measurements are performed thanks to only 2 sensors, located at ω 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 µ a changes are outlined in white.This reconstruction of µ a shows both the left ellipse and the right disk that are accurately located.Again, the µ a hints at the location of the  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 δ D disk at the center of Ω) but for the reconstruction algorithm we use a uniform D.

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 µ 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 µ 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 µ a to allow for sharp changes of the solution.A good candidate for such a regularization term would be the total variation of µ a , see [2].

Fig. 1 .
Fig. 1.Setting of the 2 domains and distribution of µ a and D.

Fig. 3 .
Fig. 3. µ a reconstruction for 4 sources distributed on the left and 4 sensors.