# Modelling Human Kinetics and Kinematics during Walking using Reinforcement
Learning

Visak Kumar

###### Abstract

In this work, we develop an automated method to generate 3D human walking
motion in simulation which is comparable to real-world human motion. At the
core, our work leverages the ability of deep reinforcement learning methods to
learn high-dimensional motor skills while being robust to variations in the
environment dynamics. Our approach iterates between policy learning and
parameter identification to match the real-world bio-mechanical human data. We
present a thorough evaluation of the kinematics, kinetics and ground reaction
forces generated by our learned virtual human agent. We also show that the
method generalizes well across human-subjects with different kinematic
structure and gait-characteristics.

##  I INTRODUCTION

Assistive devices such as exoskeletons have garnered a lot of research
interest because they can bring numerous benefits to injury/disability
rehabilitation, elderly care, human-augmentation for people doing heavy
lifting at work, exercise and military applications. However, there a few
challenges we need to overcome before assistive devices can become more
prevalent in our society. Exoskeletons work by directly applying forces on the
joints, hence safety during operation becomes paramount when designing control
policies for such devices. Often, testing directly on a human-subject is not
recommended.

Simulation has played a key role in offloading the burden of testing control
policies in the real-world with human subjects. However, simulating the real-
world accurately has been a long-standing challenge in robotics. Modelling
human motion is especially difficult because of the high-dimensional nature of
the state space and the natural variability that exists between different
human-subjects.

Our focus is on modelling human motion during locomotion tasks such as walking
and running. Humanoid locomotion is governed by highly non-linear and
discontinuous dynamics and large number of unknown quantities such as joint
damping, segment Center of Mass (COM), segment inertia, friction coefficient,
muscle activation and muscle-tendon interaction. These are extremely
challenging to measure. Traditional methods to identify simulation models,
commonly known as system identification, are often insufficient when the
system in question is high-dimensional. System identification works by
collecting data on the real-system and then inferring model parameters
parameters based on the data collected. For a high-dimensional system,
collecting data in the task-relevant state space can be quite challenging.

The difficulty in modelling human locomotion has led some researchers to focus
on reduced-order models, such as Linear-inverted pendulum models (LIPM) (or
even 3D LIPM) [1], to explain human locomotion. Although computational more
efficient, these models come at the cost of accuracy and hence are limited in
their application.

Reinforcement learning (RL) and imitation learning have taken center stage in
recent years in developing control algorithms to imitate human motion. In
particular, Peng et al [2] and Yu et al [3] learn complex motor skills like
walking, running, jumping and even performing a backflip. In addition to
learning complex motion, policies trained using RL are also robust to
variations in the underlying dynamics : variations caused by different
dynamical parameters [4] as well as shape and size of the human-subject
[Won:2019].

In this work, we leverage the following two advantages of RL, (1) Ability to
learn complex motions in high-dimensional space and (2) robustness to
variability in underlying dynamics, to learn control policies for human
walking that are bio-mechanically accurate. We validate our approach by
comparing the joint angles, joint moments and ground reaction forces generated
by the policy to real-world data collected in human-subject experiments.

##  II Related Work

Since this work aims to leverage tools and data which are outcomes of research
in two different fields - Deep Reinforcement Learning (DRL) and Biomechanics.
The literature review is organized into two sections, each describes the
current state of the art in each field and what aspects needs improvement.

###  II-A Deep Reinforcement Learning

DRL has seen remarkable success in recent years for learning complex tasks
ranging from video games to controlling robots [5, 6, 7]. In particular, DRL
has been successful in the field of developing locomotion controllers [8, 9,
4, 2]. However, these methods are seldom validated by comparing it to real
human movement data generated by experiments. In Peng et al [2], a policy
optimized using DRL algorithm was able to learn exceptional skills like
walking, running , jumping and even doing a backflip with the help of motion
capture data. But, the resulting policy, while visually pleasing, is not
validated with real-world human data. Similarly, in Yu et al [9], several
locomotion skills like walking and running at different speeds was learned
from scratch using a curriculum learning approach. However, as was the issue
with [2], this approach is not validated with real data. Additionally, the
above mentioned approach only adopts symmetry and low energy as metrics to
enforce in the learned walking strategy, but some research in biomechanics
points towards an existence of more complex relationships that gives rise to
walking. For example, Wang et al [10] identified that the Center-of-mass
velocity has a direct linear relationship to step-length. In Kumar et al [11],
preliminary comparison of the walking motion learned using policy optimization
and real-world human subject experiments were made. However, only joint
kinematics and foot-step lengths were used as a metric for validation. We need
a more thorough evaluation using joint kinematics, kinetics, ground reaction
forces and foot step lengths across different individuals to be more confident
about the simulation model.

###  II-B Biomechanics

To study human motion during walking and recovery, biomechanics researchers
often adopt an experimental approach. First, data is collected in the real-
world, then control policies are synthesized using the real-world data.
Winters et al [12] was among the first in the field to study human gait, and
the data published in this work remains relevant to this day. Wang et al [10]
and Hof et al [13] performed perturbation experiments and identified important
relationships between COM velocity, step-lengths, center of pressure, stepping
vs ankle strategy, etc.. We aim to leverage the finding of this research to
validate some of the results. In Joshi et al [1], a balance recovery
controller was derived using the results reported in [10], however, they use a
3D Linear Inverted pendulum model to approximate the human dynamics. A 3D LIPD
does not capture the dynamics fully, for example, angular momentum about the
center of mass. Most relevant to our work, Antoine et al [14], used a direct-
collocation trajectory optimization to synthesize a walking controller for a
3D musculo-skeletal model in OpenSim (OpenSim gait2392 model). The gait
generated by the controller closely matched experimental data. The proposed
method, relies on understanding the basic principles that lead to walking,
such as minimizing metabolic cost, muscle activations,etc. However, the
proposed solution enforces left-right symmetry, which works for walking, but
is not ideal for disturbance recovery. Hence its unclear how well this
approach will perform when there is an external disturbance to the human.

##  III Method

We propose a framework to automate the process of developing bio-mechanically
accurate 3D Human walking policies in simulation. The difference between the
dynamics of a simulated human agent and a real human subject is caused by
multiple factors. Some quantities such as mass, height, etc.. can be easily
measured. However other quantities such as dynamical parameters: joint
damping, ground friction, lower-limb joint axis location and the accurate
segment lengths can be challenging to estimate. Note that while the segment
lengths can be easily measured on a real human, the bio mechanical data sets
usually do not provide this information.

We present an iterative approach to develop accurate walking model in which,
first, we use Deep Reinforcement Learning (DRL) to learn a walking policy
using a nominal dynamical model. During training, domain randomization is used
to ensure the policy is robust to dynamical and kinematic parameters. Second,
once we have a policy to generate walking motion, we perform an optimization
step to identify the optimal parameters that explain real-world walking
motion, these two steps are repeated until convergence. To validate our
simulated models, we compare the gait characteristics such as joint
kinematics, kinetics and ground reaction forces generated by our policy to
real-world data collected with 5 human participants.

We formulate this problem of learning human walking as a Markov Decision
Processes (MDPs), $(\mathcal{S},\mathcal{A},\mathcal{T},r,p_{0},\gamma)$,
where $\mathcal{S}$ is the state space, $\mathcal{A}$ is the action space,
$\mathcal{T}$ is the transition function, $r$ is the reward function, $p_{0}$
is the initial state distribution and $\gamma$ is a discount factor. We take
the approach of model-free reinforcement learning to find a policy $\pi$, such
that it maximizes the accumulated reward:

$J(\pi)=\mathbb{E}_{\mathbf{s}_{0},\mathbf{a}_{0},\dots,\mathbf{s}_{T}}\sum_{t=0}^{T}\gamma^{t}r(\mathbf{s}_{t},\mathbf{a}_{t}),$

where $\mathbf{s}_{0}\sim p_{0}$, $\mathbf{a}_{t}\sim\pi(\mathbf{s}_{t})$ and
$\mathbf{s}_{t+1}=\mathcal{T}(\mathbf{s}_{t},\mathbf{a}_{t})$.

We denote the human walking policy as $\pi_{h}(\mathbf{a}_{h}|\mathbf{s}_{h})$
where $\mathbf{s}_{h}$, $\mathbf{a}_{h}$ represent the states and actions,
respectively.

##  IV Data and Simulation Environment

Fukuchi et al [15] published open-source biomechanics data of 38 human
subjects walking overground and on treadmill. We use data of 10 subjects (for
proof-of-concept) with different physical characteristics such as mass,
height, leg-length, gender and speed of walking. For each subject, this
dataset provides motion capture marker data, and individual gait
characteristics like joint angles, joint moments and ground reaction forces.

We use DART physics [16] engine as our simulation environment in which the
virtual human agent learns to walk. We first create 10 virtual agents with the
same physical characteristics as the human subjects shown in table I. We
choose hunt-crossley contact model to generate contact forces between two
rigid bodies such as feet and ground. Hunt-crossley model is widely used in
biomechanics research and is also one of the contact models used in OpenSim
simulator [17].

Subject | Mass (kg) | Height (cm) | Age | Speed (m/s)
---|---|---|---|---
1 (M) | 74 | 172.40 | 25 | 1.25
2 (M) | 52.9 | 166.80 | 22 | 1.40
3 (F) | 48.8 | 158 | 24 | 1.15
4 (M) | 61.5 | 180.70 | 22 | 1.28
5 (F) | 153 | 64 | 31 | 1.09
6 (M) | 69.85 | 155 | 38 | 1.3
7 (F) | 64.6 | 151.50 | 57 | 0.91
8 (M) | 63.3 | 175 | 71 | 0.57
9 (F) | 46.05 | 149.20 | 63 | 0.80
10 (M) | 66.35 | 155.50 | 84 | 0.60
TABLE I: Human subject physical characteristics. Figure 1:  We model a
31-Degree of Freedom(DoF) humanoid in PyDart.

###  IV-A Human Walking Policy

We take a model-free reinforcement learning approach to develop a human
locomotion policy $\pi_{h}(\mathbf{a}_{h}|\mathbf{s}_{h})$. To achieve natural
walking behaviors, we train a policy that imitates the human walking reference
motion (motion capture data from the open-source data base) similar to Peng et
al. [2]. The human 3D model (agent) consists of $25$ actuated joints with a
floating base as shown in Figure 1. This gives rise to a $71$ dimensional
state space
$\mathbf{s}_{h}=[\mathbf{q},\mathbf{\dot{q}},\mathbf{v}_{com},\boldsymbol{\omega}_{com},\phi,\mu]$,
including joint positions, joint velocities, linear and angular velocities of
the center of mass (COM), and a phase variable($\phi$) that indicates the
target frame in the motion clip. Vector $\mu$ ($R^{1}4$) described as

$\mu\ =\ [\beta,\sigma,f_{l},s_{l},t_{l},f_{j},h_{j},s_{j}]$

where $\beta$ is joint damping, $\sigma$ is friction coefficient between the
foot and the ground. $f_{l},s_{l},t_{l}$ and foot, shin and thigh segment
lengths respectively and $f_{j},h_{j},s_{j}$ are their corresponding joint
axis location expressed in the robot root coordinate frame. These quantities
are randomized during training so that the learned policy can be robust to
variations in these quantities.

The action determines the target joint angles $\mathbf{q}_{t}^{target}$ of the
proportional-derivative (PD) controllers, deviating from the joint angles in
the reference motion:

$\mathbf{q}^{target}_{t}=\hat{\mathbf{q}}_{t}(\phi)+\mathbf{a}_{t},$ (1)

where $\hat{\mathbf{q}}_{t}(\phi)$ is the corresponding joint position in the
reference motion at the given phase $\phi$. Our reward function is designed to
imitate the reference motion:

$r_{h}(\mathbf{s}_{h},\mathbf{a}_{h})=w_{q}(\mathbf{q}-\hat{\mathbf{q}}(\phi))+w_{c}(\mathbf{c}-\hat{\mathbf{c}}(\phi))\\\
+w_{e}(\mathbf{e}-\hat{\mathbf{e}}(\phi))-w_{\tau}||\boldsymbol{\tau}||^{2}$
(2)

where $\hat{\mathbf{q}}$, $\hat{\mathbf{c}}$, and $\hat{\mathbf{e}}$ are the
desired joint positions, COM positions, and end-effector positions from the
reference motion data, respectively. The reward function also penalizes the
magnitude of torque $\boldsymbol{\tau}$. We use the same weight $w_{q}=5.0$,
$w_{c}=2.0$, $w_{e}=0.5$, and $w_{\tau}=0.005$ for all experiments. We also
use early termination of the rollouts, if the agent’s pelvis drops below a
certain height or if the base rotates about any axis beyond a threshold, we
end the rollout and re-initialize the state.

We exert random forces to the agent during policy training. Each random force
has a magnitude uniformly sampled from $[0,800]\ N$ and a direction uniformly
sampled from [-$\pi/2$,$\pi/2$], applied for $50$ milliseconds on the agent’s
pelvis in parallel to the ground. The maximum force magnitude induces a
velocity change of roughly $0.6$m/sec. This magnitude of change in velocity is
comparable to experiments found in literature such as [10],[18] and [13]. We
also randomize the time when the force is applied within a gait cycle.
Training in such a stochastic environment is crucial for reproducing the human
motion. We represent a human policy as a multi-layered perceptron (MLP) neural
network with two hidden layers of $128$ neurons each. The formulated MDP is
trained with Proximal Policy Optimization (PPO) [6].

###  IV-B Parameter Identification

As described in the previous section, the policy $\pi_{h}$ is a function of
vector $\mu$. So, once the policy is trained we can search over the parameter
$\mu$ to identify the values which best explains reference motion data. To do
this, we perform an optimization step using CMA-ES algorithm [19]. We use the
same reward function as described in equation 2. In this optimization process,
at the beginning of each trajectory we set the value of $\mu$, then generate
the motion by executing the policy $\pi_{h}$, at the end of each generated
trajectory the reward is computed using the equation 3. The goal is to
maximize the sum of reward accumulated along the trajectory generated by the
policy $\pi_{h}$.

$\mu^{*}=\arg\max_{\mu}^{(s,a)~{}\pi}\sum_{t=0}^{T}(r_{h}^{t}(s_{h}^{t},a_{h}^{t}))$
(3)

###  IV-C Error Threshold

We iterate between policy optimization and parameter identification until the
error, defined in equation 4, is less than a predefined threshold $\kappa$.

$\epsilon=\sum_{t=0}^{T}||q-\hat{q}||+||\tau-\hat{tau}||+||GRF-\hat{GRF}||)$
(4)

Here, $T$ is the time taken to complete one gait cycle, $q$ is the joint
angles generated by the policy and $\hat{q}$ is the ground truth joint angles
for that particular human subject. Similarly, $\tau$ is the joint moments and
$GRF$ is the ground reaction forces (vector of dimension $R^{3}$).
$\hat{\tau}$ and $\hat{GRF}$ are the ground truth joint moments and ground
reaction forces measured in the real-world.

The procedure is outlined in algorithm 1.

1: for human subject i:1…n do

2: Input: Dataset $D_{i}$ \- Motion capture trajectory of Human-subject $i$,
$\kappa$ \- error threshold

3: Initialize CMA optimizer with generation size 8

4: Initialize $\pi_{i}(a|s,\mu)$, $V_{i}(s)$

5: while $\epsilon>\kappa$ do

6: Optimize policy $\pi$ using reward function 2

7: Optimize $\mu^{*}$ using equation 3

8: Compute error $\epsilon$ using 4

9: end while

10: end for

11: return $\mu^{*}$ and $\pi$

Algorithm 1 Algorithm

##  V Results

We trained walking policies for 5 human subjects (first 5 in table I). We do
the following evaluation to validate our approach,

  1. 1.

Ablation study : How well does the joint angles, joint moments and ground
reaction forces generated by our method compare to a baseline method which
only uses RL? In other words, how much of a change does the second
optimization step benefit the accuracy? We use Root-mean squared error as the
metric for comparison.

  2. 2.

We highlight preliminary analysis on the generalizability of the approach by
comparing the joint angles of 5 different human subjects.

###  V-A Ablation Study

This ablation study is done only for one subject. In figure 2, we illustrate
the benefit of our method, the additional optimization step has a clear effect
on the RMSE error. Further, in figure 3, we also show the lower-limb joint
angles before and after the CMA-ES optimization step. The ankle joint in
particular matches the real-world data more closely. This highlights the
ability of control policies trained using RL to be robust to variation in the
dynamical parameters.

Figure 2: RMSE error of our method compared to baseline approach which uses
just RL to train the walking policy. We use 6 different initialization seeds
for the policy, the bar plot indicates the RMSE error mean and standard
deviation for the 6 policies for the same subject.

Figure 3: Comparison of lower-limb joint angles during one gait cycle. Top:
Joint angles after RL step, the ankle joints are noticeably different although
the curve profile is similar. Bottom: After CMA-ES optimization step, the
joint angles match better with real-world data.

###  V-B Joint moments, ground reaction forces and work loops

We also compare the joint moments and ground reaction forces generated by the
policy to ground truth data by our method. The comparisons are illustrated in
figures 5 and 6 respectively.

In addition to this, we also compare the torque loop at the hip joint and
compare it to the bio-mechanical data reported in [12] (a long-standing gold
standard). Torque loop is a plot of the torques generated at the hip joint in
the y-axis and the joint angle in the x-axis during one-gait cycle. The arrows
indicate the direction of movement and the points $[0,50,100]$% of the gait
cycle match fairly well.

Figure 4: Comparison of torque loops of a typical trajectory generated by our
policy and human data reported by [12] at the hip of stance leg during a gait
cycle. The green dots indicate the start and the black dots indicate 50% of
the gait cycle. The arrows show the progression of the gait from 0% to 100%.
Figure 5: Joint moments vs ground truth data.

Figure 6: Ground reaction forces Top: Tangential ground reaction force, its
the force that propels a person forward while walking Bottom: Vertical ground
reaction force, this is the normal force a person experiences while walking.

###  V-C Results for five subjects

Our approach generalizes well to different human subjects. We apply our
approach to 5 different human subjects walking at 5 different speeds. Figure 7
(top) illustrates the ground truth knee joint angle profiles for these 5
different human subjects. This shows the varied nature of the joint
trajectories depending on speed and individual characteristics such as leg-
length and mass. Our approach can capture these differences, as shown in
figure 7 (bottom), the joint angles generated for the 5 different policies
match well with the ground truth data.

Figure 7: Comparison of recovery performance when perturbation is applied at
four different phases. Top: Comparison of stability region. Bottom: Comparison
of COM velocity across five gait cycles. Perturbation is applied during the
gait cycle ’p’. The increasing velocity after perturbation indicates that our
policy is least effective at recovering when the perturbation occurs later in
the swing phase.

##  VI Discussion and Conclusion

In this work, we have shown that our proposed method can generate 3D human
walking motion that matches well with real-world data. Our ablation study
indicates there is an improvement over existing RL methods to learn walking
motion. There are plenty of questions we can try to answer in future steps:

  * •

First, we would like to test our algorithm on all 42 subjects in the open
source dataset [15].

  * •

How well can we model a walking gait generated by a person with disability?

  * •

How well can we model a walking gait generated by a person wearing an
assistive device such as hip/ankle exoskeletons?

  * •

Can we leverage this human model to learn interesting assistive strategies for
an exoskeleton such as assistive walking to reduce metabolic cost of walking
or to assist recovery when there is an external perturbation such as trips or
slips?

## References

  * [1] V. Joshi and M. Srinivasan, “A controller for walking derived from how humans recover from perturbations,” 2019.
  * [2] X. B. Peng, P. Abbeel, S. Levine, and M. van de Panne, “Deepmimic: Example-guided deep reinforcement learning of physics-based character skills,” _ACM Transactions on Graphics (Proc. SIGGRAPH 2018)_ , 2018.
  * [3] W. Yu, G. Turk, and C. K. Liu, “Learning symmetric and low-energy locomotion,” _ACM Transactions on Graphics_ , vol. 37, no. 4, p. 1–12, Jul 2018. [Online]. Available: http://dx.doi.org/10.1145/3197517.3201397
  * [4] W. Yu, V. C. V. Kumar, G. Turk, and C. K. Liu, “Sim-to-real transfer for biped locomotion,” _CoRR_ , vol. abs/1903.01390, 2019. [Online]. Available: http://arxiv.org/abs/1903.01390
  * [5] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International Conference on Machine Learning_ , 2015, pp. 1889–1897.
  * [6] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [7] H. van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double q-learning,” _CoRR_ , vol. abs/1509.06461, 2015. [Online]. Available: http://arxiv.org/abs/1509.06461
  * [8] X. B. Peng, G. Berseth, K. Yin, and M. Van De Panne, “Deeploco: Dynamic locomotion skills using hierarchical deep reinforcement learning,” _ACM Transactions on Graphics (TOG)_ , vol. 36, no. 4, p. 41, 2017.
  * [9] W. Yu, G. Turk, and C. K. Liu, “Learning symmetric and low-energy locomotion,” _ACM Transactions on Graphics (Proc. SIGGRAPH 2018)_ , vol. 37, no. 4, 2018.
  * [10] Y. Wang and M. Srinivasan, “Stepping in the direction of the fall: The next foot placement can be predicted from current upper body state in steady-state walking,” _Biology Letters_ , vol. 10, no. 9, 2014.
  * [11] V. C. V. Kumar, S. Ha, G. Sawicki, and C. K. Liu, “Learning a control policy for fall prevention on an assistive walking device,” _Internation Conference of Robotics and Animation_ , 2019.
  * [12] D. A. Winter, _Biomechanics and motor control of human gait: normal, elderly and pathological_. Waterloo Biomechanics, 1991.
  * [13] A. Hof, S. Vermerris, and W. Gjaltema, “Balance responses to lateral perturbations in human treadmill walking,” _Journal of Experimental Biology_ , vol. 213, no. 15, pp. 2655–2664, 2010.
  * [14] F.Antoine, S. Gil, D. Christopher, G. Joris, J. Ilse, and F. Groote, “Rapid predictive simulations with complex musculoskeletal models suggest that diverse healthy and pathological human gaits can emerge from similar control strategies,” 2019.
  * [15] R. F. Claudine A Fukuchi and M. Duarte, “A public dataset of overground and treadmill walking kinematics and kinetics in healthy individuals,” 2016. [Online]. Available: https://doi.org/10.1080/00207540701450013
  * [16] S. Ha, “Pydart2,” 2016. [Online]. Available: https://github.com/sehoonha/pydart2
  * [17] OpenSim, “Hunt-Crossley model,” 2011.
  * [18] D. Martelli, V. Vashista, S. Micera, and S. K. Agrawal, “Direction-dependent adaptation of dynamic gait stability following waist-pull perturbations,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 24, no. 12, pp. 1304–1313, Dec 2016.
  * [19] N. Hansen, “The CMA evolution strategy: A tutorial,” _CoRR_ , vol. abs/1604.00772, 2016. [Online]. Available: http://arxiv.org/abs/1604.00772


# A high contrast and resolution reconstruction algorithm in quantitative
photoacoustic tomography††thanks: Partially supported by the BMBF-Project
iDeLIVER : Intelligent MR Diagnosis of the Liver by Linking Model and Data-
driven Processes, US National Science Foundation Grant No: 2309491, and
University of Texas at Arlington, Research Enhancement Program Grant No:
2022-605.

Anwesa Dey<EMAIL_ADDRESS>Department of Mathematics, University of
Utah, John A. Widtsoe Building 155 S 1400 E, Salt Lake City, 84112, Utah, USA.
Alfio Borzì<EMAIL_ADDRESS>Institut für Mathematik,
Universität Würzburg, Emil-Fischer-Strasse 30, 97074 Würzburg, Germany. Souvik
Roy<EMAIL_ADDRESS>Department of Mathematics, University of Texas at
Arlington, 411 S. Nedderman Drive, Arlington, 76019, Texas, USA.

###### Abstract

A framework for reconstruction of optical diffusion and absorption
coefficients in quantitative photoacoustic tomography is presented. This
framework is based on a Tikhonov-type functional with a regularization term
promoting sparsity of the absorption coefficient and a prior involving a
Kubelka-Munk absorption-diffusion relation that allows to obtain superior
reconstructions. The reconstruction problem is formulated as the minimization
of this functional subject to the differential constraint given by a photon-
propagation model. The solution of this problem is obtained by a fast and
robust sequential quadratic hamiltonian algorithm based on the Pontryagin
maximum principle. Results of several numerical experiments demonstrate that
the proposed computational strategy is able to obtain reconstructions of the
optical coefficients with high contrast and resolution for a wide variety of
objects.

Keywords: Quantitative photoacoustic tomography, nonsmooth PDE optimization,
Pontryagin maximum principle, sequential quadratic hamiltonian algorithm.

MSC: 35R30, 49J20, 49K20, 65M08, 82C31

##  1 Introduction

Photoacoustic tomography (PAT) is an emerging hybrid imaging technique that
combines two imaging modalities: optical tomography, with a large contrast of
optical parameters in an object, and ultrasound imaging, which provides high
resolution images. The PAT technique is based on the photoacoustic effect
where an object is exposed to a short-pulse optical radiation that propagates
through it and is, in part, partially absorbed. This phenomenon leads to
heating and thermal expansion of the region where the propagation occurs, so
that acoustic pulses are emitted that travel back to the boundary of the
object, where they can be measured [61, 63]. From these measured acoustic
pulses, one can obtain the initial pressure distribution and, subsequently,
optical coefficients like the absorption and diffusion, through the
formulation of inverse problems.

The reconstruction of the initial pressure distribution denoted with
$\mathcal{H}(x),~{}x\in\Omega\subset\mathbb{R}^{n}$, from the measured
acoustic pulses on the boundary of the object $\partial\Omega$, can be
formulated as an inverse source problem for the wave operator as follows: Let
$p(x,t)$ and $\mathcal{C}(x)$ denote the pressure and the speed of the
acoustic wave, respectively. Then $p$ satisfies the following equation

$\displaystyle p_{tt}-\mathcal{C}^{2}(x)\,\Delta
p=0,\quad(x,t)\in\Omega\times[0,{T}]$ (1) $\displaystyle
p(x,0)=\mathcal{H}(x),\quad x\in\Omega$ $\displaystyle\dfrac{\partial
p}{\partial t}(x,0)=0,\quad x\in\Omega$

where ${T}$ represents the time period of propagation of the acoustic waves
starting with the initial pressure field $\mathcal{H}$. Then, the time-
dependent acoustic wave signals of the form

$p(x,t),~{}\dfrac{\partial{p}}{\partial
n}(x,t)\quad(x,t)\in\partial\Omega\times[0,{T}]$ (2)

are measured at the acquisition boundary $\partial\Omega$. The inverse problem
is to recover $\mathcal{H}$ given the data in (2). This is also known as the
classical PAT inverse problem [10, 11, 12, 13, 25, 31]

The aforementioned inverse problem has been very well studied; see, e.g., [25,
31], and the solution to this problem gives the initial pressure distribution
$\mathcal{H}$ profile in the entire domain $\Omega$. However, this procedure
does not give information about the more specific properties of the object
like the optical parameters of absorption and diffusion coefficient. Thus, to
reconstruct the optical parameters from the initial pressure distribution, a
second inverse problem is solved, which is known as quantitative photoacoustic
tomography (QPAT). One of the major applications of QPAT is to provide
accurate images of the heart for the study of congenital heart diseases [32].
Other applications of QPAT include detecting cancerous tissues, thus
facilitating biomedical studies [34, 51]. Notice that there are other optical
imaging methods that have a structure similar to QPAT in the sense that these
methods make use of the photoacoustic effect to determine specific internal
features in an object. For example, the fluorescence optical diffuse
tomography (FDOT), which is used to detect fluorescing targets in tissues [46,
53].

However, in QPAT and similar modalities, the accurate reconstruction of the
optical parameters poses several theoretical and computational challenges that
are under investigation by many research groups. Our purpose is to contribute
to this research work with a novel QPAT reconstruction framework with high
resolution and high contrast, which is of paramount importance especially in
the context of biomedical imaging for detection of cancerous tissues.

Mathematically, the initial pressure distribution $\mathcal{H}({x})$ is a
function of the absorption coefficient $\sigma_{a}({x})$, of the intensity of
radiation $u({x})$, and of the Grüneisen coefficient $\Gamma({x})$, which
measures the amount of ultrasound generated by the absorbed radiation at the
point ${x}$. In general, the intensity of radiation $u$ is a function of space
and time and its evolution is governed by the radiative transfer equation
(RTE) (see e.g., [9, Section 3.2], [38, Chapter 7], [51]). However, assuming
the use of short light pulses such that the radiation propagates for a very
short time at the scale of small acoustic wave lengths, we have that the
scattering is significantly larger than absorption (see [38, Chapter 9], [9,
Section 3.3]). Thus, the photons are scattered almost uniformly in all
directions leading to an almost uniform angular distribution. This phenomenon
is prevalent inside biological tissues [30, 51] and blood [38, Chapter 9]. In
this case, it has been shown in [9, Sec 3.3] that the photon current satisfies
a steady state assumption, which results in the diffusion approximation of the
RTE equation. Therefore, in this setting, the intensity of radiation $u({x})$
is modelled by the following diffusion equation [11, 12, 13]:

$\displaystyle-\nabla\cdot D({x})\,\nabla u({x})+\sigma_{a}({x})\,u({x})$
$\displaystyle=0,\qquad\mbox{ in }\Omega,\,$ $\displaystyle u({x})$
$\displaystyle=g({x}),\quad\mbox{ on }\partial\Omega,$

where $g$ is the known intensity of the radiation that is applied at the
boundary of the object. In this framework, the initial pressure distribution
function is given by

$\mathcal{H}({x})=\Gamma({x})\,\sigma_{a}({x})\,u({x}),\quad{x}\in\Omega.$

The QPAT inverse problem requires to recover the absorption coefficient
$\sigma_{a}({x})$ and the diffusion coefficient $D({x})$, given $\mathcal{H}$
and $g$.

There are several works related to theoretical and computational frameworks
for reconstructions in QPAT; see, e.g.,[11, 15, 33, 36, 64]. In the work [13],
it has been shown that, if $\Gamma$ is given, then the coefficients
$(D,\sigma_{a})$ can be uniquely and stably reconstructed based on two
different sets of initial pressure distribution data generated by two well-
chosen boundary conditions $g_{1}$ and $g_{2}$ for two copies of the diffusion
equation given above. Also in [13], we find a method for reconstruction that
relies on the solution of a set of transport and elliptic equations. Further,
in [11, 12] it has been shown that all three parameters
$(D,\sigma_{a},\Gamma)$ cannot be reconstructed simultaneously, irrespective
of the number of illuminations unless one considers the QPAT inversion problem
at different wavelengths. The correspondence between the QPAT inversion
problem and parameter estimation with PDE’s was illustrated in [14]. Other
works on QPAT include: the multiple-source framework in [65], and its analysis
[54], a PDE-optimization framework using the radiative transfer equation [1],
the recovery of absorption coefficient of a small absorber [8], the
reconstruction of absorption and scattering coefficients [59], and the
analysis of the linearised QPAT problem [36].

In addition to the theoretical investigation of the QPAT problem, there are
several numerical reconstruction algorithms for solving related inversion
problems. In [24], the authors develop an iterative method using the radiation
map to reconstruct the optical parameters. A non-linear model-based inversion
scheme was formulated in [43] for determining chromophore concentrations in
2D. In [52], a gradient-based scheme for QPAT using the radiative transfer
equation is discussed. A new iterative reconstruction algorithm using a
different form of the radiation function and Robin-boundary data was proposed
in [47]. It appears from the simulation results in these works that existing
algorithms for QPAT have difficulties in reconstructing accurately both $D$
and $\sigma_{a}$. Furthermore, these algorithms do not perform well to
reconstruct the optical coefficients with high contrast, which is an inherent
property of tissues, thus limiting their range of applicability in biomedical
problems.

One of the major reason for the poor reconstruction of $D$ is due to the fact
that the initial pressure distribution data does not directly contain
information on $D$. On the other hand, it has been shown that there is a
relation between $D$ and $\sigma_{a}$ and the reduced scattering coefficient
$\sigma_{s}$, which is given by [6, 27]

$D=\dfrac{1}{3(\sigma_{a}+\sigma_{s})}.$ (3)

This fact shows a deficiency in the use of the diffusion approximation for
reconstructing $D$, which would require to employ the radiative transfer (or
transport) equation (RTE) [38] as the adequate model for light transport in
tissues. Therefore, although the main point in photo acoustic is a
thermodynamics phenomenon that explains how absorbed energy is transformed in
a pressure wave, also scattering is relevant for a successful reconstruction
of the diffusion coefficient, as it is illustrated by results in, e.g., [59].

However, in some applications a proportionality between $\sigma_{a}$ and
$\sigma_{s}$ appears. In particular, the following relation between
$\sigma_{a}$ and $\sigma_{s}$ has for long been described by the Kubelka-Munk
theory in the diffusive regime [41, 37, 62, 28, 48, 29]. We have

$\sigma_{s}=\dfrac{2R_{\infty}}{(1-R_{\infty})^{2}}\sigma_{a},$ (4)

where $R_{\infty}$ is the diffuse reflectance of the biological object. In the
near infrared wavelength region, ${\mathbb{R}}_{\infty}$ has been
experimentally found to be lying somewhere in the range of 0.5-1.25 (see
Figure 4 in [29]). Furthermore, a similar linear relation between $\sigma_{a}$
and $\sigma_{s}$ has also been found to be experimentally true for the tissues
of brain, breast, heart, lungs, bone, prostate and tumor; see, e.g., [26, 51],
and Table 1, where we report measured values of the absorption and reduced
scattering coefficients for various human tissues.

Table 1: Values of the absorption and the reduced scattering coefficients for various tissues at the near-infrared wavelengths [51]. Tissue | Wavelength (nm) |  $\sigma_{a}$ (cm-1) |  $\sigma_{s}$ (cm-1)
---|---|---|---
Brain | 760 | 0.11-0.17 | 4.0-10.5
Breast | 760 | 0.031-0.1 | 8.3-12.0
Heart | 661 | 0.12-0.18 | 5.22-90.80
Lungs | 661 | 0.49-0.88 | 21.14-22.52
Bone | 760 | 0.07-0.09 | 11.9-14.1
Prostate | 732 | 0.09-0.72 | 3.37-29.8
Tumor | 795 | 0.068-0.102 | 12.4-13.1

These results suggest that the reconstruction of the diffusion coefficient in
the framework of the diffusion approximation can be improved if we include our
qualitative prior knowledge represented by (3) and (4) in the formulation of
our inverse problem.

For this purpose, based on the results given in Table 1 and supported by the
Kubelka-Munk relation, we construct our prior based on the following tentative
relation between the absorption and the reduced scattering coefficients

$\sigma_{s}\approx c~{}\sigma_{a}.$ (5)

This choice is based on results of measurements indicating that the value of
$c$ varies among different biological tissues in the approximate range
$(30,800)$. Notice that in our approach we interpret (5) as a mean approximate
relation between $\sigma_{a}$ and $\sigma_{s}$ in the tissues (with the value
of $c$ arbitrarily chosen in a given interval determined by the application)
and aim at using this information for better reconstruction of the QPAT
coefficients.

We demonstrate that the consistent use of (3) and (5) embedded as a prior in
our nonsmooth optimization problem results in a novel reconstruction strategy
that allows to obtain reconstructions of $D$ and $\sigma_{a}$ with high
resolution and high contrast. However, it is important to emphasize that the
high-quality reconstructions that we obtain do not necessarily satisfy the
relations (3) and (5). On the other hand, we show that our computational
framework is certainly valid for biomedical applications with near-infrared
light illumination of tissues in order to accurately detect the cancerous ones
inside the organs mentioned in Table 1, where the Kubelka-Munk relation (5)
holds approximately true.

In our problem formulation, we consider the minimization of a Tikhonov-type
least squares functional with regularization terms including the relation
between $D$ and $\sigma_{a}$ mentioned above. Further, we include terms that
promote sparsity in the reconstruction, in order to accommodate the fact that,
in biomedical application, different tissues appear distributed as almost non-
overlapping bounded regions of the domain under consideration. However, while
the introduction of nonsmooth penalization terms in the functional
considerably contribute to the ability of our method to provide sharp
reconstructions as in [2, 49], they require the development of sophisticated
optimization techniques that apply in the nonsmooth case. We refer to [56] for
a work that pioneered the use of $L^{1}$ penalization terms to promote
sparsity in PDE optimization and the development of semi-smooth Newton (SSN)
methods to solve the corresponding problems. On the other hand, the
implementation of a SSN scheme requires the application of the Hessian of the
optimization problem and, for this purpose, in addition to solving the model
equation and its adjoint, one needs to solve other two linearized equations.
This fact and further implementation issues make the development of the SSN
scheme difficult to practitioners.

To overcome these difficulties, we extend the sequential quadratic hamiltonian
(SQH) algorithm, proposed in [18, 19, 20] for nonsmooth PDE optimal control
problems, to solve our inverse problem. This algorithm is based on the
Pontryagin maximum principle (PMP), which is well-known in the theory of
optimal control of dynamical systems. Our focus on the SQH algorithm is
motivated by the fact that it is easy to implement, it appears very efficient
and robust, and it is wellposed as an iterative scheme [18, 19]. A specific
feature of the iterative SQH method is that it relies on a pointwise
optimization procedure.

Several numerical experiments are performed to demonstrate the computational
effectiveness of our methodology. The results of these experiments show that
our QPAT inverse problem, which includes a prior relation on $D$ and
$\sigma_{a}$, and the SQH algorithm solving this problem, result in a fast and
accurate reconstruction procedure. These are the main novelties of our work.

This paper is organized as follows. In the next section, we discuss our
optimization formulation of the QPAT inverse problem, and illustrate our use
of a relation between $D$ and $\sigma$. Further, we discuss $L^{2}$ and
$L^{1}$ penalization terms. In Section 3, we provide theoretical consideration
on solutions to our minimization problem and discuss their characterization in
the PMP framework. Section 4 is devoted to illustrating all details of our
numerical optimization procedure that includes wellposedness of the SQH
algorithm and the numerical approximation of our PDE governing model and of
its optimization adjoint. In Section 5, results of numerical experiments are
presented that successfully validate our QPAT framework. For this purpose,
different test cases with synthetic data are considered that allow to
investigate the role of the optimization parameters and the ability of our
computational framework to provide reconstructions of the coefficients $D$ and
$\sigma$ with high contrast and resolution, also in the case of data affected
by noise. A section of conclusions completes this work.

##  2 Formulation of a QPAT optimization problem

We start discussing the governing model of our QPAT inverse problem. For this
purpose, we specify the dimension $n=2$, and refer to the set of coordinates
$(x,y)$. Let $u(x,y)$ be the total number of photons at a point $(x,y)$ inside
an object represented by the domain $\Omega$. In our discussion, $\Omega$ is
either a rectangle or a convex subset of ${\mathbb{R}}^{2}$ with boundary
$\partial\Omega$ of class $C^{2}$. We assume that $u$ satisfies the following
photon-propagation problem

$\displaystyle-\nabla\cdot(D(x,y)\,\nabla u(x,y))+\sigma_{a}(x,y)\,u(x,y)$
$\displaystyle=0,\qquad(x,y)\in\Omega$ (6) $\displaystyle u(x,y)$
$\displaystyle=g(x,y),\quad(x,y)\in\partial\Omega,$

where $D(x,y)>0$ represents the diffusion coefficient at $(x,y)$. The non-
negative absorption coefficient $\sigma_{a}$ is given by

$\sigma_{a}(x,y)=\sigma(x,y)+\sigma_{b},$

where $\sigma_{b}$ is a known background absorption coefficient and $\sigma$
represents the deviation from this value in correspondence of organs and sites
of tumor. We remark that even though $\sigma$ is a signed function, however
later on, we will impose bounds on $\sigma$ such that $\sigma_{a}>0$. Further,
the function $g$ represents the known intensity of radiation at the boundary
of the object. For the reconstruction of $D$ and $\sigma_{a}$, we require the
knowledge, by measurements, of the optical energy given by

$\mathcal{H}(x,y,\sigma,u_{i})=\Gamma\,(\sigma(x,y)+\sigma_{b})\,u_{i}(x,y),\qquad
i=1,2,$ (7)

for two sets of boundary radiation functions $g_{i}$, where $\Gamma$ is the
known Grüneisen coefficient that is assumed to be a known constant.

Our QPAT inverse problem consists in obtaining the unknown pair $(D,\sigma)$,
given $H$ and $g$. For this purpose, we formulate a minimization problem with
the following cost functional

$\displaystyle J(D,\sigma,u_{1},u_{2})$
$\displaystyle=\frac{\alpha}{2}\sum_{i=1}^{2}\int_{\Omega}(\mathcal{H}(x,y,\sigma,u_{i})-G_{i}^{\delta}(x,y))^{2}~{}dxdy+\frac{\xi_{1}}{2}\int_{\Omega}\sigma^{2}(x,y)~{}dxdy$
(8)
$\displaystyle+\frac{\xi_{2}}{2}\int_{\Omega}\left(D(x,y)-\frac{1}{3c~{}(\sigma(x,y)+\sigma_{b})}\right)^{2}~{}dxdy+\gamma\int_{\Omega}|\sigma(x,y)|~{}dxdy,$

where $G_{i}^{\delta}\in L^{\infty}(\Omega)$ is a given noisy measurement
function of the initial pressure distribution corresponding to the boundary
radiation functions $g_{i}$, $i=1,2$. Therefore the first term in the
objective functional $J$ represents a least-squares data fit of the initial
pressure distribution measurements. We assume that the background value of the
absorption coefficient $\sigma_{b}$ is known and therefore $\sigma$ equals
zero in the background. In this setting, it is reasonable to introduce a
$L^{2}$ and $L^{1}$ regularization of $\sigma$, which result in the second and
last terms in (8). In particular, the $L^{1}$ term promotes sparsity of
$\sigma$, thus enhancing contrast.

One novelty of our work is the prior given by the third term in (8). To
explain this term, we recall (3) and the related discussion. Hence, by
assuming $c\in[30,800]$, we can make the following approximation

$3\,(\sigma_{a}+\sigma_{s})\approx 3\,(\sigma_{a}+c\sigma_{a})\approx
3c\sigma_{a}.$

We remark that even though we will be choosing $c\in[30,800]$ in our numerical
experiments, the choice of $c$ is not essential. Thus, we employ the following
function as prior for the diffusion coefficient $D$

$D(x,y)\approx\bar{D}:=\frac{1}{3c~{}(\sigma(x,y)+\sigma_{b})}.$ (9)

Therefore the third term in $J$ measures the deviation of $D$ in the tissues
and tumors from $\bar{D}$.

Our QPAT optimization problem is stated as follows

$\displaystyle\min_{D,\sigma\in
D_{ad}\times\Sigma_{ad}}~{}J(D,\sigma,u_{1},u_{2})$ (10) s.t.
$\displaystyle\mathcal{L}(u_{1},D,\sigma,g_{1})=0,$
$\displaystyle\mathcal{L}(u_{2},D,\sigma,g_{2})=0,$

(‘s.t.’ means ‘subject to’) where $\mathcal{L}(u_{i},D,\sigma,g_{i})=0$,
$i=1,2$, denote (6) for the boundary radiation functions $g_{i}$, $i=1,2$,
respectively, and our aim is to recover $D$ and $\sigma$ in the following
admissible sets

$\displaystyle D_{ad}=\\{D\in L^{2}(\Omega):0<D_{l}\leq D(x,y)\leq
D_{r},\mbox{ a.e. in }\Omega\\},$ (11) $\displaystyle\Sigma_{ad}=\\{\sigma\in
L^{2}(\Omega):\sigma_{l}\leq\sigma(x,y)\leq\sigma_{r},\mbox{ a.e. in
}\Omega,\\},$

where $\sigma_{l}$ is chosen such that the total absorption coefficient
$\sigma_{a}=\sigma+\sigma_{b}>0$, e.g., one can choose
$\sigma_{l}=-\sigma_{b}+\sigma_{\epsilon}$, with $\sigma_{\epsilon}$ a small
positive number.

##  3 Analysis of the QPAT inverse problem

We remark that our diffusion model (6), with the given bounds on the
coefficients and Dirichlet boundary conditions, is a well-studied problem and
we refer to [35, 42] for two classical references. Specifically, one can
verify that with $D\in D_{ad}$ and $\sigma\in\Sigma_{ad}$ our operator is
uniformly elliptic and the boundedness conditions on the coefficients required
in, e.g., [35] are satisfied. Then, assuming that the Dirichlet data $g$
corresponds to the trace of a $H^{1}(\Omega)$ function, the following
proposition holds whose proof can be found e.g., in [42].

######  Proposition 1.

Let $D\in D_{ad}$, $\sigma\in\Sigma_{ad}$. Then there exists a unique solution
$u\in H^{1}(\Omega)$ to (6).

Additional regularity of the solution, that is, $u\in H^{1}(\Omega)\cap
H^{2}(\Omega)$, is attained at the interior of $\Omega$ if $D$ is assumed
Lipschitz continuous, and this regularity holds globally if $\partial\Omega$
is of class $C^{2}$; see [35]. Further results can be found in, e.g., [3]. In
particular, we refer to [22, 45] for the problem of establishing boundedness
of the gradient of solutions to the diffusion model.

With Proposition 1, we have that the map $(D,\sigma)\mapsto u=u(D,\sigma)$,
where $u$ is the solution to (6) for given $(D,\sigma)\in
D_{ad}\times\Sigma_{ad}\subset L^{\infty}(\Omega)\times L^{\infty}(\Omega)$,
is well defined. Furthermore, this map is continuous as shown below; see also
[40].

######  Lemma 3.1.

Let $\Omega\subset{\mathbb{R}}^{2}$ and consider any two pairs
$(u_{1},D_{1},\sigma_{1}),\,(u_{2},D_{2},\sigma_{2})\in H^{1}(\Omega)\times
D_{ad}\times\Sigma_{ad}$ that satisfy (6) with $g$ the trace of a
$H^{1}(\Omega)$ function. Then the following estimate holds

$\|u_{1}-u_{2}\|_{H^{1}(\Omega)}\leq
C_{1}\,\big{(}\|D_{1}-D_{2}\|_{L^{\infty}(\Omega)}+\|\sigma_{1}-\sigma_{2}\|_{L^{\infty}(\Omega)}\big{)}.$

Furthermore, if for the given pairs the gradient of $u_{1}$ (respectively
$u_{2}$), is essentially bounded, then the following estimate holds

$\|u_{1}-u_{2}\|_{H^{1}(\Omega)}\leq
C_{2}\,\big{(}\|D_{1}-D_{2}\|_{L^{2}(\Omega)}+\|\sigma_{1}-\sigma_{2}\|_{L^{2}(\Omega)}\big{)}.$

###### Proof.

The two pairs satisfy both the diffusion equation in weak form. We have

$\left(D_{i}\,\nabla u_{i},\nabla
v\right)+\left(\tilde{\sigma}_{i}\,u_{i},v\right)=0,\qquad v\in
H_{0}^{1}(\Omega),$

where $\tilde{\sigma}_{i}=\sigma_{i}+\sigma_{b}$ and $u_{i}=g$ on the
boundary, $i=1,2$.

Therefore we have

$\left(D_{1}\,\nabla u_{1},\nabla
v\right)+\left(\tilde{\sigma}_{1}\,u_{1},v\right)=\left(D_{2}\,\nabla
u_{2},\nabla v\right)+\left(\tilde{\sigma}_{2}\,u_{2},v\right).$

Thus, by algebraic manipulation, we obtain

$\left(D_{1}\,\left[\nabla u_{1}-\nabla u_{2}\right],\nabla
v\right)+\left(\tilde{\sigma}_{1}\,\left[u_{1}-u_{2}\right],v\right)=\left(\left[D_{2}-D_{1}\right]\,\nabla
u_{2},\nabla
v\right)+\left(\left[\sigma_{2}-\sigma_{1}\right]\,u_{2},v\right).$

In particular, taking $v=u_{1}-u_{2}$, we have

$\displaystyle\left(D_{1}\,\nabla\left[u_{1}-u_{2}\right],\nabla\left[u_{1}-u_{2}\right]\right)+\left(\tilde{\sigma}_{1}\,\left[u_{1}-u_{2}\right],\left[u_{1}-u_{2}\right]\right)$
$\displaystyle=\left(\left[D_{2}-D_{1}\right]\,\nabla
u_{2},\nabla\left[u_{1}-u_{2}\right]\right)+\left(\left[\sigma_{2}-\sigma_{1}\right]\,u_{2},\left[u_{1}-u_{2}\right]\right).$

Since $D_{i}\in D_{ad}$ and $\sigma_{i}\in\Sigma_{ad}$, we have

$c_{1}\,\|u_{1}-u_{2}\|_{H^{1}(\Omega)}^{2}\leq\left(\left[D_{2}-D_{1}\right]\,\nabla
u_{2},\nabla\left[u_{1}-u_{2}\right]\right)+\left(\left[\sigma_{2}-\sigma_{1}\right]\,u_{2},\left[u_{1}-u_{2}\right]\right).$

where $c_{1}=\min\\{D_{l},\sigma_{l}+\sigma_{b}\\}$. Hence, considering that
both $D$ and $\sigma$ are bounded, we have the estimate

$c_{1}\,\|u_{1}-u_{2}\|_{H^{1}(\Omega)}\leq
c_{2}\,\big{(}\|D_{1}-D_{2}\|_{L^{\infty}(\Omega)}+\|\sigma_{1}-\sigma_{2}\|_{L^{\infty}(\Omega)}\big{)},$

where $c_{2}=\max\\{\|u_{2}\|_{L^{2}(\Omega)},\|\nabla
u_{2}\|_{L^{2}(\Omega)}\\}$. Hence, the first statement of the lemma is proved
with $C_{1}=c_{2}/c_{1}$.

Further, by the assumption that $\|\nabla u_{2}\|_{L^{\infty}(\Omega)}\leq M$
(resp. $\|\nabla u_{1}\|_{L^{\infty}(\Omega)}\leq M$), for some positive
constant $M$, we have

$\displaystyle c_{1}\,\|u_{1}-u_{2}\|_{H^{1}(\Omega)}^{2}$
$\displaystyle\leq\|\nabla
u_{2}\|_{L^{\infty}(\Omega)}\,\|D_{2}-D_{1}\|_{L^{2}(\Omega)}\,\|\nabla\left[u_{1}-u_{2}\right]\|_{L^{2}(\Omega)}$
$\displaystyle+\|u_{2}\|_{L^{\infty}(\Omega)}\,\|\sigma_{2}-\sigma_{1}\|_{L^{2}(\Omega)}\,\|u_{1}-u_{2}\|_{L^{2}(\Omega)}$
$\displaystyle\leq\big{(}\|\nabla
u_{2}\|_{L^{\infty}(\Omega)}\,\|D_{2}-D_{1}\|_{L^{2}(\Omega)}+\|u_{2}\|_{L^{\infty}(\Omega)}\,\|\sigma_{2}-\sigma_{1}\|_{L^{2}(\Omega)}\big{)}\,\|u_{1}-u_{2}\|_{H^{1}(\Omega)}$

Therefore the claim of the lemma follows from the following result

$c_{1}\,\|u_{1}-u_{2}\|_{H^{1}(\Omega)}\leq
c_{3}\,\big{(}\|D_{1}-D_{2}\|_{L^{2}(\Omega)}+\|\sigma_{1}-\sigma_{2}\|_{L^{2}(\Omega)}\big{)},$

where $c_{3}=\max\\{\|u_{2}\|_{L^{\infty}(\Omega)},\|\nabla
u_{2}\|_{L^{\infty}(\Omega)}\\}$. Thus, the second statement is proved with
$C_{2}=c_{3}/c_{1}$. ∎

With this preparation, we can introduce the following reduced cost functional

$\hat{J}(D,\sigma)=J(D,\sigma,u_{1}(D,\sigma),u_{2}(D,\sigma)).$ (12)

Notice that, with this functional, our QPAT optimization problem can be
equivalently formulated as the following unconstrained minimization problem

$\min_{(D,\sigma)\in D_{ad}\times\Sigma_{ad}}\hat{J}(D,\sigma).$ (13)

The proof of existence of a minimizer for (13) by using the classical direct
method is a delicate issue since for the minimizing sequences
$\\{D^{k}\\}\subset L^{\infty}(\Omega)$ we do not have the strong convergence
required to prove fulfilment of the differential constraint (6) in the limit
of the minimizing process; see, e.g., Chapter VI of [16] for a related
framework. However, this would be possible assuming a compact set of
admissible controls, e.g., $\tilde{D}_{ad}=D_{ad}\cap V_{l}$ where
$V_{l}:=\\{v\in H^{1}(\Omega)\,|\,\|v\|_{H^{1}\left(\Omega\right)}\leq l\\}$
where $l$ is a given positive constant; see, e.g., [44] for a discussion on
this setting. Nevertheless, since these are sufficient and not necessary
conditions for existence of an optimal pair $(D^{*},\sigma^{*})$, we prefer to
work with the admissible sets defined in (11) that appear convenient in
applications, and assume existence of the optimal optical coefficients
$(D^{*},\sigma^{*})$. Our focus is the construction of a method that
efficiently computes these coefficients.

Now, concerning the characterization of solutions to our QPAT optimization
problem, we can use subdifferential calculus to derive first-order necessary
optimality conditions [23, 60]. In particular, if $\hat{J}(D,\sigma)$ is
Fréchet differentiable, then the optimal pair $(D^{*},\sigma^{*})$ must
satisfy

$\left\langle\nabla\hat{J}(D^{*},\sigma^{*}),(D,\sigma)-(D^{*},\sigma^{*})\right\rangle\geq
0,\qquad(D,\sigma)\in D_{ad}\times\Sigma_{ad},$

where $\langle\cdot,\cdot\rangle$ denotes the duality product in $L^{2}$.

Alternatively, we can use the Pontryagin maximum principle that can be proven
in our case provided that the gradient of the adjoint variable is essentially
bounded. However, we can also apply a linearized version of the PMP as
discussed in [58]. In both cases, the starting point of the PMP approach is
the construction of the Hamiltonian-Pontryagin (HP) function
$H:\mathbb{R}^{2}\times[D_{l},D_{r}]\times[\sigma_{l},\sigma_{r}]\times\mathbb{R}^{4}$
that is defined as follows (see [18, 19, 20])

$\displaystyle
H(x,y,D(x,y),\sigma(x),u_{1}(x,y),u_{2}(x,y),q_{1}(x,y),q_{2}(x,y))$ (14)
$\displaystyle=\frac{\alpha}{2}\sum_{i=1}^{2}(\mathcal{H}(x,y,D(x,y),\sigma(x,y),u_{i}(x,y))-G_{i}^{\delta}(x,y))^{2}+\frac{\xi_{1}}{2}\sigma^{2}(x,y)$
$\displaystyle+\frac{\xi_{2}}{2}\left(D(x,y)-\bar{D}(x,y)\right)^{2}+\gamma|\sigma(x)|$
$\displaystyle+\sum_{i=1}^{2}D(x,y)\nabla u_{i}(x,y)\cdot\nabla
q_{i}(x,y)+\sum_{i=1}^{2}(\sigma(x,y)+\sigma_{b})\,u_{i}(x,y)\,q_{i}(x,y),$

where $q_{i},~{}i=1,2$ solve the following adjoint equations, which are
derived from the Lagrangian and Hamiltonian framework [4]

$\displaystyle-\nabla\cdot(D(x,y)\,\nabla
q_{i}(x,y)))+(\sigma(x,y)+\sigma_{b})\,q_{i}(x,y)$
$\displaystyle=-\alpha\Gamma~{}(\sigma(x,y)+\sigma_{b})\,[\mathcal{H}(x,y,\sigma,u_{i})-G_{i}^{\delta}],~{}\mbox{
in }\Omega$ (15) $\displaystyle q_{i}(x,y)$ $\displaystyle=0,\qquad\mbox{ on
}\partial\Omega,$

Notice that similar consideration, as for the diffusion model, concerning
existence and regularity of solutions apply to these elliptic problems.

Then subject to appropriate conditions the PMP characterization of a minimum
of the minimization problem (10) is given by the following theorem (expressed
in the equivalent form of a minimum principle)

######  Theorem 3.1.

Let $(D^{*},\sigma^{*},u_{1}^{*},u_{2}^{*})$ be a solution to (10). Then

$\displaystyle
H(x,y,D^{*}(x,y),\sigma^{*}(x,y),u_{1}^{*}(x,y),u_{2}^{*}(x,y),q_{1}^{*}(x,y),q_{2}^{*}(x,y))$
(16)
$\displaystyle=\mathop{\min_{(v_{1},v_{2})~{}\in[D_{l},D_{r}]\times[\sigma_{l},\sigma_{r}]}}H(x,y,v_{1},v_{2},u_{1}^{*}(x,y),u_{2}^{*}(x,y),q_{1}^{*}(x,y),q_{2}^{*}(x,y)),$

for almost all $(x,y)\in\Omega$, and $u_{i}^{*},~{}i=1,2$ and
$q_{i}^{*},~{}i=1,2$ solve (6) and (15), respectively, with $D=D^{*}$ and
$\sigma=\sigma^{*}$.

We can rewrite condition (16) as follows

$(D^{*}(x,y),\sigma^{*}(x,y))=\mathop{\operatorname*{arg\,min}_{(v_{1},v_{2})~{}\in[D_{l},D_{r}]\times[\sigma_{l},\sigma_{r}]}}H(x,y,v_{1},v_{2},u_{1}^{*}(x,y),u_{2}^{*}(x,y),q_{1}^{*}(x,y),q_{2}^{*}(x,y)).$
(17)

We note that (17) represents a pointwise optimality criteria at each $(x,y)$
and does not involve any derivatives with respect to the optimization
variables $D$ and $\sigma$. In the next section, we discuss an efficient and
robust numerical scheme based on (17) that solves (10) (resp. (13)) to obtain
$D^{*}$ and $\sigma^{*}$.

##  4 The SQH algorithm for solving the QPAT inverse problem

The sequential quadratic hamiltonian (SQH) algorithm has been proposed in [18,
20], as a new variant of the successive approximations methods; see [21] for
an early review. It includes an adaptive quadratic penalty of the updates that
results in an augmented Hamiltonian as first proposed in [50]. In this way, an
efficient and robust iterative procedure is obtained that is able to solve
nonsmooth optimization problems governed by PDEs.

In our SQH implementation, the augmented HP function is given by

$H_{\epsilon}(x,y,D,\sigma,\tilde{D},\tilde{\sigma},u_{1},u_{2},q_{1},q_{2})=H(x,y,D,\sigma,u_{1},u_{2},q_{1},q_{2})+\epsilon\,[(D-\tilde{D})^{2}+(\sigma-\tilde{\sigma})^{2}],$
(18)

where $\epsilon>0$ is a penalization parameter that is chosen adaptively at
each step of the SQH iteration process: a larger value of $\epsilon$ is chosen
if a sufficient decrease of the functional $J$ is not observed while a smaller
value of $\epsilon$ is chosen if $J$ decreases sufficiently. Specifically, if
$\tilde{D}$ and $\tilde{\sigma}$ denote a previous approximation to the
diffusion and absorption coefficients sought, then the purpose of quadratic
term $\epsilon\,[(D-\tilde{D})^{2}+(\sigma-\tilde{\sigma})^{2}]$ is to have
the pointwise minimizer of $H_{\epsilon}$, and thus an update to $D$ and
$\sigma$, that is close to the previous values $\tilde{D}$ and
$\tilde{\sigma}$ as much as $\epsilon$ is large. We also remark that during a
step of the optimization problem, that is, for a minimization sweep on all
grid points of the $(x,y)$ mesh and, during this process, the values of $u$
and $q$ are those obtained in the previous iteration. The SQH algorithm is
given as a pseudocode as follows.

######  Algorithm 4.1 (SQH algorithm).

  1. 1.

Choose $\epsilon>0$, $\kappa>0$, $\lambda>1$, $\zeta\in(0,1)$,
$\rho\in(0,\infty)$, $(D^{0},\sigma^{0})\in D_{ad}\times\Sigma_{ad}$

  2. 2.

Compute $u_{i}^{0}(D^{0},\sigma^{0})$ and
$q_{i}^{0}(D^{0},\sigma^{0},u_{i}^{0})$ for $i=1,2$. Set $k=0$.

  3. 3.

Find $(D,\sigma)\in D_{ad}\times\Sigma_{ad}$ such that

$H_{\epsilon}(x,y,D,\sigma,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})=\mathop{\min_{w\in[D_{l},D_{r}]}}_{z\in[\sigma_{l},\sigma_{r}]}H_{\epsilon}(x,y,w,z,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})$

for all $(x,y)\in\Omega.$

  4. 4.

Calculate $u_{i}(D,\sigma)~{},i=1,2$, and set
$\tau=\|D-D^{k}\|^{2}_{L^{2}(\Omega)}+\|\sigma-\sigma^{k}\|^{2}_{L^{2}(\Omega)}.$

  5. 5.

If
$J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})>-\rho\,\tau$,
choose $\epsilon=\lambda\,\epsilon$,
else choose $\epsilon=\zeta\,\epsilon$, $u_{i}^{k+1}=u_{i}$, $D^{k+1}=D$,
$\sigma^{k+1}=\sigma$, and calculate the adjoint variables $q_{i}^{k+1}$,
$i=1,2$, corresponding to $D^{k+1}$, $\sigma^{k+1}$, $u_{i}^{k+1}$, $i=1,2$.

  6. 6.

Set $k=k+1$

  7. 7.

If $\tau<\kappa$, STOP and return optimal $(D^{k},\sigma^{k})$ else go to Step
3.

  8. 8.

end

The pointwise minimization in Step 3 of this algorithm is carried out by a
direct search process, i.e., we divide the domain
$[D_{l},D_{r}]\times[\sigma_{l},\sigma_{r}]$ into a two-dimensional discrete
grid, evaluate $H_{\epsilon}$ on this grid, and then choose a grid point
$(w,z)$ at which $H_{\epsilon}$ is minimum. We remark that in Step 5 of this
algorithm, if the inequality
$J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})>-\rho\,\tau$
is true this means that no sufficient decrease of the value of the objective
functional $J$ as been achieved. In this case, a larger value of $\epsilon$ is
taken (since $\lambda>1$) and the optimization procedure in Step 3. is
repeated with the corresponding new augmented HP function. On the other hand,
if the inequality above is false, the decrease has met the required criteria
for reduction of the value of $J$, and the new values $D^{k+1}=D$,
$\sigma^{k+1}=\sigma$ are obtained. Correspondingly, the updates $u_{i}^{k+1}$
and $q_{i}^{k+1}$, $i=1,2$, are computed. In this case, the value of
$\epsilon$ is reduced by a factor $\zeta<1$.

The fundamental result required to guarantee wellposedness of the SQH
algorithm applied to our QPAT inverse problem is to prove that, if not already
at the optimal point, it is always possible to find a value of $\epsilon>0$
such that the minimization of the corresponding $H_{\epsilon}$ results in new
values of $D$ and $\sigma$ for which the value of the functional $J$ is
reduced. We present this result in Theorem 4.1 below and refer to [18, 19] for
the use of this theorem for further results on the convergence of the SQH
algorithm in a finite number of steps.

######  Theorem 4.1.

Let $(u_{1},u_{2},D,\sigma)$ and $(u_{1}^{k},u_{2}^{k},D^{k},\sigma^{k})$ be
as in Algorithm 4.1, $k\in\mathbb{N}_{0}$, and denote $\delta
D:=D-D^{k},~{}\delta\sigma:=\sigma-\sigma^{k}$, $\delta
u_{i}:=u_{i}-u_{i}^{k}$, $i=1,2$. If the assumptions of Lemma 3.1 are
satisfied, then there is a $\theta>0$, independent of $\epsilon$, $k$, and
$D^{k}$, $\sigma^{k}$, such that for the $\epsilon$ currently chosen by the
SQH Algorithm 4.1, it holds that

$J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})<-\left(\epsilon-\theta\right)\,\big{(}\|D-D^{k}\|_{L^{2}(\Omega)}^{2}+\|\sigma-\sigma^{k}\|_{L^{2}(\Omega)}^{2}\big{)}.$

###### Proof.

See the Appendix. ∎

For implementing Algorithm 4.1, we need to solve the photon propagation
equation (6) for different boundary conditions $g_{i},~{}i=1,2$, and the
corresponding adjoint equations (15). We use a cell-nodal finite-difference
scheme, proposed in [7], to discretize (6) and (15). To illustrate this
scheme, we consider the generic form of the equations (6) and (15) defined in
the domain $\Omega=(a,b)^{2}\subset\mathbb{R}^{2}$ as follows

$\displaystyle-\nabla\cdot({D(x,y)}\nabla{u(x,y)})+(\sigma(x,y)+\sigma_{b})\,u(x,y)$
$\displaystyle=f$ $\displaystyle\qquad\mbox{ in }\Omega,$ (19) $\displaystyle
u(x,y)$ $\displaystyle=f_{D}(x,y)$ $\displaystyle\qquad\mbox{ on
}\partial\Omega,$

where $f$ and $f_{D}$ are assumed continuous in $\Omega$ and on
$\partial\Omega$, respectively. Consider a sequence of uniform grids
$\\{\Omega_{h}\\}_{h>0}$ given by

$\Omega_{h}=\\{{(x_{1}^{i},x_{2}^{j})\in\mathbb{R}^{2}:(x_{1}^{i},x_{2}^{j})=(a+ih,a+jh),~{}(i,j)\in\\{{0,\ldots,N}\\}^{2}}\\}\cap\Omega,$
(20)

where $N$ represents the number of cells in each direction and
$h=\dfrac{(b-a)}{N}$ is the mesh size. Then the corresponding cell-nodal
discretization places the unknowns $u,D,\sigma$ at the nodal points
$(x^{i},y^{j})$ of the grid that results in the following approximation of
(19) at $(x^{i},y^{j})$

$\displaystyle\dfrac{1}{h^{2}}\Bigg{\\{}({D_{i+1/2,j}}+{D_{i-1/2,j}}+{D_{i,j+1/2}}+{D_{i,j-1/2}})u_{i,j}$
(21)
$\displaystyle-{D_{i+1/2,j}}u_{i+1,j}-{D_{i-1/2,j}}u_{i-1,j}-{D_{i,j+1/2}}u_{i,j+1}-{D_{i,j-1/2}}u_{i,j-1}\Bigg{\\}}$
$\displaystyle+(\sigma_{i,j}+\sigma_{b})u_{i,j}=f_{i,j},\qquad 1\leq i,j\leq
N-1.$

We denote $D_{i\pm 1,j}=D(x^{i}\pm h,y^{j})$, $D_{i,j\pm 1}=D(x^{j},y^{j}\pm
h)$, and $f_{i,j}=f(x^{i},y^{j})$. The required intermediate values of $D$ are
computed as follows

${D_{i\pm 1/2,j}}=\dfrac{1}{2}\Bigg{(}{D_{i\pm
1,j}}+{D_{i,j}}\Bigg{)}\quad\mbox{ and }\quad{D_{i,j\pm
1/2}}=\dfrac{1}{2}\Bigg{(}{D_{i,j\pm 1}}+{D_{i,j}}\Bigg{)}.$

The Dirichlet boundary data $f_{D}$ is included in the standard way in the
right-hand side of the algebraic equation. We use one-sided finite difference
discretizations to evaluate the right-hand sides of the adjoint equations.
With this discretization scheme in (21), we solve a matrix-vector equation to
obtain the unknowns $u_{i,j}$.

##  5 Numerical results

We present results of numerical experiments to validate the ability of our
QPAT optimization framework in the reconstruction of the diffusion coefficient
$D$ and the absorption coefficient $\sigma$ profiles with high accuracy. These
also demonstrate the efficiency and robustness of our SQH-QPAT methodology. We
consider a widely varying set of phantoms for $D,\sigma$. In Test Case 1, we
consider a disk phantom with different intensities for $D,\sigma$ and
investigate reconstruction with the different regularization terms considered
in our cost functional $J$. Next, we discuss two experiments with the heart
lung phantom: one in which the relation (9) is satisfied between ${D}$ and
$\sigma$, and the other in which (9) is not satisfied, making it a very
generic set of phantoms for $D,\sigma$. We demonstrate that in both cases, by
incorporating the aprior assumption (9) in the functional $J$, we can obtain
high quality reconstructions. We conclude our set of experiments with the more
challenging set of two Shepp-Logan phantoms: the standard one, and another
augmented with tumor structures. In the test cases with the heart and lung
phantom and the Shepp-Logan phantom, we also implement our method on data
containing additive Gaussian noise in order to demonstrate the robustness of
our framework.

For the experiments below, we first remark that the units of the spatial
variable $x$ is $cm$ and of the absorption coefficient is $cm^{-1}$. Starting
with an experimental square domain of side length 5 cm, we first scale and
translate the diffusion equation (6), such that the domain
$\Omega=(-1,1)\times(-1,1)$, is discretized into 100 equally spaced points in
both the coordinate directions. The corresponding scaling is 1 experimental
unit corresponds to 2.5 cm. Thus, based on Table 1, an experimental range of
0-1 for $\sigma_{a}$ corresponds to 0-0.4 $cm^{-1}$. The two boundary
radiation functions are $g_{1}(x,y)=e^{x},~{}g_{2}(x,y)=e^{y}.$ The weights of
the functional $J$ given in (8) are chosen as
$\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0.01$ and $\Gamma=1.0$. The large
value of $\xi_{2}$ is chosen to emphasize the need to satisfy the relation
between the diffusion and the absorption coefficient, given by (9), in the
numerical experiments. The value of $c$ in the Kubelka-Munk relation (9) is
chosen to be 100/3, though as we have mentioned earlier, this choice is not
essential. The initial value of $\epsilon=10$ and the value of the stopping
criteria for convergence is chosen as $\kappa=10^{-6}$. To generate the data
$G_{i}^{\delta},~{}i=1,2$, we first solve for $u_{i}$ in (6) with given test
values of $D,\sigma$ and boundary illumination data $g_{i}$ on a finer mesh
with $N=400$ using the cell-nodal finite difference model described in Section
4. We next compute $G_{i}^{\delta}$ on the finer mesh using the values of
$\sigma,\sigma_{b},u_{i}$ from (7). We finally restrict the resulting
$G_{i}^{\delta}$ onto the coarser mesh with $N=100$ and take this as our given
data.

To quantify the comparisons of the reconstructions obtained in terms of their
contrast (function values) and resolution (discontinuities or edges), with and
without noise in the data, we use the following quantitative figures of merit
[39]

$\displaystyle\mbox{Relative mean square error (RMSE)
percentage}=\frac{\|f_{\mbox{rec}}-f_{\mbox{ex}}\|_{2}}{\|f_{\mbox{ex}}\|_{2}}*100\%,$
(22) $\displaystyle\mbox{Peak signal-to-noise ratio
(PSNR)}=10\log_{10}\left(\frac{\max{f_{\mbox{ex}}}}{\|f_{\mbox{rec}}-f_{\mbox{ex}}\|^{2}_{2}}\right),$

where $f_{ex}$ is the exact phantom and $f_{rec}$ is the reconstructed
phantom.

Test Case 1: we consider a phantom represented by a disk centered at
$(0.25,0.25)$ and having radius 0.25. The value of $\sigma$ inside the disk is
1 and outside is 0. The background value $\sigma_{b}$ is chosen to be 0.16.
The corresponding value of D is chosen as 0.003 inside the disk and 0.02
outside. The plots of the actual phantoms for $D$ and the total absorption
coefficient $\sigma_{a}$ are shown in Figure 1.

(a) Exact $D$

(b) Exact $\sigma_{a}$

Figure 1: Test Case 1: Exact $D$ and $\sigma$ represented by disk phantoms.

In Figure 2, we present the reconstruction of $D$ and $\sigma_{a}$ using
various choices of the regularization terms in $J$. The first column of images
depicts reconstructions of $D$ and $\sigma_{a}$ by setting the weight
$\xi_{2}=0$ of the regularization term $\left\|D-\bar{D}\right\|_{2}^{2}$,
which is the term that measures the deviation of $D$ from the Kubelka-Munk
approximation $\bar{D}$, given in (9). We see that the reconstruction of $D$
has a very poor resolution and contrast and in turn, this causes loss of
contrast in the reconstruction of $\sigma_{a}$. In the second column, we use
our regularization term $\left\|D-\bar{D}\right\|_{2}^{2}$ but setting the
weight $\gamma=0$ of the sparsity that promotes the $L^{1}$ regularization of
$\sigma$. We now see a reconstructed $\sigma_{a}$ with a poor contrast with
respect to the background $\sigma_{b}$. The reconstructed $D$ is now far more
superior than in the previous cases, but there is still a significant loss of
contrast.

(a) Reconstruction of $D$ with $\alpha=1,\xi_{1}=0.01,\xi_{2}=0,\gamma=0.01$

(b) Reconstruction of $D$ with $\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0$

(c) Reconstruction of $\sigma_{a}$ with
$\alpha=1,\xi_{1}=0.01,\xi_{2}=0,\gamma=0.01$

(d) Reconstruction of $\sigma_{a}$ with
$\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0$

Figure 2: Test Case 1: The reconstructed disk phantoms with different choices
of the values of the regularization weights.

Next, we show that in order to enhance contrast and resolution in the
reconstruction, we need to have the $L^{2}$ regularization term
$\left\|D-\bar{D}\right\|_{2}^{2}$ for $D$ and the sparsity promoting $L^{1}$
regularization term for $\sigma$. In this setting, Figure 3 shows the
reconstructions of $D$ and $\sigma_{a}$ with our framework using the SQH
algorithm. We see that the reconstructed functions are of high resolution and
high contrast when compared to the actual phantoms in Figure 1. We also
compute the RMSE percentages and PSNR values for the reconstructions in Table
2. We observe that without the Kubelka-Munk regularization term for $D$ and
without the $L^{1}$ regularization term for $\sigma$, we have very high values
of RMSE % and low values of PSNR values. With those regularization terms,
using the SQH algorithm, we obtain RMSE % of around 1 and high PSNR values,
reaffirming the superiority of the reconstructions.

(a) Reconstruction of $D$ with $\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0.01$

(b) Reconstruction of $\sigma_{a}$ with
$\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0.01$

Figure 3: Test Case 1: The reconstructed disk phantoms using the weights $\alpha=1,\xi_{1}=0.01,\xi_{2}=20,\gamma=0.01$.  Method | RMSE % ($D$) | RMSE % ($\sigma_{a}$) | PSNR ($D$) | PSNR ($\sigma_{a}$)
---|---|---|---|---
Without Kubelka-Munk ($\xi_{2}=0$) | 225.7 | 9.24 | 25.09 | 32.37
Without $L^{1}$ ($\gamma=0$) | 218 | 30.47 | 27.06 | 23.51
60mmWith Kubelka-Munk |  |  |  |
$(\xi_{2}=20)$ and $L^{1}$ ($\gamma=0.01)$ | 0.08 | 1.41 | 94.33 | 50.29
Table 2: Test Case 1: RMSE % and PSNR values for reconstructions of the disk
phantom with different regularization weights

Test Case 2: we demonstrate the robustness of our framework for reconstructing
multiple objects. In this test case, we consider a heart lung phantom for
$\sigma$. The underlying value of the phantom is 0 that is perturbed into two
ellipses that represent the lungs with value 1, and a disk representing the
heart with value 0.5. The value of $\sigma_{b}$ is chosen to be 0.03. The
value of $D$ is chosen as 0.003 inside the ellipses, 0.006 inside the disk and
0.1 elsewhere. The plots of the exact and the reconstructed phantoms are shown
in Figure 4.

(a) Exact $D$

(b) Exact $\sigma_{a}$

(c) Reconstructed $D$

(d) Reconstructed $D$ with 5% noise

(e) Reconstructed $D$ with 10% noise

(f) Reconstructed $\sigma_{a}$

(g) Reconstructed $\sigma_{a}$ with 5% noise

(h) Reconstructed $\sigma_{a}$ with 10% noise

Figure 4: Test Case 2: The actual and the reconstructed phantoms

We see from the reconstructed images of $D$ and $\sigma_{a}$ in Figures 4c and
4f that they are of high contrast and resolution. Now, in order to test the
robustness of our method with noisy data, we introduce 5% and 10%
multiplicative Gaussian noise to the interior data $\mathcal{H}$. We also
modify the $L^{2}$ and $L^{1}$ regularization weights for $\sigma$,
$\xi_{1}=0.1$ and $\gamma=0.1$. This ensures removal of artifacts from the
reconstruction of $\sigma$ due to the noisy data, which in turn helps to
remove some artifacts in the reconstruction of $D$. Figures 4d, 4e, 4g and 4h
show the corresponding reconstructions. We see that the reconstructions of $D$
and $\sigma_{a}$ involve a few artifacts near the bottom, more pronounced in
the case of 10% noise. But overall, both reconstructions still preserve obtain
high contrast and resolution near and at the regions of the heart and lungs.
This is also validated analytically from the low RMSE % and high PSNR values
of the reconstructions presented in Table 3.

Noise % | RMSE % ($D$) | RMSE % ($\sigma_{a}$) | PSNR ($D$) | PSNR ($\sigma_{a}$)
---|---|---|---|---
0 | 2.32 | 3.42 | 59.96 | 34.55
5 | 2.94 | 3.46 | 57.08 | 33.89
10 | 2.86 | 7.56 | 56.24 | 28.05
Table 3: Test Case 2: RMSE % and PSNR values for reconstructions of the heart
and lung phantom with different noise levels in the data

Test Case 3: we investigate the robustness of our SQH-QPAT algorithm in
reconstructing $D$ in the case where the Kubelka-Munk relation (9) is by far
not satisfied; i.e. we have factors different from 100. For this test case, in
regard to $\sigma$, we again consider the setting for heart and lung phantom
as described in test Case 2. However, for the value of $D$, we choose
$D=0.0006$ inside the right ellipse, $D=0.003$ inside the left ellipse,
$D=0.009$ inside the disk and $D=0.1$ for the background. The value of
$\sigma_{b}$ is again chosen to be 0.03.

(a) Exact $D$

(b) Exact $\sigma_{a}$

(c) Reconstructed $D$

(d) Reconstructed $D$ with 5% noise

(e) Reconstructed $D$ with 10% noise

(f) Reconstructed $\sigma_{a}$

(g) Reconstructed $\sigma_{a}$ with 5% noise

(h) Reconstructed $\sigma_{a}$ with 10% noise

Figure 5: Test Case 3: The actual and the reconstructed phantoms

With this challenging setting, we perform our reconstructions of $D$ and
$\sigma_{a}$ and depict the results in Figure 5 with and without noise in the
data. We again observe high contrast and resolution reconstructions of $D$ and
$\sigma_{a}$, even in the presence of 5% and 10% noise in the data. The
corresponding RMSE % and PSNR values are presented in Table 4, which are
similar to the previous test case.

Noise % | RMSE % ($D$) | RMSE % ($\sigma_{a}$) | PSNR ($D$) | PSNR ($\sigma_{a}$)
---|---|---|---|---
0 | 1.47 | 3.25 | 56.69 | 35.88
5 | 1.43 | 3.57 | 54.18 | 34.57
10 | 2.13 | 4.61 | 47.70 | 30.05
Table 4: Test Case 3: RMSE % and PSNR values for reconstructions of the second
heart and lung phantom with different noise levels in the data

Test Case 4: we perform another challenging test using the Shepp-Logan phantom
that represents the structure of a human head, which comprises ellipses, disks
of varying radii, and an elliptical annulus. In this case, $\sigma$ has the
structure of Shepp-Logan phantom as given in [55]. The phantom is described as
follows: The values of $\sigma$ and $D$, respectively, inside two big ellipses
centered at (0.22,0) and (-0.22,0) are -0.3 and 0.02, inside the big disk
centered at (0,0.35) are 0.5 and 0.003, inside the 2 small disks centered at
(0,0.1) and (0,-0.1) are 0.5 and 0.003, inside the small ellipses centered at
(-0.08,-0.605), (0,-0.606) and (0.06, -0.605) are 0.5 and 0.003, inside the
elliptical annulus are 1 and 0.002, and elsewhere 0 and 0.006. The background
$\sigma_{b}$ is chosen to be 0.5 in this case.

(a) Exact $D$

(b) Exact $\sigma_{a}$

(c) Reconstructed $D$

(d) Reconstructed $D$ with 5% noise

(e) Reconstructed $D$ with 10% noise

(f) Reconstructed $\sigma_{a}$

(g) Reconstructed $\sigma_{a}$ with 5% noise

(h) Reconstructed $\sigma_{a}$ with 10% noise

Figure 6: Test Case 4: The actual and the reconstructed phantoms

For test Case 4, Figure 6 shows the exact and the reconstructed values of $D$
and $\sigma_{a}$ with high resolution and contrast. Figures 6d, 6e and 6g, 6h
depict the reconstruction with 5% and 10% multiplicative Gaussian noise in the
interior data, respectively. We obtain high quality reconstruction with our
SQH-QPAT framework also in the case of objects with holes and inclusions. The
corresponding RMSE % and PSNR values are presented in Table 5. We do observe
that the RMSE % is higher than in the case of the heart and lung phantom,
which is expected since there are far more finer features in the Shepp-Logan
phantom. However, the PSNR values are still at a comparable level to the
previous test cases that indicate that the reconstructions have high quality.

Noise % | RMSE % ($D$) | RMSE % ($\sigma_{a}$) | PSNR ($D$) | PSNR ($\sigma_{a}$)
---|---|---|---|---
0 | 16.21 | 6.54 | 55.89 | 38.57
5 | 16.36 | 6.68 | 55.80 | 38.02
10 | 16.58 | 6.71 | 55.47 | 36.31
Table 5: Test Case 4: RMSE % and PSNR values for reconstructions of the Shepp-
Logan phantom with different noise levels in the data

Test Case 5: while in the previous test cases the phantoms to be reconstructed
represented organs in a human being, in this test case, we present results of
detecting a tumor along with organs using our SQH-QPAT framework. We consider
$\sigma$ that represents the organs as the same Shepp-Logan phantom of test
Case 4. The $\sigma$ for the tumors are represented by two disks: one of them
having a large contrast with value 1.5 and the other having a contrast similar
to the some of the organs with value 0.2. The corresponding value of $D$ is
0.001 and 0.004, respectively. The background value $\sigma_{b}=0.5$.

(a) Exact $D$

(b) Exact $\sigma_{a}$

(c) Convergence history of SQH-QPAT iteration

(d) Reconstructed $D$

(e) Reconstructed $D$ with 5% noise

(f) Reconstructed $D$ with 10% noise

(g) Reconstructed $\sigma_{a}$

(h) Reconstructed $\sigma_{a}$ with 5% noise

(i) Reconstructed $\sigma_{a}$ with 10% noise

Figure 7: Test Case 5: The actual and the reconstructed phantoms

For test Case 5, Figure 7 shows the exact and the reconstructed values of $D$
and $\sigma_{a}$ with and without noise in the interior data. We see that our
SQH-QPAT framework is able to reconstruct very well the tumors alongside the
organs, which further demonstrates the wide applicability of our
reconstruction framework. Figure 7c shows the convergence history of the SQH-
QPAT scheme. Notice the strict monotonic decrease in the functional values,
with stopping criteria fulfilled in 8 iterations. From Table 6, we also
observe similar RMSE % and PSNR values, as in the previous test case,
demonstrating high quality reconstructions.

Noise % | RMSE % ($D$) | RMSE % ($\sigma_{a}$) | PSNR ($D$) | PSNR ($\sigma_{a}$)
---|---|---|---|---
0 | 16.50 | 7.19 | 55.89 | 36.62
5 | 16.69 | 7.40 | 55.79 | 36.25
10 | 16.90 | 7.40 | 55.45 | 35.12
Table 6: Test Case 5: RMSE % and PSNR values for reconstructions of the Shepp-
Logan phantom with different noise levels in the data

##  6 Conclusion

A new optimization framework for reconstruction with high contrast and
resolution of optical diffusion and absorption coefficients in quantitative
photoacoustic tomography (QPAT) was presented. The resulting QPAT inverse
problem was formulated as the minimization of a best-fit functional with a
tentative prior relating the diffusion and absorption coefficients and $L^{2}$
and $L^{1}$ regularization terms, and subject to the constraint of the steady
photon-propagation equation.

For the solution of this problem, a robust and fast sequential quadratic
hamiltonian algorithm based on the Pontryagin maximum principle was developed
and analysed. Results of numerical experiments were presented that demonstrate
the ability of the proposed framework to efficiently compute reconstructions
of the optical coefficients with high contrast and resolution.

## Acknowledgments

This work was partially supported by the German BMBF-Project iDeLIVER :
Intelligent MR Diagnosis of the Liver by Linking Model and Data-Driven
Processes, and by the US National Science Foundation Grant No: 2309491, and
the University of Texas at Arlington, Research Enhancement Program Grant No:
2022-605.

##  Appendix A Proof of Theorem 4.1

We have

$H_{\epsilon}(x,y,D,\sigma,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})\leq
H_{\epsilon}(x,y,w,z,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})$

for all $w\in[D_{l},D_{r}],~{}z\in[\sigma_{l},\sigma_{r}]$. Hence, it follows
that

$\displaystyle
H_{\epsilon}(x,y,D,\sigma,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})\leq$
$\displaystyle
H_{\epsilon}(x,y,D^{k},\sigma^{k},D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})$
(23) $\displaystyle=$ $\displaystyle
H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})$

for all $(x,y)\in\Omega$. Let

$\displaystyle L(x,y,D,\sigma,u_{1},u_{2})=$
$\displaystyle\frac{\alpha}{2}\sum_{i=1}^{2}(\mathcal{H}(x,y,\sigma,u_{i})-G_{i}^{\delta}(x,y))^{2}+\frac{\xi_{1}}{2}\sigma^{2}(x,y)$
$\displaystyle+\frac{\xi_{2}}{2}\left(D(x,y)-\frac{1}{100~{}(\sigma(x,y)+\sigma_{b})}\right)^{2}+\gamma|\sigma(x,y)|.$

Then, we have

$J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})=\int_{\Omega}L(x,y,D,\sigma,u_{1},u_{2})-L(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})~{}dxdy.$
(24)

Adding and subtracting $D^{k}~{}\nabla q_{i}^{k}\cdot\nabla
u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}$ and $D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}$
to the right-hand side of (24), we obtain

$\displaystyle
J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})$
$\displaystyle=$
$\displaystyle\int_{\Omega}L(x,y,D,\sigma,u_{1},u_{2})~{}dxdy+\int_{\Omega}\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)~{}dxdy$
$\displaystyle-\int_{\Omega}\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)~{}dxdy$
$\displaystyle\int_{\Omega}$
$\displaystyle-L(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})~{}dxdy-\int_{\Omega}\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)~{}dxdy$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)~{}dxdy.$

This gives us

$\displaystyle
J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})$
$\displaystyle=$
$\displaystyle\int_{\Omega}H(x,y,D,\sigma,u_{1},u_{2},q_{1}^{k},q_{2}^{k})-H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})-\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)$
$\displaystyle+\sum_{i=1}^{2}\left(D^{k}~{}\nabla q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)~{}dxdy.$
$\displaystyle=$ $\displaystyle\int_{\Omega}H(x,y,D,\sigma,u_{1}^{k}+\delta
u_{1},u_{2}^{k}+\delta
u_{2},q_{1}^{k},q_{2}^{k})-H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})$
$\displaystyle-\sum_{i=1}^{2}\left(D^{k}~{}\nabla q_{i}^{k}\cdot\nabla
u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)+\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)~{}dxdy.$
$\displaystyle=$
$\displaystyle\int_{\Omega}H(x,y,D,\sigma,u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})+\Gamma\alpha\sum_{i=1}^{2}\left[-(\sigma^{k}+\sigma_{b})[\mathcal{H}(x,y,\sigma^{k},u_{i}^{k})-G^{\delta}_{i}]\delta
u_{i}+\dfrac{1}{2}(\sigma^{k}+\sigma_{b})(\delta u_{i})^{2}\right]$
$\displaystyle-H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})-\sum_{i=1}^{2}\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)$
$\displaystyle+\sum_{i=1}^{2}\left(D^{k}~{}\nabla q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)~{}dxdy$

$\displaystyle=$
$\displaystyle\int_{\Omega}H(x,y,D,\sigma,u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})+\epsilon[(\delta
D)^{2}+(\delta\sigma)^{2}]-H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})~{}dxdy$
$\displaystyle+\int_{\Omega}\sum_{i=1}^{2}\Bigg{\\{}\Gamma\alpha(\sigma^{k}+\sigma_{b})[\mathcal{H}(x,y,\sigma^{k},u_{i}^{k})-G^{\delta}_{i}]u_{i}+\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}\right)\Bigg{\\}}~{}dxdy$
$\displaystyle-\int_{\Omega}\sum_{i=1}^{2}\Bigg{\\{}\Gamma\alpha(\sigma^{k}+\sigma_{b})[\mathcal{H}(x,y,\sigma^{k},u_{i}^{k})-G^{\delta}_{i}]u_{i}^{k}+\left(D^{k}~{}\nabla
q_{i}^{k}\cdot\nabla
u_{i}^{k}+(\sigma^{k}+\sigma_{b})~{}q_{i}^{k}u_{i}^{k}\right)\Bigg{\\}}~{}dxdy$
$\displaystyle+\sum_{i=1}^{2}\int_{\Omega}\dfrac{1}{2}(\sigma^{k}+\sigma_{b})(\delta
u_{i})^{2}-\epsilon\,[(\delta D)^{2}+(\delta\sigma)^{2}]~{}dxdy$
$\displaystyle=$
$\displaystyle\int_{\Omega}H_{\epsilon}(x,y,D,\sigma,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})-H(x,y,D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k},q_{1}^{k},q_{2}^{k})~{}dxdy$
$\displaystyle+\sum_{i=1}^{2}\int_{\Omega}\dfrac{1}{2}(\sigma^{k}+\sigma_{b})(\delta
u_{i})^{2}~{}dxdy-\int_{\Omega}\epsilon\,[(\delta
D)^{2}+(\delta\sigma)^{2}]~{}dxdy,\qquad\mbox{ (using }\eqref{eq:DA_adj})$
$\displaystyle\leq$
$\displaystyle\sum_{i=1}^{2}\int_{\Omega}\dfrac{1}{2}(\sigma^{k}+\sigma_{b})(\delta
u_{i})^{2}~{}dxdy-\int_{\Omega}\epsilon\,[(\delta
D)^{2}+(\delta\sigma)^{2}]~{}dxdy,\qquad\mbox{ (using }\eqref{sqh6_6})$

Now, we recall the result of Lemma 3.1 that gives the following

$\|\delta u_{i}\|_{L^{2}(\Omega)}^{2}\leq\bar{C}_{1}\,\big{(}\|\delta
D\|_{L^{2}(\Omega)}^{2}+\|\delta\sigma\|_{L^{2}(\Omega)}^{2}\big{)},$

where $\bar{C}_{1}=2\,c_{3}^{2}/c_{1}^{2}$. Hence, we proceed as follows

$\displaystyle
J(D,\sigma,u_{1},u_{2})-J(D^{k},\sigma^{k},u_{1}^{k},u_{2}^{k})\leq$
$\displaystyle\int_{\Omega}-\epsilon\left[(\delta
D)^{2}+(\delta\sigma)^{2}\right]~{}dxdy+\sum_{i=1}^{2}\int_{\Omega}\dfrac{1}{2}(\sigma^{k}+\sigma_{b})(\delta
u_{i})^{2}\,dx$ $\displaystyle\leq$ $\displaystyle-\epsilon\left[\|\delta
D\|_{L^{2}(\Omega)}^{2}+\|\delta\sigma\|_{L^{2}(\Omega)}^{2}\right]+\sum_{i=1}^{2}\frac{\bar{\sigma}}{2}\|\delta
u_{i}\|_{L^{2}(\Omega)}^{2}$ $\displaystyle\leq$
$\displaystyle-\epsilon\left[\|\delta
D\|_{L^{2}(\Omega)}^{2}+\|\delta\sigma\|_{L^{2}(\Omega)}^{2}\right]+\bar{C}_{1}\,\bar{\sigma}\,\left[\|\delta
D\|_{L^{2}(\Omega)}^{2}+\|\delta\sigma\|_{L^{2}(\Omega)}^{2}\right]$
$\displaystyle=$
$\displaystyle-\left(\epsilon-\bar{C}_{1}\,\bar{\sigma}\right)\,\left[\|\delta
D\|_{L^{2}(\Omega)}^{2}+\|\delta\sigma\|_{L^{2}(\Omega)}^{2}\right],$

where $\bar{\sigma}=\sigma_{r}+\sigma_{b}$. Thus, the theorem is proved with
$\theta=\bar{C}_{1}\,\bar{\sigma}$.

## References

  * [1] G. S. Abdoulaev, K. Ren and A. H. Hielscher. Optical tomography as a PDE-constrained optimization problem, Inverse Problems, 21:1507-1530, 2005.
  * [2] B. Adesokan, K. Knudsen, V. P. Krishnan and S. Roy. A fully non-linear optimization approach to acousto-electric tomography. Inverse Problems, 34:104004, 2018.
  * [3] G. S. Alberti and Y. Capdeboscq. Lectures on Elliptic Methods for Hybrid Inverse Problems, Société Mathématique de France, Paris, 2018.
  * [4] A. Borzì and V. Schulz. Computational optimization of systems governed by partial differential equations, SIAM Computational Science and Engineering, 2012.
  * [5] H. W. Alt. Linear Functional Analysis. An Application-Oriented Introduction, Springer, 2016.
  * [6] R. Aronson and N. Corngold. Photon diffusion coefficient in an absorbing medium, Journal of Optical Society of America, 16(5):1066-1071, 1999.
  * [7] U. M. Ascher and E. Haber, A multigrid method for distributed parameter estimation problems, ETNA, 18 (2003), 1-18.
  * [8] H. Ammari, E. Bossy, V. Jugnon and H. Kang. Reconstruction of the optical absorption coefficient of a small absorber from the absorbed energy density, SIAM Journal of Applied Mathematics, 71:676-693, 2011.
  * [9] S. R. Arridge. Optical tomography in medical imaging, Inverse Problems, R41-R93, 1999.
  * [10] G. Bal, A. Jollivet and V. Jugnon. Inverse transport theory of photoacoustics, Inverse Problems, 26(2):025011, 2010.
  * [11] G. Bal and K. Ren. Multi-source quantitative photoacoustic tomography in a diffusive regime, Inverse Problems, 27:075003, 2011.
  * [12] G. Bal and K. Ren. On multi-spectral quantitative photoacoustic tomography in diffusive regime, Inverse Problems, 28:025010, 2012.
  * [13] G. Bal and G. Uhlmann. Inverse diffusion theory of photoacoustics, Inverse Problems, 26:085010, 2010.
  * [14] G. Bal and G. Uhlmann. Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions, Communications in Pure and Applied Mathematics, 66:1692-1752, 2013.
  * [15] B. Banerjee, S. Bagchi, R. M. Vasu and D. Roy. Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map, Journal of Optical Society of America A, 25(9):2347-2356, 2008.
  * [16] H.T. Banks and K. Kunisch. Estimation Techniques for Distributed Parameter Systems, Birkhäuser, Basel, 1989.
  * [17] N. Bosschaart, G. J. Edelman, M. C. G. Aalders, T. G. van Leeuwen and D. J. Faber. A literature review and novel theoretical approach on the optical properties of whole blood, Lasers in Medical Science, 29:453–479, 2014.
  * [18] T. Breitenbach and A. Borzì. A sequential quadratic Hamiltonian method for solving parabolic optimal control problems with discontinuous cost functionals, Journal of Dynamical and Control Systems, 25(3):1-33, 2018.
  * [19] T. Breitenbach and A. Borzì. On the SQH Scheme to solve nonsmooth PDE optimal control problems, Numerical Functional Analysis and Optimization, 40(13): 1489-1531, 2019.
  * [20] T. Breitenbach and A. Borzì. A sequential quadratic Hamiltonian scheme for solving non-smooth quantum control problems with sparsity, Journal of Computational and Applied Mathematics, 369:112583, 2019.
  * [21] F. L. Chernous’ko, A. A. Lyubushin. Method of successive approximations for solution of optimal control problems, Optimal Control Applications and Methods, 3(2):101 - 114, 2007.
  * [22] A. Cianchi and V. Maz’ya. Global gradient estimates in elliptic problems under minimal data and domain regularity. Communications on Pure & Applied Analysis, 14:285, 2015.
  * [23] F. H. Clarke. Optimization and Nonsmooth Analysis. John Wiley & Sons, New York, 1983.
  * [24] B. T. Cox, S. R. Arrdige, P. Köstli and P. C. Beard. Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method, Applied Optics, 45:1866-1875, 2006.
  * [25] B. T. Cox, S. R. Arrdige and P. C. Beard. Estimating chromophore distributions from multiwavelength phtoacoustic images, Journal of Optical Society of America, 26:443-455, 2009.
  * [26] A. Dimofte, J. Finlay and T. Zhu. A method for determination of the absorption and scattering properties interstitially in turbid media, Physics in Medicine and Biology, 50(10): 2291–2311, 2005.
  * [27] D. J. Durian. The diffusion coefficient depends on absorption, Optics Letters, 23(19):1502-1504, 1998.
  * [28] V. Dzimbeg-Malcic, Zeljka Barbaric-Mikocevic and Katarina Itric. Kubelka-Munk theory in describing optical properties of paper (I), Tehnicki Vjesnik, 19(1):191-196, 2012.
  * [29] F. Fanjul-Vélez, S. Pampín-Suárez and J. L. Arce-Diego. Application of classification algorithms to diffuse reflectance spectroscopy measurements for ex vivo characterization of biological tissues, Entropy, 22:736, 2020.
  * [30] A diffusion theory model of spatially resolved, steady-state diffuse reflectance for the noninvasive determination of tissue optical properties in vivo, Medical Physics, 19(4):879-888, 1992.
  * [31] D. Finch and Rakesh. Recovering a function from its spherical mean values in two and three dimensions, Photoacoustic Imaging and Spectroscopy, ed. L. H. Wang, CRC Press, 2009.
  * [32] D. Gao, Y. Liu, Y. Wang and Z. Yuan. Protein-modified ultra-small gold clusters for dual-modal in vivo fluorescence/photoacoustic imaging, Quantitative Imaging in Medicine and Surgery, 8(3):326-332, 2018.
  * [33] H. Gao, J. Feng and L. Song. Limited-view multi-source quantitative photoacoustic tomography, Inverse Problems, 31:065004, 2015.
  * [34] H. Gao, S. Osher and H. Zhao. Quantitative photoacoustic tomography, Mathematical Modeling in Biomedical Imaging II, Lecture Notes in Mathematics 2035, Springer, Heidelberg, 131-158, 2012.
  * [35] D. Gilbarg and N. Trudinger. Elliptic Partial Differential Equations of Second Order. Springer-Verlag Berlin Heidelberg New York, 1998.
  * [36] M. Haltmeier, L. Neumann, L. Nguyen and S. Rabanser. Analysis of the linearized problem of quantitative photoacoustic tomography, SIAM Journal of Applied Mathematics, 78(1):457-478, 2018.
  * [37] H. G. Hecht. The interpretation of diffuse reflectance spectra, Journal of Research of the National Bureau of Standards-A. Physics and Chemistry, 80A(4):567-583, 1976.
  * [38] A. Ishimaru. Wave Propagation and Scattering in Random Media , New York, Academic Press, 1978.
  * [39] P. E. Kinahan and J. S. Karp. Figures of merit for comparing reconstruction algorithms with a volume-imaging PET scanner, Physics in Medicine and Biology, 39(3):631-642, 1994.
  * [40] P. Kuchment and D. Steinhauer. Stabilizing inverse problems by internal data, Inverse Problems, 28(8):084007, 2012.
  * [41] P. Kubelka and F. Munk. Ein Beitrag zur Optik der Farbanstriche, Z. Tech. Phys., 12:593-601, 1931.
  * [42] O. A. Ladyzhenskaya and N. N. Ural’tseva. Linear and Quasilinear Elliptic Equations, Academic Press, New York and London, 1968.
  * [43] J. Laufer, B. Cox, E. Zhang and P. Beard. Quantitative determination of chromophore concentrations from 2d photoacoustic images using a nonlinear model-based inversion scheme, Applied Optics, 49:1219-1233, 2010.
  * [44] R. Leitao and E. Teixeira. Regularity and geometric estimates for minima of discontinuous functionals, Rev. Mat. Iberoam, 31(1):69-108, 2015.
  * [45] G. Mingione. Regularity of minima: An invitation to the dark side of the calculus of variations, Applications of Mathematics, 51:355-426, 2006.
  * [46] V. Ntziachristos, Fluorescence molecular imaging, Annual Review of Biomedical Engineering, 8:1-33, 2006.
  * [47] K. Ren, H. Gao and H. Zhao. A hybrid reconstruction method for quantitative PAT, SIAM Journal of Imaging Sciences, 6:32-55, 2013.
  * [48] A. Roy, R. Ramasubramaniam and H. A. Gaonkarb. Empirical relationship between Kubelka–Munk and radiative transfer coefficients for extracting optical parameters of tissues in diffusive and nondiffusive regimes, Journal of Biomedical Optics, 17(11):115006, 2012.
  * [49] S. Roy and A. Borzì. A new optimisation approach to sparse reconstruction of log-conductivity in acousto-electric tomography, SIAM Journal on Imaging Sciences, 11(2):1759–1784, 2018.
  * [50] Y. Sakawa and Y. Shindo. On global convergence of an algorithm for optimal control. IEEE Transactions on Automatic Control, 25(6):1149–1153, 1980.
  * [51] J. L. Sandell and T. C. Zhu. A review of in-vivo optical properties of human tissues and its impact on PDT, Journal of Biophotonics, 4(11-12): 773–787, 2011.
  * [52] T. Saratoon, T. Tarvainen, B. T. Cox and S. R. Arrdige. A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation, Inverse Problems, 29:075006, 2013.
  * [53] J. C. Schotland, Direct reconstruction methods in optical tomography, Lecture Notes in Mathematics, Vol. 2035, Springer-Verlag, Berlin, 1-29, 2011.
  * [54] P. Shao, B. Cox and R. J. Zemp. Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography, Applied Optics, 50:3145-3154, 2011.
  * [55] L. A. Shepp and B. F. Logan. The Fourier reconstruction of a head section.  IEEE Transactions on Nuclear Science, 21(3):21–43, 1974.
  * [56] G. Stadler. Elliptic optimal control problems with $L^{1}$-control cost and applications for the placement of control devices, J Comp Opt and App, 44:159-181, 2007.
  * [57] Light scattering methods for tissue diagnosis, Optica, 6(4):479-489, 2019.
  * [58] M.I. Sumin. Optimal control of objects described by quasilinear elliptic equations, Differ. Equations, 25: 1004-1012, 1989.
  * [59] T. Tarvainen, B. T. Cox, J. P. Kaipio and S. R. Arridge. Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography, Inverse Problems, 28:084009, 2012.
  * [60] M. Ulbrich. Semismooth Newton Methods for Variational Inequalities and Constrained Optimization Problems in Function Spaces, SIAM, Philadelphia, 2011.
  * [61] L. V. Wang. Multiscale photoacoustic microscopy and computed tomography, Nature Photonics, 3:503-509, 2009.
  * [62] H-J. Wei, D. Xing, G-Y. Wu, Y. Jin and H-M. Gu. Optical properties of human normal small intestine tissue determined by Kubelka-Munk method in vitro, World Journal of Gastroenterology, 9(9):2068-2072, 2003.
  * [63] J. Xia, J. Yao, and L. V. Wang. Photoacoustic tomography: principles and advances, Electromagnetic Waves, 147:1-22, 2014.
  * [64] L. Yao, Y. Sun and H. Jiang. Transport-based quantitative photoacoustic tomography: simulations and experiments, Physics in Medicine and Biology, 55:1917-1934, 2010.
  * [65] R. J. Zemp. Quantitative photoacoustic tomography with multiple optical sources, Applied Optics, 49:3566-3572, 2010.