# Quantitative Photoacoustic Tomography Using Iteratively Refined Wavefield Reconstruction Inversion: A Simulation Study

Seyed Mohsen Ranjbaran, Hossein Aghamiry, Ali Gholami, Stephane Operto, Kamran Avanaki

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

► To cite this version:

Seyed Mohsen Ranjbaran, Hossein Aghamiry, Ali Gholami, Stephane Operto, Kamran Avanaki. Quantitative Photoacoustic Tomography Using Iteratively Refined Wavefield Reconstruction Inversion: A Simulation Study. IEEE Transactions on Medical Imaging, 2024, 43 (2), pp.874-885. 10.1109/TMI.2023.3324922. hal-04922071

HAL Id: hal-04922071

https://hal.science/hal-04922071v1

Submitted on 31 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Quantitative Photoacoustic Tomography Using Iteratively Refined Wavefield Reconstruction Inversion: A Simulation Study

Seyed Mohsen Ranjbaran<sup>D</sup>, Hossein S. Aghamiry<sup>D</sup>, Member, IEEE, Ali Gholami<sup>D</sup>, Stéphane Operto, and Kamran Avanaki<sup>D</sup>

**Abstract**— The ultimate goal of photoacoustic tomography is to accurately map the absorption coefficient throughout the imaged tissue. Most studies either assume that acoustic properties of biological tissues such as speed of sound (SOS) and acoustic attenuation are homogeneous or fluence is uniform throughout the entire tissue. These assumptions reduce the accuracy of estimations of derived absorption coefficients (DeACs). Our quantitative photoacoustic tomography (qPAT) method estimates DeACs using iteratively refined wavefield reconstruction inversion (IR-WRI) which incorporates the alternating direction method of multipliers to solve the cycle skipping challenge associated with full wave inversion algorithms. Our method compensates for SOS inhomogeneity, fluence decay, and acoustic attenuation. We evaluate the performance of our method on a neonatal head digital phantom.

**Index Terms**— Iteratively refined wavefield reconstruction inversion, alternating direction method of multipliers, fluence compensation, quantitative photoacoustic tomography, speed of sound estimation.

## I. INTRODUCTION

PHOTOACOUSTIC imaging (PAI) is an emerging hybrid imaging modality, combining the advantages of optical contrast with acoustic penetration. PAI has been used extensively in preclinical [1], [2], and more recently, clinical studies [3], [4]. As shown in Figure 1, upon nanosecond

Manuscript received 13 April 2023; revised 30 August 2023; accepted 2 October 2023. Date of publication 17 October 2023; date of current version 2 February 2024. The work of Hossein S. Aghamiry was supported by the WIND Consortium (https://www.geoazur.fr/WIND). The work of Kamran Avanaki was supported by the National Institutes of Health under Grant R01EB027769 and Grant R01EB028661. (Corresponding author: Kamran Avanaki.)

Seyed Mohsen Ranjbaran is with the Richard and Loan Hill Department of Biomedical Engineering, University of Illinois at Chicago, Chicago, IL 60607 USA (e-mail: sranjb3@uic.edu).

Hossein S. Aghamiry is with the Geoazur, CNRS, IRD, OCA, University Côte d'Azur, 06560 Valbonne, France, and also with the Center for Biomedicine, Charité—Universitätsmedizin, 10117 Berlin, Germany (e-mail: Hossein.aghamiry@charite.de).

Ali Gholami is with the Institute of Geophysics, Polish Academy of Sciences, 01-452 Warsaw, Poland (e-mail: agholami@igf.edu.pl).

Stéphane Operto is with the Geoazur, CNRS, IRD, OCA, University Côte d'Azur, 06560 Valbonne, France (e-mail: operto@geoazur.unice.fr).

Kamran Avanaki is with the Richard and Loan Hill Department of Biomedical Engineering, the Department of Dermatology, and the Department of Pediatrics, University of Illinois at Chicago, Chicago, IL 60607 USA (e-mail: avanaki@uic.edu).

Digital Object Identifier 10.1109/TMI.2023.3324922

<figure><img src="image_3.png" /><figcaption>Fig. 1. In photoacoustic imaging, a pulsed laser illuminates tissue, and acoustic pressure waves (PA waves) are generated by chromophores in the tissue. These PA waves propagate through the tissue and are measured by ultrasound (US) transducers.</figcaption></figure>

laser irradiation of tissue, chromophores such as oxyhemoglobin (HbO2) and deoxy-hemoglobin (Hb) absorb energy and generate photoacoustic waves through the thermoelastic effect, if thermal and stress confinements are met [5], [6]. The generated acoustic waves are detected by ultrasound (US) transducers and processed by an image reconstruction algorithm (e.g. back projection, time reversal, beamforming, Fourier transform [7], [8], [9], [10], [11]) to restore the induced initial pressure or equivalently, an optical absorption map.

Generally, photoacoustic (PA) image reconstruction algorithms rely on the assumption that the imaging area is homogeneous [12], [13], [14], which implies that the speed of sound (SOS) is constant throughout the imaging region. This is an invalid assumption in most biological tissue (for example, brain tissue has many regions with different SOS) [15], [16], leading to distortion in the received PA signals, and consequently in reconstructed PA images. There have been several studies that aim to implement SOS estimation [12], [14], [17], [18], [19]. Xu and Wang [12] used an iterative algorithm to compensate for SOS variation in breast tissue. Deán-Ben et al. [13] used prior information on the position of acoustic deformities, which is inferred from usual anatomical data or another imaging method to improve image quality. Zhang and Wang [19] proposed an algorithm based on the correlation of PA recorded signals to neutralize heterogeneity in acoustical properties. Treeby et al. [14] used the autofocus
strategy to automatically estimate SOS in the medium by maximizing image sharpness. In Jose et al. [17], the reconstruction algorithm relies on refracted ray paths rather than direct ray paths to obtain more realistic SOS maps. Recently, Yin et al. [18] developed an algorithm that utilizes an iterative procedure to estimate an SOS map relying on the fact that the recorded signals by two neighboring ultrasonic transducers are very similar in the received waveforms. Sun and Sun [20] developed a method to estimate both the optical absorption coefficient and the SOS map simultaneously.

In contrast to methods that either use part of the recorded data, e.g., travel time or amplitude, or make unrealistic assumptions about the medium (e.g. constant speed of sound, homogeneity, etc.), full waveform inversion (FWI) exploits the full information contained in the recorded wavefields to estimate the constitutive properties of a medium without any assumption about its level of heterogeneity or requirement of a level of homogeneity [21], [22]. This enriched information guarantees superior performance of wave-based tomography relative to travel time-based approaches in terms of ability to estimate several constitutive properties such as SOS, density, and acoustic attenuation. Considering emerging high-performance computing facilities and the current state of theoretical understanding of this methodology, it is reasonable to anticipate that wave-based tomography will rapidly become the preferred method to perform accurate tomography [7]. Algorithms using full waveform information are capable of producing images with higher resolution (on the order of one wavelength) than algorithms using ray-based methods [23].

FWI for brain imaging is challenging due to both large variability in SOS within the brain tissue and the high SOS difference between the soft brain tissue and the surrounding skull. Classical FWI, which is based on the iterative least-squares minimization of the differences between the measured data and the numerically simulated counterparts, requires prior information about the skull to reconstruct an accurate SOS map of the brain. To relax this requirement, several approaches extend the search space of classical FWI by incorporating more convex distances between the recorded and simulated data. One such approach is a so-called adaptive waveform inversion (AWI) which has been used to reconstruct the SOS map of the human brain without resorting to prior information about the skull [24].

There have been several attempts to use FWI in PA imaging to compute PA image or the absorption coefficient (referred to in this study as derived absorption coefficient or DeAC) [25], [26], [27], [28], [29]. Computing the DeAC (referred to as quantitative photoacoustic tomography or qPAT) is challenging, particularly for biological tissues where optical fluence distribution is unknown [30], [31], [32]. In qPAT, the optical inverse problem needs to be addressed [33], [34]. This problem is typically tackled using iterative approaches based on mathematical models, as described in [35], [36], [37], [38], [39], [40], and [41]. The model employed in these approaches is based on the radiative transfer equation (RTE), which accurately simulates the propagation of light to calculate optical fluence [37], [42], [43]. The optical fluence distribution

throughout the tissue is calculated based on optical properties of different tissue regions as well as specifications and geometrical characteristics of the laser [44], [45], [46]. In some qPAT studies, a simplified version of the RTE called the diffusion approximation is applied [35], [36], [37], [47], [48], [49]. This approximation is generally satisfactory when dealing with highly scattering media and nearly uniform scattering in all directions [50]. However, these assumptions are not applicable to studies involving the human head. To ensure uniqueness in the inversion process, three different approaches are commonly utilized [48]: solving for multiple optical wavelengths [47], [51], [52], assuming prior knowledge of the scattering coefficient [38], [53], [54], or employing multiple optical sources [36], [55], [56]. In the current study, we assume knowledge of the scattering coefficients and anisotropy factors are available. Testing several fluence compensation methods [34], [57], we observed that the correction works only at shallow depths, i.e., superficial layers of the target.

In this paper, we propose a qPAT method that calculates DeACs by estimating SOS, acoustic attenuation, and fluence maps. Our method uses iteratively refined wavefield reconstruction inversion (IR-WRI), which incorporates the alternating direction method of multipliers (ADMM) to solve the cycle skipping challenge associated with FWI algorithms. We evaluate the performance of our qPAT method on a digital phantom of a neonatal head.

## II. NOTATION

The mathematical symbols adopted in this paper are as follows. We use italics for scalar quantities, boldface lowercase letters for vectors, and boldface capital letters for matrices and tensors. We use the superscript *T* to denote the adjoint of an operator. The *i*th component of the column vector **x** is shown by $x_i$ and its absolute value is returned by $|x_i|$. For the real-valued *n*-length column vectors **x** and **y**, the dot product is defined by $\langle x, y \rangle = x^T y = \sum_{i=1}^{n} x_i y_i$ and their Hadamard product, denoted by $x \circ y$, is another vector made up of their component-wise products, i.e. $(x \circ y)_i = x_i y_i$. The $\ell$2-and $\ell$1-norms of **x** are, respectively, defined by $||x||_2 = \sqrt{\langle x, x \rangle} = \sqrt{\sum_{i=1}^{n} x_i^2}$ and $||x||_1 = \sum_{i=1}^{n} |x_i|$, and, $|| \cdot ||_F$ denotes the Frobenius norm. Finally, for a complex number $z = \Re(z) + i\Im(z)$, $\Re(z)$ and $\Im(z)$ refer to the real and imaginary parts, and $i = \sqrt{-1}$.

## III. METHOD

### A. Digital Phantom of Neonate Head and Simulation Setup

We created a digital phantom of a 30-week-old neonatal head by averaging 324 MRI atlases made by Brigadoi et al. [58]. The model is composed of extra-cranial tissue (ECT), cerebrospinal fluid (CSF), gray matter (GM), and white matter (WM). Assigned optical properties for all layers were taken from the literature [59], [60], [61]. ECT includes scalp and skull; the averaged optical properties of scalp and skull were assigned to this layer. To fully mimic the neonatal head, we included complex vasculature as the primary imaging target of the phantom. The size of the computational
<figure><img src="image_4.png" /><figcaption>Fig. 2. Configuration of the simulated US/PA tomography system. US data collection proceeds as follows. Each of the 128 transducers is activated singly, with the other 127 transducers in receive mode. ECT: extra-cranial tissue, CSF: cerebrospinal fluid, GM: gray matter, WM: white matter.</figcaption></figure>

grid for this head model is 1000 × 1000 pixels, in which pixel size is 120 microns. In this work, 64 optical fibers with size of 0.6 mm are equally interspersed between 128 US transducers (active aperture = 1 mm) to form a ring (diameter of 120 mm) around the digital phantom. Fig. 2 shows the regions in the digital phantom of the neonatal head as well as the locations and arrangement of the optical fibers and US transducers.

## B. SOS and Acoustic Attenuation Estimation

We implemented FWI in the frequency domain as opposed to the more conventional time domain formulation. The frequency domain decreases the nonlinearity of the problem, which helps to bypass the local minima (cycle-skipping) problem. It is also computationally efficient when the inversion can be limited to a few discrete frequencies. This is the case when the acquisition device (in our case, a ring array) provides a redundancy of information due to abundant illumination of the target in both the US and PA setups. The frequency domain FWI can be formulated as the following multivariate constrained optimization problem [62]

$$ \begin{array}{ll} \text{minimize} & \mathcal{R}(\mathbf{m}) \\ \text{subject to} & \left\{ \begin{array}{l} \mathbf{A}(\mathbf{m}, \omega)\mathbf{U}(\omega) = \mathbf{B}(\omega), \\ \mathbf{P}\mathbf{U}(\omega) = \mathbf{D}(\omega), \\ \mathbf{U}(\omega), \mathbf{m} \in \mathcal{M} \quad \omega = \omega_1, \omega_2, \dots, \omega_{n_f} \end{array} \right. \end{array} \quad (1) $$

where $\mathcal{R}(\mathbf{m})$ is a regularization function and the constraints are the wave equation $\mathbf{A}(\mathbf{m}, \omega)\mathbf{U}(\omega) = \mathbf{B}(\omega)$ and the observation equation $\mathbf{P}\mathbf{U}(\omega) = \mathbf{D}(\omega)$. Moreover, $\mathbf{m} \in \mathbb{C}^{N \times 1}$ denotes complex-valued model parameters, constructed by a mapping between SOS and acoustic attenuation. The mathematical model describing the attenuating medium is given in Appendix A. In addition, $N$ is the number of discretized points of the medium, $\omega$ is the angular frequency, $n_f$ is the number of frequencies, $\mathbf{U}(\omega) = [u_1(\omega), u_2(\omega), \dots, u_{n_s}(\omega)] \in \mathbb{C}^{N \times n_s}$ is a matrix, whose columns represent the wavefields

triggered by each source, $\mathbf{A}(\mathbf{m}, \omega) = \Delta + \omega^2 \text{diag}(\mathbf{m}) \in \mathbb{C}^{N \times N}$ is the impedance matrix in the acoustic environment, $\Delta$ is the discretized Laplace matrix, and $\text{diag}(\cdot)$ is a diagonal matrix with $\cdot$ in its diagonal. In addition, $\mathbf{B}(\omega) = [\mathbf{b}_1(\omega), \mathbf{b}_2(\omega), \dots, \mathbf{b}_{n_s}(\omega)] \in \mathbb{C}^{N \times n_s}$ is a matrix, the columns of which contain the $n_s$ sources, $\mathbf{D}(\omega) = [\mathbf{d}_1(\omega), \mathbf{d}_2(\omega), \dots, \mathbf{d}_{n_s}(\omega)] \in \mathbb{C}^{n_r \times n_s}$ is the matrix of recorded data, $n_r$ is the number of receivers and $\mathbf{P} \in \mathbb{R}^{n_r \times N}$ is the linear observation (sampling) operator that samples the wavefields $\mathbf{U}$ at receiver positions. Finally, $\mathcal{M}$ is a convex set defined according to our prior knowledge of $\mathbf{m}$. For example, if we know the lower and upper bounds on $\mathbf{m}$ then

$$ \mathcal{M} = \{ \mathbf{m} \mid \mathbf{m}_{\min} \le \mathbf{m} \le \mathbf{m}_{\max} \}. \quad (2) $$

The level of non-linearity of FWI in the optimization problem (1) is moderate [63], which means that local optimization methods, like Newton-family solvers [64] even with an inaccurate starting model, converge to a reasonable solution. To solve the optimization problem (1) in the full parameter space, we should deal simultaneously with all the optimization variables $(\mathbf{U}(\omega), \mathbf{m})$. This approach which is called all-at-once [65] is not feasible for practical applications of FWI because of massive computational burden and memory requirements. The classical formulation of FWI was implemented on the reduced parameter space to keep the computational burden and memory requirements reasonable [22], [66]. In classical FWI, the wave-equation constraint is solved for $\mathbf{U}(\omega)$, i.e. $\mathbf{U}(\omega) = \mathbf{A}(\mathbf{m}, \omega)^{-1}\mathbf{B}(\omega)$, and implemented in the second constraint, i.e. $\mathbf{P}\mathbf{A}(\mathbf{m}, \omega)^{-1}\mathbf{B}(\omega) = \mathbf{D}(\omega)$, which implies that the optimization variables reduce down to the constitutive parameters of the medium. In this case, the inverse problem is highly nonlinear due to dependency on $\mathbf{A}(\mathbf{m}, \omega)^{-1}$, meaning that local optimization methods require a good starting model to converge to the global minimum. Without a good starting model, the problem is trapped in a local minimum, a process which is called cycle skipping [67]. In order to solve the cycle skipping challenge of FWI, we recently proposed IR-WRI [68], [69] which uses an all-at-once approach for (1) in which the optimization variables are updated in an alternating mode based on ADMM [70] instead of simultaneous update of all optimization variables. IR-WRI mitigates cycle skipping by keeping the wavefield as an independent optimization variable. In this method, first the wavefields that jointly satisfy the observation equation and the wave equation are computed in a least-squares sense. Then, the constitutive parameters of the medium are updated by minimizing the wave-equation errors. Then, the Lagrange multipliers, which are required for solving a constraint optimization, are updated. We solve (1) with an augmented Lagrangian method [64] whose objective function is given by:

$$ \begin{aligned} \mathcal{L}_A(\mathbf{m}, \mathbf{U}(\omega), \mathbf{V}(\omega), \mathbf{W}(\omega)) & = \mathcal{R}(\mathbf{m}) - \sum_{\omega=\omega_1}^{\omega_{n_f}} \left\langle \mathbf{V}(\omega)^T, \mathbf{A}(\mathbf{m}, \omega)\mathbf{U}(\omega) - \mathbf{B}(\omega) \right\rangle \\ & \quad - \sum_{\omega=\omega_1}^{\omega_{n_f}} \left\langle \mathbf{W}(\omega)^T, \mathbf{P}\mathbf{U}(\omega) - \mathbf{D}(\omega) \right\rangle \end{aligned} $$
$$ 
\begin{aligned}
&+ \frac{\lambda_1}{2} \sum_{\omega=\omega_1}^{\omega_{n_f}} ||\mathbf{A}(\mathbf{m}, \omega)\mathbf{U}(\omega) - \mathbf{B}(\omega)||_F^2 \\
&+ \frac{\lambda_2}{2} \sum_{\omega=\omega_1}^{\omega_{n_f}} ||\mathbf{P}\mathbf{U}(\omega) - \mathbf{D}(\omega)||_F^2,
\end{aligned}
\quad (3) $$

where $\mathbf{V}(\omega) \in \mathbb{C}^{N \times n_s}$ and $\mathbf{W}(\omega) \in \mathbb{C}^{n_r \times n_s}$ are Lagrange multipliers and the scalars $\lambda_1, \lambda_2 > 0$ are the penalty parameters assigned to the wave equation and the observation equation constraints, respectively. The primal variables $\mathbf{m}$ and $\mathbf{U}(\omega)$ and the dual variables (i.e., the Lagrange multipliers) $\mathbf{V}(\omega)$ and $\mathbf{W}(\omega)$ can be updated in alternating mode within the framework of ADMM to break the full multivariate problem into several manageable subproblems [68]. The stationary point of (3) is determined iteratively as

$$ \mathbf{U}(\omega)^{k+1} = \arg \min_{\mathbf{U}(\omega)} \mathcal{L}_A(\mathbf{m}^k, \mathbf{U}(\omega), \mathbf{V}(\omega)^k, \mathbf{W}(\omega)^k), \quad (4a) $$

$$ \mathbf{m}^{k+1} = \arg \min_{\mathbf{m}} \sum_{\omega} \mathcal{L}_A(\mathbf{m}, \mathbf{U}(\omega)^{k+1}, \mathbf{V}(\omega)^k, \mathbf{W}(\omega)^k), \quad (4b) $$

$$ \mathbf{V}(\omega)^{k+1} = \mathbf{V}(\omega)^k + \lambda_1[\mathbf{B}(\omega) - \mathbf{A}(\mathbf{m}^{k+1}, \omega)\mathbf{U}(\omega)^{k+1}], \quad (4c) $$

$$ \mathbf{W}(\omega)^{k+1} = \mathbf{W}(\omega)^k + \lambda_2[\mathbf{D}(\omega) - \mathbf{P}\mathbf{U}(\omega)^{k+1}], \quad (4d) $$

where $\omega = \omega_1, \omega_2, \dots, \omega_{n_f}$. The wavefield $\mathbf{U}(\omega)^{k+1}$ in (4a) can be obtained in a closed-form expression as given in [68]. Once the wavefields are reconstructed, the model parameters are updated by minimizing the wave-equation errors keeping the wavefields fixed (4b). The ill-posedness of any FWI method, including IR-WRI, also requires regularization to remove the undesired features in the reconstructed model that may be due to inhomogeneous illumination, sparse acquisition, or parameter cross-talk during multiparameter reconstruction [71].

A suitable regularization renders the solution unique, increases its stability, and prevents data over-fitting that may cause coherent noise in the reconstructed model. The regularization should be able to (mathematically) describe the solution while being easy to implement [72]. Conventional regularizations such as those promoting smoothness, sparseness, and blockiness are derived according to some prior assumptions about the targeted model and are independent of the data. In contrast, plug-and-play regularizations (such as those based on nonlocal means filters and block matching 3D filters (BM3D) [73]) are more flexible for regularizing complex models. Another rationale for using BM3D is that proximal methods, such as ADMM [70], provide a versatile and computationally efficient approach to implement non-differentiable and hybrid regularization functions (such as BM3D) through variable splitting, as reviewed in [71]. Here, we use BM3D regularization for SOS and acoustic attenuation estimation. SOS and acoustic attenuation can be easily represented by one complex-valued frequency-dependent term (see Appendix A). In [69], Aghamiry et al. presents a method to process this complex-valued parameter as the optimization variable in IR-WRI, instead of processing SOS and attenuation as two

<table><caption>TABLE I
INPUT PARAMETERS FOR MCX SIMULATION</caption><thead><tr><td>No. Optical fibers</td><td>Numerical aperture (NA)</td><td>Fiber diameter (mm)</td></tr></thead><tbody><tr><td>64</td><td>0.5</td><td>0.6</td></tr></tbody></table>

real optimization variables (as is typically implemented in conventional methods).

### C. Derived Absorption Coefficient Calculations

1) **Image Segmentation:** In order to calculate the fluence map and the DeACs inside a neonate head, we start by identifying the topological structure including the boundaries of the ECT, CSF, GM, WM, and the vessel regions. These regions are specified by applying a segmentation method, the nearest neighbor algorithm described in [74], on the IR-WRI derived SOS map. As published in the literature [75] and data available on the website (https://itis.swiss), the following reference SOS values are used: 2300 m/s, 1600 m/s, 1650 m/s, 1700 m/s, 1510 m/s for ECT, CSF, GM, WM, and vessel regions, respectively. With the segmentation algorithm, each pixel is assigned to the closest reference SOS. This results in identification of five segmented regions. As an example, if one pixel of the estimated SOS was 1655m/s, the reference value of 1650m/s is the closest value and therefore GM is assigned to that pixel. IR-WRI could theoretically estimate a DeAC for each of the segments described above. However, in reality, the tissue is optically heterogeneous. Determining only five DeACs, one for each of the five segmented head regions, is a rough estimation. To better estimate a real-life scenario of varying DeAC in the brain tissue, we further divide the main regions (i.e., ECT, CSF, GM, WM, and vessel regions) into subregions ($n_g$ = total number of subregions or segments) by overlaying a grid structure and using regional divisions to subdivide each gridded square (see Fig. 5(b)). We then determine DeAC for each subregion separately.

2) **Fluence Calculation:** Fluence refers to the amount of optical energy delivered to tissue per unit area. MCX, also known as Monte Carlo eXtreme, is a software specifically developed for simulating the transport of photons in turbid media using the Monte Carlo method. It specializes in time-resolved photon transport simulations [76]. In our study, we used MCX (with the input parameters in Table I) to calculate fluence. Iteratively, fluence was calculated based on scattering coefficients ($\mu_s$) and anisotropy values (g) extracted from literature (referenced as [77] and listed in Table II), and DeACs (generated by IR-WRI in each iteration) for various subregions. The simulation involved 640 million photon packets and utilized a grid size of 120 microns.

3) **DeAC Estimation:** Generally, the imaging problem is very sensitive to variations in the source term $\mathbf{b}$, which is the element-wise product of fluence and DeAC as defined in (5). However, there is a significant trade-off between $\mathbf{b}$ and the properties of the medium when both quantities need to be estimated. This ill-posedness prevents assignment of a degree of freedom to each sample of the medium for DeAC. Instead, a parsimonious parametrization of DeAC is preferred. Thus, in order to reduce the dimensionality of search space, we use
<table><caption>TABLE II
OPTICAL SCATTERING COEFFICIENT ($\mu_s$) AND ANISOTROPY FACTOR (g) OF THE NEONATE'S HEAD AT 800 NM [77]</caption><thead><tr><th>region</th><th>μ<sub>s</sub> (1/mm)</th><th>g</th></tr></thead><tbody><tr><td>ECT</td><td>12.57</td><td>0.89</td></tr><tr><td>CSF</td><td>3.0</td><td>0.9</td></tr><tr><td>GM</td><td>12.57</td><td>0.9</td></tr><tr><td>WM</td><td>12.57</td><td>0.9</td></tr><tr><td>Vessel</td><td>6.25</td><td>0.9</td></tr></tbody></table>

the segmented image (section III-C.1) based on the SOS map obtained from IR-WRI (section III-B). The PA source or initial pressure $\mathbf{b} \in \mathbb{C}^{N \times 1}$ can be decomposed as

$$ \mathbf{b} = \mathbf{FMx}, \quad (5) $$

where $\mathbf{F} \in \mathbb{R}^{N \times N}$ is the diagonal form of the fluence that was estimated in section III-C.2, and the columns of $\mathbf{M} \in \mathbb{R}^{N \times ng}$ contain the location of the segments according to the information extracted in section III-C.1. Also, $\mathbf{x} \in \mathbb{R}^{ng \times 1}$ is a vector of attenuation coefficient values for different segments. We direct interested readers to [78], [79], and [80] for further information on the discretized Helmholtz equation (1) and source model (5). Next, we demonstrate the IR-WRI method to extract values of DeACs in this parsimonious parametrization. Frequency-domain FWI with known model parameters $\mathbf{m}^0$ and unknown PA source $\mathbf{b}$ can be written as:

$$ \begin{array}{ll} \text{minimize} & \mathcal{R}(\mathbf{x}) \\ \text{subject to} & \left\{ \begin{array}{l} \mathbf{A}(\omega)\mathbf{u}(\omega) - \mathbf{FMx} = 0, \\ \mathbf{Pu}(\omega) = \mathbf{d}(\omega), \\ \mathbf{u}(\omega), \mathbf{x} \quad \tilde{\omega} = \tilde{\omega}_1, \tilde{\omega}_2, \dots, \tilde{\omega}_{n_f} \end{array} \right. \end{array} \quad (6) $$

where $\mathcal{R}(\mathbf{x})$ is an appropriate regularization function of $\mathbf{x}$, $\mathbf{u}(\omega) \in \mathbb{C}^{N \times 1}$ and $\mathbf{d}(\omega) \in \mathbb{C}^{N \times 1}$ denote the wavefield and the recorded data for frequency $\omega$, respectively, and $\tilde{\omega}_1, \tilde{\omega}_2, \dots, \tilde{\omega}_{n_f}$ indicate the frequency content of the PA signal. Also, $\mathbf{A}(\omega) \equiv \mathbf{A}(\mathbf{m}^0, \omega)$, and the dependency to the passive SOS model is removed (for the sake of simplicity). Determination of the optimum solution of (6) is extremely difficult and ill-posed, requiring sophisticated regularizations. In this study, we have used the total variation (TV) regularization [81] for $\mathbf{x}$, i.e. $\mathcal{R}(\mathbf{x}) = \| \nabla \mathbf{x} \|_1$ where $\nabla$ is the first-order difference operator. The augmented Lagrangian function associated with the problem in (6) is given by

$$ \begin{aligned} \tilde{\mathcal{L}}_A(\mathbf{u}(\omega), \mathbf{x}, \mathbf{v}(\omega), \mathbf{w}(\omega)) & = \|\nabla \mathbf{x}\|_1 - \sum_{\omega=\omega_1}^{\omega_{n_f}} \left\langle \mathbf{v}(\omega)^T, \mathbf{A}(\omega)\mathbf{u}(\omega) - \mathbf{FMx} \right\rangle \\ & \quad - \sum_{\omega=\omega_1}^{\omega_{n_f}} \left\langle \mathbf{w}(\omega)^T, \mathbf{Pu}(\omega) - \mathbf{d}(\omega) \right\rangle \\ & + \frac{\gamma_1}{2} \sum_{\omega=\omega_1}^{\omega_{n_f}} \|\mathbf{A}(\omega)\mathbf{u}(\omega) - \mathbf{FMx}\|_F^2 \\ & + \frac{\gamma_2}{2} \sum_{\omega=\omega_1}^{\omega_{n_f}} \|\mathbf{Pu}(\omega) - \mathbf{d}(\omega)\|_F^2, \end{aligned} \quad (7) $$

where the scalars $\gamma_1, \gamma_2 > 0$ are the penalty parameters assigned to the wave equation and the observation equation constraints, respectively, and $\mathbf{v}(\omega) \in \mathbb{C}^{N \times 1}$ and $\mathbf{w}(\omega) \in \mathbb{C}^{n_f \times 1}$ are the Lagrange multipliers. Beginning with the model $\mathbf{m}^0$ obtained in section III-B, $\mathbf{v}^0(\omega) = 0, \mathbf{w}^0(\omega) = 0 \forall \omega$ and initial value of absorption coefficients $\mathbf{x}^0$, ADMM solves the multivariate optimization problem (7) iteratively as follows:

$$ \mathbf{u}(\omega)^{k+1} = \arg \min_{\mathbf{u}(\omega)} \tilde{\mathcal{L}}_A(\mathbf{u}(\omega), \mathbf{x}^k, \mathbf{v}(\omega)^k, \mathbf{w}(\omega)^k), \quad (8a) $$

$$ \mathbf{x}^{k+1} = \arg \min_{\mathbf{x}} \tilde{\mathcal{L}}_A(\mathbf{u}(\omega)^{k+1}, \mathbf{x}, \mathbf{v}(\omega)^k, \mathbf{w}(\omega)^k), \quad (8b) $$

$$ \mathbf{v}(\omega)^{k+1} = \mathbf{v}(\omega)^k + \gamma_1 [\mathbf{FMx}^{k+1} - \mathbf{A}(\omega)\mathbf{u}(\omega)^{k+1}], \quad (8c) $$

$$ \mathbf{w}(\omega)^{k+1} = \mathbf{w}(\omega)^k + \gamma_2[\mathbf{d}(\omega) - \mathbf{Pu}(\omega)^{k+1}], \quad (8d) $$

where $\omega = \tilde{\omega}_1, \tilde{\omega}_2, \dots, \tilde{\omega}_{n_f}$.

Problem (8a) is another wavefield reconstruction problem (a problem that jointly satisfies the wave equation and the observation equation in a least-squares sense), whose solution can be obtained in a closed form similar to the one of (4a). The optimization problem (8b) reads

$$ \begin{aligned} \mathbf{x}^{k+1} = & \arg \min_{\mathbf{x}} \|\nabla \mathbf{x}\|_1 \\ & + \sum_{\omega=\tilde{\omega}_1}^{\tilde{\omega}_{n_f}} \left\langle \mathbf{w}(\omega)^k, \mathbf{A}(\omega)\mathbf{u}(\omega)^k - \mathbf{FMx} \right\rangle \\ & + \frac{\gamma_1}{2} \sum_{\omega=\tilde{\omega}_1}^{\tilde{\omega}_{n_f}} \|\mathbf{A}(\omega)\mathbf{u}(\omega)^k - \mathbf{FMx}\|_2^2. \end{aligned} \quad (9) $$

By adding and subtracting the term $\|\mathbf{w}(\omega)^k\|_2^2$ to (9), we have

$$ \begin{aligned} \mathbf{x}^{k+1} = & \arg \min_{\mathbf{x}} \|\nabla \mathbf{x}\|_1 + \frac{\gamma_1}{2} \sum_{\omega=\tilde{\omega}_1}^{\tilde{\omega}_{n_f}} \|\mathbf{A}(\omega)\mathbf{u}(\omega)^k - \mathbf{FMx} + \frac{1}{\gamma_1} \mathbf{w}(\omega)^k\|_2^2. \end{aligned} \quad (10) $$

where we have ignored the $-\|\mathbf{w}(\omega)^k\|_2^2$-term as it does not impact the optimization result. Equation (10) is a one dimension TV regularization and thus can be solved efficiently [82]. TV regularization has been used to control the large null space of the inversion. As previously mentioned, this algorithm requires an initial guess of DeACs ($\mathbf{x}^0$). It is estimated jointly with a set of monochromatic wavefields $\mathbf{u}^0(\omega)$. To this end, we may neglect the regularization term in (7) and solve the optimization problem jointly for $\mathbf{u}(\omega)$ and $\mathbf{x}$ for a few frequencies. The closed-form expression of this system for $r$ frequencies, i.e., $\tilde{\omega}_1, \dots, \tilde{\omega}_r$ can be written as:

$$ \begin{bmatrix} \sqrt{\gamma_1} \mathbf{A}(\tilde{\omega}_1) & 0 & 0 & 0 & -\mathbf{FM} \\ 0 & \sqrt{\gamma_1} \mathbf{A}(\tilde{\omega}_2) & 0 & 0 & -\mathbf{FM} \\ 0 & 0 & \dots & 0 & -\mathbf{FM} \\ 0 & 0 & 0 & \sqrt{\gamma_1} \mathbf{A}(\tilde{\omega}_r) & -\mathbf{FM} \\ \sqrt{\gamma_2} \mathbf{P} & 0 & 0 & 0 & 0 \\ 0 & \sqrt{\gamma_2} \mathbf{P} & 0 & 0 & 0 \\ 0 & 0 & \dots & 0 & \vdots \\ 0 & 0 & 0 & \sqrt{\gamma_2} \mathbf{P} & 0 \end{bmatrix} $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Block diagram of our qPAT method.</figcaption></figure>

$$ \begin{bmatrix} \mathbf{u}(\tilde{\omega}_1)^0 \\ \mathbf{u}(\tilde{\omega}_2)^0 \\ \vdots \\ \mathbf{u}(\tilde{\omega}_r)^0 \\ \mathbf{x}^0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ \sqrt{\gamma_2} \mathbf{d}(\tilde{\omega}_1) \\ \sqrt{\gamma_2} \mathbf{d}(\tilde{\omega}_2) \\ \vdots \\ \sqrt{\gamma_2} \mathbf{d}(\tilde{\omega}_r) \end{bmatrix} . \quad (11) $$

After solving (11) for $\mathbf{x}^0$, it, along with the extracted SOS and acoustic attenuation from section III-B can be used as inputs to subproblems (8) to extract accurate DeACs. In summary, our method involves the following steps:

* IR-WRI [68] is employed to find the SOS and acoustic attenuation maps based on US data.

* Pixels of the estimated SOS map are segmented into five major regions of the head.

* Finally, IR-WRI is employed to estimate DeAC from the MCX-updated fluence map and the acquired PA data in which the SOS and acoustic attenuation maps.

Fig. 3 shows the block diagram explaining the methodology we have used.

## IV. RESULTS

For all of the numerical computations, a computer with an Intel Xeon Platinum 8260 24C/35.75M/2.40G/165W as a CPU and 16x 64GB DDR4-2933 RDIMM as RAM was used. Also, the coefficients optimized nine-point stencil finite-difference method implemented with anti-lumped mass and perfectly matched layer (PML) absorbing boundary conditions [78] was used to perform forward modeling for all the tests.

### A. Estimation of SOS and Acoustic Attenuation

The acquisition device is comprised of 128 single element transducers of size 1 mm which collect data in the frequency band from 0.3 to 1 MHz. The direction of each transducer

<figure><img src="image_6.png" /><figcaption>Fig. 4. SOS and acoustic attenuation estimation by IR-WRI. (a) True SOS map, (b) initial guess for the SOS map, (c) estimated SOS map, (d-f), same as (a-c) but for acoustic attenuation, (g) cross sections of the SOS maps in a, b, and c at the location of the dashed line, (h) error graph for the SOS map (blue) and acoustic attenuation map (red).</figcaption></figure>

points towards the center of the ring. We use a set of frequencies, $f = [0.3 \ 0.4 \ 0.7 \ 0.8 \ 1]$ MHz, for mono-frequency IR-WRI followed by a multiscale frequency continuation strategy. We perform three passes through the frequencies, using the final model of one pass as the initial model for the next one. The starting and finishing frequencies of the three passes are [0.3, 0.4], [0.3, 0.4], and [0.3, 1] MHz, respectively. The stopping criterion with respect to iterations is IR-WRI must stop after a maximum of 30, 25, and 15 iterations for the first, second, and third passes, respectively. The hyperparameters ($\lambda_1$ and $\lambda_2$) are tuned based on a small fraction, e.g. 0.001, of the highest eigenvalue of $A(m^0, \omega_1)^{-T}P^{-T}PA(m^0, \omega_1)$. The readers can refer to appendix C of [71] for more details. Fig. 4(a) shows the true SOS model. Figs. 4(b) and 4(c) show the initial guess and estimated SOS maps after 185 iterations, respectively. Figs. 4(d)-4(f) are the same as Figs. 4(a)-4(c), but for acoustic attenuation. Cross-section profiles of the true, initial, and estimated SOS maps at the white dashed lines shown in Figs. 4(a)-4(c) are plotted in Fig. 4(g). The SOS percent error (SOS Error) defined as

$$ 100 \sum_{pixel} \frac{|\text{true SOS} - \text{estimated SOS}|}{|\text{true SOS}|}, \quad (12) $$
is shown in Fig. 4(h) at different iterations in blue, and the same plot for the estimated acoustic attenuation map (Acoustic Attenuation Error) is shown in red. The computational time for 185 iterations was 25 hours. As seen, IR-WRI is able to provide a highly accurate SOS map.

## B. DeAC Calculation

1) **Image Segmentation:** As explained in section III-C.1, a segmented image of the head was obtained using the SOS map estimated by IR-WRI. Fig. 5(a) shows the five regions (ECT, CSF, GM, WM, vessels) extracted from the SOS map. Fig. 5(b) shows the segmented model described in III-C.1. Results were used to calculate fluence maps as well as DeACs.

2) **Fluence Calculation:** The fluence map was calculated iteratively based on the specifications of the optical fibers and the optical properties of the neonatal head phantom at a wavelength of 800 nm (refer to Tables I and II in the Methods section) and the DeAC values were acquired from IR-WRI. Fig. 5(c) illustrates the calculated fluence map at the first iteration, while Fig. 5(e) represents the fluence map at the final iteration. It is evident from these figures that the optical scattering effect in biological tissue results in a significant decrease in the optical energy deposited in the central part of the head. To provide a clearer visualization of the fluence distribution, particularly in deeper areas, Fig. 5(d) and 5(f) display the logarithmic scale representation of the fluence map at the first and last iterations, respectively.

3) **DeAC Estimation of Subregions:** One hundred frequencies selected evenly over 0.3 to 1 MHz are used. To begin, the first 20 frequencies are used to solve (11) for extracting an initial DeAC ($x^0$). Then, all 100 frequencies are used simultaneously for solving subproblems (8) using the SOS and acoustic attenuation estimation maps shown in Fig. 4(c) and 4(f) and the fluence map estimation shown in Fig. 5(c). Tuning hyperparameters like $\gamma_1$ and $\gamma_2$ in (7) were challenging. We set them according to the following recipe, but of note, it is not the optimal one. The hyperparameter $\gamma_1$ is tuned as a small fraction, e.g. 0.001, of the highest eigenvalue of $M^T F^T F M$. Then $\gamma_2$ is tuned such that $\gamma = \gamma_1/\gamma_2$ will be a small fraction of the highest eigenvalue $\xi$ of $A(\tilde{\omega}_1)^{-T} P^T P A(\tilde{\omega}_1)$. We set $\gamma = 10^{-4}\xi$ for noiseless data and increase it when analyzing noisy data to prevent data overfitting. For example, we used $\gamma = 0.01\xi$ for the noisy data.

The true DeAC map used in the simulation is shown in Fig. 5(g). Next, an estimated DeAC map based on the SOS, acoustic attenuation and fluence maps that undergoes updates in each iteration is computed. Since the model we used in this study had only 5 regions with different DeACs, and real-case scenarios, because of the heterogeneity of tissue, will contain more regions, we decided to initially assume the possibility of the existence of up to 200 different DeACs and test to see if this process will correctly converge on the five DeACs in the true DeAC map. In order to quantify the accuracy of the reconstructed DeACs, we define DeAC percentage error.

$$100 \sum_{\text{pixel}} \frac{|\text{true DeAC} - \text{estimated DeAC}|}{|\text{true DeAC}|}, \quad (13)$$

<figure><img src="image_7.png" /><figcaption>Fig. 5. (a)-(b) Segmentation of the estimated SOS map with: (a) five regions, and (b) 200 regions (example regions are shown with r1, r2, ..., and r7). (c)-(f) Fluence map estimation in: (c) linear scale at 1st iteration, and (d) logarithmic scale at 1st iteration, (e) linear scale at 252th iteration, and (f) logarithmic scale at 252th iteration. (g)-(j) DeACs estimation by IR-WRI. (h) True DeAC, (i) DeACs at the 1st iteration with 200 regions, and (j) Error graph for the DeAC map.</figcaption></figure>

The DeAC percentage error for the full image is 9.41% and it is shown for each iteration in Fig. 5(j). The DeAC percentage error for each of the five regions is ECT=9.81%, CSF=9.62%, GM=9.01%, WM=10.40%, and vessel=10.25%. The results
<figure><img src="image_8.png" /><figcaption>Fig. 6. DeACs obtained from iterative time reversal (100 iterations). (a) With homogeneous SOS equal to 1500 m/s. (b) with SOS and acoustic attenuation maps extracted using IR-WRI (Figs. 4(c) and 4(f)).</figcaption></figure>

were obtained by running iterations until the rate of change in the error was < 1%: these conditions were satisfied beyond the 252^nd^ iteration. The results for the 1^st^ and 252^nd^ iterations are shown in Figs. 5(h)-5(i). The computational time for DeAC calculation for 450 iterations was ~49 hours.

For further assessment of the proposed method, we compare our results with the least squares (iterative) time reversal reconstruction algorithm [83], [84] which selected for comparison because it and the IR-WRI algorithms use the same wave equations and similar starting data. The time reversal algorithm was run for 100 iterations and stopped when the error did not decrease. We reconstructed the PA image and then divided it by the fluence map to generate a time-reversal (TR) DeAC map. Fig. 6(a) utilizes a homogeneous SOS = 1500 m/s and there is no adjustment for acoustic attenuation. Fig. 6(b) utilizes the SOS and acoustic attenuation maps derived previously using IR-WRI (Figs. 4(c) and 4(f)). The TR DeAC percentage error for the model shown in Fig. 6(b) is 71.98% and for each of the five regions is ECT=82.51%, CSF=86.3%, GM=80.6%, WM=70.3%, and vessel=40.2%.

Comparing the DeAC map from our proposed algorithm (Fig. 5(i)) and those extracted using the iterative time reversal algorithm (Fig. 6(b)), shows our proposed algorithm can recover the DeAC map nearly perfectly, while the iterative time reversal algorithm, using nearly the same starting information, poorly recovers the DeAC map, especially in deeper regions.

We next assess the robustness of our method against random noise in the data. We added different amounts of noise level to the raw PA data (noise level = 1 to 5 (Fig. 7(d)) where noise level 5 represents similar amplitude for noise and signal, and noise levels 1-4 represent 20% to 80% the amplitude of noise level 5. We next performed IR-WRI, and DeAC error varied from 11.6% to 56.1% (Fig. 7(c)). The DeAC maps for noise level 1 and noise level 5 are shown in Fig. 7(a) and 7(b).

V. DISCUSSION

Many commonly-used ray-based PA image reconstruction algorithms (e.g., universal back-projection [85], delay and sum [86] and multiply delay and sum [9]) assume that SOS and fluence maps are homogeneous across the entire imaging area and there is no acoustic attenuation. Therefore, to estimate DeACs, insufficient information is available, leading to large

<figure><img src="image_9.png" /><figcaption>Fig. 7. DeACs estimation for (a) noise level 1, (b) noise level 5, (c) DeAC error for different noise levels.</figcaption></figure>

error. Prior information from a different imaging modality can be used to decrease DeAC estimation error [30], [87], [88]. Because of the use of US transducers in PA data collection, US-PA dual modality is a natural choice to obtain such prior information. On the other hand, in wave-based methods, the full information contained in the recorded wavefields, i.e., amplitudes and phases, is considered without any assumption about level of heterogeneity [21], [22]. Therefore, these techniques can be effective for estimating accurate SOS and acoustic attenuation maps from US data. Several algorithms have been proposed for this very purpose [89], [90], [91]. Some limitations of these US methods are sensitivity to the selection of the initial model [90], massive computational burden [92], and suboptimal reconstruction of the acoustic attenuation map [93].

Several studies have attempted to estimate the SOS, the DeAC, or optical absorbance density [94], [95], [96], [97], [98], [99]. However, a few limitations associated with these methods (but more recently resolved) are: utilizing the Born approximation implies assuming SOS(x) = SOS₀(x) + SOSₚerturbation(x), where the magnitude of SOSₚerturbation/ SOS₀ is significantly smaller than 1 [97]. However, this assumption does not hold in our case due to the presence of the skull, where the SOS is approximately 2400 m/s while SOS₀ is around 1500 m/s. Another limitation is a requirement for incorporating prior knowledge about the geometry of the SOS map [99].

To estimate DeACs, several qPAT algorithms have been proposed [38], [40], [41], [43], [53], [100], [101]. Some historical limitations of these qPAT methods are: The assumption that the transducers have a point-like nature and lack directivity,
with a constrained aperture size [40]; uniform distribution of acoustic properties within the tissue, which may not always be the case in practice [38]; prior knowledge of the speed of sound in the medium is required [43]; use of diffusion approximation [100] to calculate fluence, which is not a correct approximation in the head due to the presence of cerebrospinal fluid (CSF) region, which acts as a non-scattering medium and poses a challenge for this approach. In fact, some older qPAT methods eschew fluence modeling entirely, and assume that the fluence is uniform throughout the tissue [7], [27], [29]. Moreover, recently, several research studies have made use of minimization-based methods [38], [39], [40], [101], [102], [103]. Nonetheless, the primary drawback associated with these newer approaches is their tendency to exhibit escalating computational complexity, ultimately resulting in substantial computation time requirements [104]. The time it takes for calculations depends on various factors, including the number of photon packets traced, the specifications of the hardware used, the size of the computational area, the quantity of transducers, and the complexity of the phantom being studied. For example, in [38], during the inversion process, 10⁷ photon packets were utilized to compute radiance and fluence, while 5 × 10⁶ photons were employed to calculate the corresponding adjoint quantities. On a high-end consumer GPU such as the NVIDIA GeForce Titan X Pascal, a single iteration of the inversion process, including the adjoint model, took approximately 84 seconds. On a different note, as described in [40], computations were performed on a computer equipped with two Intel(R) Xenon(R) Gold 6136 processors, possessing a total of 24 cores and 256 GB of memory. Each reconstruction step took around 150 minutes, with approximately 20 to 30 minutes allocated for tasks other than the primary computation.

Classical FWI can also be used for DeAC estimation. In this case the wave-equation constraint in (6) is solved for $u(\omega)$, i.e. $u(\omega) = A(\omega)^{-1}FMx$, and substituted into the second constraint in (1). The resulting optimization problem will be a regularized least-squares time-reversal. This optimization problem cannot be solved easily because of the presence of $A(\omega)^{-1}$, which is a very large and dense matrix that cannot be stored or directly inverted to extract $x$ [105]. Using the second IR-WRI algorithm to iteratively optimize parameters to include for fluence correction can be implemented with no initial guess or boundary assumption for the fluence, or preconceptions for the iteration process. This task would be infeasible for existing FWI-based methods because the nonlinearity of fluence compensation would set up another requirement for inverting a large and dense matrix.

In this paper we describe a qPAT framework to estimate DeACs. The requirements for this method are as follows:

* US data acquired by a ring of transducers when each transducer is activated as a transmitter in turn,

* Optical source locations and specifications, and

* PA data acquired by all transducers when the imaging target is irradiated by a laser beam.

The method operates by first calculating the SOS and acoustic attenuation maps using the acquired US data through the IR-WRI algorithm (with BM3D regularization).

It should be noted that the calculated acoustic attenuation map (error < 14%) is less accurate than the calculated SOS map (error < 5%). This is because the US recorded data is less responsive to acoustic attenuation compared to SOS. Subsequently, the IR-WRI algorithm is used again (with TV regularization) to solve a source estimation problem on the PA data (using the obtained SOS, acoustic attenuation and fluence maps). The PA problem is then solved by segmenting the SOS map to compensate for the trade-off between not knowing the pixelwise DeACs and the reduced size of the null space (11).

Some of the limitations of the proposed method are as follows. First, the segmentation method used is simplistic. Although the starting penalties used in conjunction (4) and (8) ultimately worked, these need to be optimized. The method could also be made more robust by making the DeAC estimation step less sensitive to the quality of the SOS map. Finally, computation time needs to be reduced either by more efficient coding or use of GPU.

## VI. CONCLUSION

In this paper we describe a qPAT framework to estimate DeACs. A key advantage of our qPAT method is its accuracy throughout a large imaging target due to the combined use of full wave inversion with fluence decay compensation delivered by the MCX toolbox: the estimated DeAC map is within 12% of the actual absorption map on which the simulation was based. Our results demonstrate the potential use of qPAT for transcranial imaging of neonatal brains which could be developed, for example, to detect conditions such as hemorrhage and hypoxia. Our method suffers from expensive computations; we plan to overcome this to some extent by implementing it on GPU.

## APPENDIX

### A. Standard Linear Solid (SLS) Attenuation Mechanism

Basically, m is a frequency-dependent non-linear mapping between SOS and attenuation factor at the reference frequency. [106] have compiled eight different models related to some empirical attenuation mechanisms. In this paper, we use SLS mechanism defined as

$$m = \frac{1}{\nu^2} \Re \left( \sqrt{\frac{1 + i \omega_r \tau_{\sigma}}{1 + i \omega_r \tau_{\varepsilon}}} \right)^{-2} \frac{1 + i \omega \tau_{\sigma}}{1 + i \omega \tau_{\varepsilon}}, \quad (14)$$

where $\omega_r$ is a reference frequency, $\nu$ is SOS, $\tau_\varepsilon$ and $\tau_\sigma$ are relaxation times related to the constants of the effective springs and dash-pot of the model [107] and defined as

$$\tau_{\varepsilon} = \frac{1}{\omega_r} \left( \sqrt{1 + \alpha^2} + \alpha \right), \quad (15)$$

$$\tau_{\sigma} = \frac{1}{\omega_r} \left( \sqrt{1 + \alpha^2} - \alpha \right), \quad (16)$$

where $\alpha$ is the acoustic attenuation. Also, by having complex valued m, one can extract acoustic attenuation as

$$\alpha = \frac{\omega^2 + \omega_r^2}{2\omega\omega_r \frac{\Re(1/m)}{\Im(1/m)}} \quad (17)$$
After extracting $\tau_{\varepsilon}$ and $\tau_{\sigma}$ using (15) and (16), SOS will be

$$v = \sqrt{\Re\left(\frac{1}{m}\right)\Re\left(\sqrt{\frac{1+i\omega_r\tau_{\varepsilon}}{1+i\omega_r\tau_{\sigma}}}\right) - 2\frac{1+\omega^2\tau_{\sigma}^2}{1+\omega^2\tau_{\sigma}\tau_{\varepsilon}}}. \quad (18)$$

REFERENCES

[1] H. Kye, Y. Song, T. Ninjbadgar, C. Kim, and J. Kim, "Whole-body photoacoustic imaging techniques for preclinical small animal studies," *Sensors*, vol. 22, no. 14, p. 5130, Jul. 2022.

[2] M. Xavierselvan and S. Mallidi, "LED-based functional photoacoustics—Portable and affordable solution for preclinical cancer imaging," in *LED-Based Photoacoustic Imaging: From Bench to Bedside*. Singapore: Springer, 2020, pp. 303–319.

[3] A. B. E. Attia et al., "A review of clinical photoacoustic imaging: Current and future trends," *Photoacoustics*, vol. 16, Dec. 2019, Art. no. 100144.

[4] L. Lin and L. V. Wang, "The emerging role of photoacoustic imaging in clinical oncology," *Nature Rev. Clin. Oncol.*, vol. 19, no. 6, pp. 365–384, Jun. 2022.

[5] T. J. Allen and P. C. Beard, "High power visible light emitting diodes as pulsed excitation sources for biomedical photoacoustics," *Biomed. Opt. Exp.*, vol. 7, no. 4, pp. 1260–1270, 2016.

[6] S. Gupta, M. R. Chatni, A. L. Rao, V. I. Vullev, L. V. Wang, and B. Anvari, "Virus-mimicking nano-constructs as a contrast agent for near infrared photoacoustic imaging," *Nanoscale*, vol. 5, no. 5, pp. 1772–1776, 2013.

[7] C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, "Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media," *IEEE Trans. Med. Imag.*, vol. 32, no. 6, pp. 1097–1110, Jun. 2013.

[8] K. P. Köstli and P. C. Beard, "Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response," *Appl. Opt.*, vol. 42, no. 10, pp. 1899–1908, Apr. 2003.

[9] M. Mozaffarzadeh, A. Mahlooijar, M. Orooji, S. Adabi, and M. Nasirivavanaki, "Double-stage delay multiply and sum beamforming algorithm: Application to linear-array photoacoustic imaging," *IEEE Trans. Biomed. Eng.*, vol. 65, no. 1, pp. 31–42, Jan. 2018.

[10] B. E. Treeby, E. Z. Zhang, and B. T. Cox, "Photoacoustic tomography in absorbing acoustic media using time reversal," *Inverse Problems*, vol. 26, no. 11, Nov. 2010, Art. no. 115003.

[11] M. Xu and L. V. Wang, "Universal back-projection algorithm for photoacoustic computed tomography," *Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.*, vol. 71, no. 1, Jan. 2005, Art. no. 016706.

[12] Y. Xu and L. V. Wang, "Effects of acoustic heterogeneity in breast thermoacoustic tomography," *IEEE Trans. Control*, vol. 50, no. 9, pp. 1134–1146, Sep. 2003.

[13] X. L. Deán-Ben, V. Ntziachristos, and D. Razansky, "Statistical optical-acoustic image reconstruction using a-priori knowledge on the location of acoustic distortions," *Appl. Phys. Lett.*, vol. 98, no. 17, Apr. 2011, Art. no. 171110.

[14] B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard, "Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach," *J. Biomed. Opt.*, vol. 16, no. 9, pp. 090501-1-090501-3, 2011.

[15] T. Lin, "Correlations of sound speed with tissue constituents in normal and diffuse liver disease," *Ultrason. Imag.*, vol. 9, no. 1, pp. 29–40, Jan. 1987.

[16] K. Zell, J. I. Sperl, M. W. Vogel, R. Niessner, and C. Haisch, "Acoustic properties of selected tissue phantom materials for ultrasound imaging," *Phys. Med. Biol.*, vol. 52, no. 20, pp. N475–N484, Oct. 2007.

[17] J. Jose, R. G. H. Willemink, W. Steenbergen, C. H. Slump, T. G. van Leeuwen, and S. Manohar, "Speed-of-sound compensated photoacoustic tomography for accurate imaging," *Med. Phys.*, vol. 39, no. 12, pp. 7262–7271, Dec. 2012.

[18] J. Yin, J. He, C. Tao, and X. Liu, "Enhancement of photoacoustic tomography of acoustically inhomogeneous tissue by utilizing a memory effect," *Opt. Exp.*, vol. 28, no. 8, pp. 10806–10817, 2020.

[19] C. Zhang and Y. Wang, "A reconstruction algorithm for thermoacoustic tomography with compensation for acoustic speed heterogeneity," *Phys. Med. Biol.*, vol. 53, no. 18, pp. 4971–4982, Sep. 2008.

[20] Z. Sun and L.-S. Sun, "Simultaneous reconstruction of optical absorption property and speed of sound in intravascular photoacoustic tomography," *Inverse Problems Sci. Eng.*, vol. 29, no. 12, pp. 1764–1788, Dec. 2021.

[21] R. Gerhard Pratt, L. Huang, N. Duric, and P. Littrup, "Sound-speed and attenuation imaging of breast tissue using waveform tomography of transmission ultrasound data," *Proc. SPIE*, vol. 6510, pp. 1523–1534, Mar. 2007.

[22] R. G. Pratt, C. Shin, and G. J. Hick, "Gauss–Newton and full Newton methods in frequency–space seismic waveform inversion," *Geophys. J. Int.*, vol. 133, no. 2, pp. 341–362, May 1998.

[23] M. Warner et al., "Anisotropic 3D full-waveform inversion," *Geo-physics*, vol. 78, no. 2, pp. R59–R80, 2013.

[24] L. Guasch, O. C. Agudo, M.-X. Tang, P. Nachev, and M. Warner, "Full-waveform inversion imaging of the human brain," *npj Digit. Med.*, vol. 3, no. 1, p. 28, Mar. 2020.

[25] C. Huang, K. Wang, L. Nie, L. V. Wang, and M. A. Anastasio, "Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media," *IEEE Trans. Med. Imag.*, vol. 32, no. 6, pp. 1097–1110, Jun. 2013.

[26] G. Zangerl, M. Haltmeier, L. V. Nguyen, and R. Nuster, "Full field inversion in photoacoustic tomography with variable sound speed," *Appl. Sci.*, vol. 9, no. 8, p. 1563, Apr. 2019.

[27] J. Qian, P. Stefanov, G. Uhlmann, and H. Zhao, "An efficient Neumann series–based algorithm for thermoacoustic and photoacoustic tomography with variable sound speed," *SIAM J. Imag. Sci.*, vol. 4, no. 3, pp. 850–883, Jan. 2011.

[28] L. Yao and H. Jiang, "Enhancing finite element-based photoacoustic tomography using total variation minimization," *Appl. Opt.*, vol. 50, no. 25, pp. 5031–5041, Sep. 2011.

[29] Z. Yuan and H. Jiang, "Three-dimensional finite-element-based photoacoustic tomography: Reconstruction algorithm and simulations," *Med. Phys.*, vol. 34, no. 2, pp. 538–546, Feb. 2007.

[30] S. Mahmoodkalayeh et al., "Improving vascular imaging with co-planar mutually guided photoacoustic and diffuse optical tomography: A simulation study," *Biomed. Opt. Exp.*, vol. 11, no. 8, pp. 4333–4347, 2020.

[31] S. M. Ranjbaran, K. Kratkiewicz, R. Manwar, H. Fallah, M. Hajimahmoodzadeh, and M. R. N. Avanaki, "Numerical study on tumor volume detection in the human brain using photoacoustic tomography," *Proc. SPIE*, vol. 10878, pp. 508–515, Mar. 2019.

[32] L. Zhao, M. Yang, Y. Jiang, and C. Li, "Optical fluence compensation for handheld photoacoustic probe: An in vivo human study case," *J. Innov. Opt. Health Sci.*, vol. 10, no. 4, Jul. 2017, Art. no. 1740002.

[33] F. Guerra and D. S. Dumani, "An iterative method of light fluence distribution estimation for quantitative photoacoustic imaging," *Proc. SPIE*, vol. 11642, pp. 243–250, Mar. 2021.

[34] M. Kirillin, V. Perekatova, I. Turchin, and P. Subochev, "Fluence compensation in raster-scan optoacoustic angiography," *Photoacoustics*, vol. 8, pp. 59–67, Dec. 2017.

[35] T. Tarvainen, A. Pulkkinen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, "Bayesian image reconstruction in quantitative photoacoustic tomography," *IEEE Trans. Med. Imag.*, vol. 32, no. 12, pp. 2287–2298, Dec. 2013.

[36] H. Gao, S. Osher, and H. Zhao, "Quantitative photoacoustic tomography," in *Mathematical Modeling in Biomedical Imaging II: Optical, Ultrasound, and Opto-Acoustic Tomographies*. Berlin, Germany: Springer, 2011, pp. 131–158.

[37] T. Tarvainen, B. T. Cox, J. P. Kaipio, and S. R. Arridge, "Reconstructing absorption and scattering distributions in quantitative photoacoustic tomography," *Inverse Problems*, vol. 28, no. 8, Aug. 2012, Art. no. 084009.

[38] J. Buchmann, B. A. Kaplan, S. Powell, S. Prohaska, and J. Laufer, "Three-dimensional quantitative photoacoustic tomography using an adjoint radiance Monte Carlo model and gradient descent," *J. Biomed. Opt.*, vol. 24, no. 6, p. 066001, 2019.

[39] S. Zheng, H. Yingsa, S. Meichen, and M. Qi, "Quantitative photoacoustic tomography with light fluence compensation based on radiance Monte Carlo model," *Phys. Med. Biol.*, vol. 68, no. 6, 2023, Art. no. 065009.

[40] A. A. Leino, T. Lunttila, M. Mozumder, A. Pulkkinen, and T. Tarvainen, "Perturbation Monte Carlo method for quantitative photoacoustic tomography," *IEEE Trans. Med. Imag.*, vol. 39, no. 10, pp. 2985–2995, Oct. 2020.
[41] N. Hänninen, A. Pulkkinen, S. Arridge, and T. Tarvainen, "Adaptive stochastic Gauss-Newton method with optical Monte Carlo for quantitative photoacoustic tomography," *J. Biomed. Opt.*, vol. 27, no. 08, Apr. 2022, Art. no. 083013.

[42] T. Saratoon, T. Tarvainen, B. T. Cox, and S. R. Arridge, "A gradient-based method for quantitative photoacoustic tomography using the radiative transfer equation," *Inverse Problems*, vol. 29, no. 7, Jul. 2013, Art. no. 075006.

[43] M. Haltmeier, L. Neumann, and S. Rabanser, "Single-stage reconstruction algorithm for quantitative photoacoustic tomography," *Inverse Problems*, vol. 31, no. 6, Jun. 2015, Art. no. 065005.

[44] K. Kratkiewicz et al., "Development of a stationary 3D photoacoustic imaging system using sparse single-element transducers: Phantom study," *Appl. Sci.*, vol. 9, no. 21, p. 4505, Oct. 2019.

[45] S. Mahmoodkalayeh et al., "A simulation study on co-planar mutually guided photoacoustic and diffuse optical tomography," *Proc. SPIE*, vol. 11642, Mar. 2021, Art. no. 116423X.

[46] Y. Tang and J. Yao, "3D Monte Carlo simulation of light distribution in mouse brain in quantitative photoacoustic computed tomography," *Quant. Imag. Med. Surg.*, vol. 11, no. 3, pp. 1046-1059, Mar. 2020.

[47] A. Pulkkinen, B. T. Cox, S. R. Arridge, J. P. Kaipio, and T. Tarvainen, "A Bayesian approach to spectral quantitative photoacoustic tomography," *Inverse Problems*, vol. 30, no. 6, Jun. 2014, Art. no. 065012.

[48] A. Javaherian and S. Holman, "Direct quantitative photoacoustic tomography for realistic acoustic media," *Inverse Problems*, vol. 35, no. 8, Aug. 2019, Art. no. 084004.

[49] A. Hannukainen, N. Hyvönen, H. Majander, and T. Tarvainen, "Efficient inclusion of total variation type priors in quantitative photoacoustic tomography," *SIAM J. Imag. Sci.*, vol. 9, no. 3, pp. 1132-1153, Jan. 2016.

[50] T. Tarvainen, M. Vauhkonen, V. Kolehmainen, and J. P. Kaipio, "Finite element model for the coupled radiative transfer equation and diffusion approximation," *Int. J. Numer. Methods Eng.*, vol. 65, no. 3, pp. 383-405, Jan. 2006.

[51] G. Bal and K. Ren, "On multi-spectral quantitative photoacoustic tomography in diffusive regime," *Inverse Problems*, vol. 28, no. 2, Feb. 2012, Art. no. 025010.

[52] B. Cox, J. G. Laufer, S. R. Arridge, and P. C. Beard, "Quantitative spectroscopic photoacoustic imaging: A review," *J. Biomed. Opt.*, vol. 17, no. 6, 2012, Art. no. 061202.

[53] A. Pattyn, Z. Mumm, N. Alijabbari, N. Duric, M. A. Anastasio, and M. Mehrmohammadi, "Model-based optical and acoustical compensation for photoacoustic tomography of heterogeneous mediums," *Photoacoustics*, vol. 23, Sep. 2021, Art. no. 100275.

[54] B. T. Cox, S. R. Arridge, K. 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," *Appl. Opt.*, vol. 45, no. 8, pp. 1866-1875, Mar. 2006.

[55] H. Gao, J. Feng, and L. Song, "Limited-view multi-source quantitative photoacoustic tomography," *Inverse Problems*, vol. 31, no. 6, Jun. 2015, Art. no. 065004.

[56] G. Bal and K. Ren, "Multi-source quantitative photoacoustic tomography in a diffusive regime," *Inverse Problems*, vol. 27, no. 7, Jul. 2011, Art. no. 075003.

[57] X. Zhou, N. Akhlaghi, K. A. Wear, B. S. Garra, T. J. Pfefer, and W. C. Vogt, "Evaluation of fluence correction algorithms in multispectral photoacoustic imaging," *Photoacoustics*, vol. 19, Sep. 2020, Art. no. 100181.

[58] S. Brigadoi, P. Aljabar, M. Kuklisova-Murgasova, S. R. Arridge, and R. J. Cooper, "A 4D neonatal head model for diffuse optical imaging of pre-term to term infants," *NeuroImage*, vol. 100, pp. 385-394, Oct. 2014.

[59] S. Brigadoi and R. J. Cooper, "How short is short? Optimum source-detector distance for short-separation channels in functional near-infrared spectroscopy," *Neurophotonics*, vol. 2, no. 2, May 2015, Art. no. 025005.

[60] L. Wu, Y. Lin, and T. Li, "Effect of human brain edema on light propagation: A Monte Carlo modeling based on the visible Chinese human dataset," *IEEE Photon. J.*, vol. 9, no. 5, pp. 1-10, Oct. 2017.

[61] Q. Zhao, L. Ji, F. Shi, and T. Jiang, "Distinguish activations on sensorimotor cortex using high-resolution diffuse optical tomography," *Proc. SPIE*, vol. 6534, pp. 46-53, May 2007.

[62] H. S. Aghamiry, A. Gholami, and S. Operto, "Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: Application to large contrast media," *Geophys. J. Int.*, pp. 855-872, Apr. 2019.

[63] A. Rieder, "An all-at-once approach to full waveform inversion in the viscoelastic regime," *Math. Methods Appl. Sci.*, vol. 44, no. 8, pp. 6376-6388, May 2021.

[64] J. Nocedal and S. J. Wright, *Numerical Optimization* (Springer Series in Operations Research and Financial Engineering). New York, NY, USA: Springer, 1999.

[65] B. Kaltenbacher, "All-at-once versus reduced iterative methods for time dependent inverse problems," *Inverse Problems*, vol. 33, no. 6, Jun. 2017, Art. no. 064002.

[66] A. Tarantola, "Theoretical background for the inversion of seismic waveforms, including elasticity and attenuation," in *Scattering and Attenuations of Seismic Waves, Part I*. Basel, Switzerland: Birkhäuser, 1988, pp. 365-399.

[67] J. Virieux and S. Operto, "An overview of full-waveform inversion in exploration geophysics," *Geophysics*, vol. 74, no. 6, pp. WCC1-WCC26, Nov. 2009.

[68] H. S. Aghamiry, A. Gholami, and S. Operto, "Improving full-waveform inversion by wavefield reconstruction with the alternating direction method of multipliers," *Geophysics*, vol. 84, no. 1, pp. R125-R148, Jan. 2019.

[69] H. S. Aghamiry, A. Gholami, and S. Operto, "Complex-valued imaging with total variation regularization: An application to full-waveform inversion in visco-acoustic media," *SIAM J. Imag. Sci.*, vol. 14, no. 1, pp. 58-91, Jan. 2021.

[70] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Found. Trends Mach. Learn.*, vol. 3, no. 1, pp. 1-122, 2010.

[71] H. S. Aghamiry, A. Gholami, and S. Operto, "Compound regularization of full-waveform inversion for imaging piecewise media," *IEEE Trans. Geosci. Remote Sens.*, vol. 58, no. 2, pp. 1192-1204, Feb. 2020.

[72] M. Benning and M. Burger, "Modern regularization methods for inverse problems," *Acta Numerica*, vol. 27, pp. 1-111, May 2018.

[73] K. Dabov, "Video denoising by sparse 3D transform-domain collaborative filtering," in *Proc. Eur. Signal Process. Conf.*, 2007, pp. 145-149.

[74] B. Li, Y. W. Chen, and Y. Q. Chen, "The nearest neighbor algorithm of local probability centers," *IEEE Trans. Syst. Man, Cybern. B, Cybern.*, vol. 38, no. 1, pp. 141-154, Feb. 2008.

[75] L. Palaniappan and V. Velusamy, "Ultrasonic study of human cerebrospinal fluid," *Indian J. Pure Appl. Phys.*, vol. 42, pp. 591-594, Aug. 2004.

[76] Q. Fang and D. A. Boas, "Monte Carlo Simulation of photon migration in 3D turbid media accelerated by graphics processing units," *Opt. Exp.*, vol. 17, no. 22, pp. 20178-20190, 2009.

[77] S. L. Jacques, "Optical properties of biological tissues: A review," *Phys. Med. Biol.*, vol. 58, no. 11, pp. R37-R61, Jun. 2013.

[78] Z. Chen, T. Wu, and H. Yang, "An optimal 25-point finite difference scheme for the Helmholtz equation with PML," *J. Comput. Appl. Math.*, vol. 236, no. 6, pp. 1240-1258, Oct. 2011.

[79] J. O. Robertsson, J. O. Blanch, K. Nihei, and J. Tromp, *Numerical Modeling of Seismic Wave Propagation: Gridded Two-Way Wave Equation Methods*. Tulsa, Oklahoma: Society of Exploration Geophysicists, Jan. 2012.

[80] H. S. Aghamiry, A. Gholami, L. Combe, and S. Operto, "Accurate 3D frequency-domain seismic wave modeling with the wavelength-adaptive 27-point finite-difference stencil: A tool for full-waveform inversion," *Geophysics*, vol. 87, no. 3, pp. R305-R324, May 2022.

[81] L. I. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," *Phys. D, Nonlinear Phenomena*, vol. 60, nos. 1-4, pp. 259-268, Nov. 1992.

[82] T. Goldstein and S. Osher, "The split Bregman method for L1-regularized problems," *SIAM J. Imag. Sci.*, vol. 2, no. 2, pp. 323-343, Jan. 2009.

[83] G. Montaldo, M. Taner, and M. Fink, "Revisiting iterative time reversal processing: Application to detection of multiple targets," *J. Acoust. Soc. Amer.*, vol. 115, no. 2, pp. 776-784, Feb. 2004.

[84] B. T. Taddese, T. M. Antonsen, E. Ott, and S. M. Anlage, "Iterative time reversal with tunable convergence," 2011, arXiv:1107.1425.

[85] L. V. Wang, *Photoacoustic Imaging and Spectroscopy*. Boca Raton, FL, USA: CRC Press, 2017.

[86] P. Warbal, M. Pramanik, and R. K. Saha, "A robust modified delay-and-sum algorithm for photoacoustic tomography imaging with apodized transducers," in *Proc. Eur. Conf. Biomed. Opt.*, 2019, pp. 1107713-1-107713-8, Paper 11077_38.

[87] V. Sorgato et al., "Validation of optical properties quantification with a dual-step technique for biological tissue analysis," *J. Biomed. Opt.*, vol. 23, no. 9, 2018, Art. no. 096002.
[88] S. Zhang et al., "MRI information-based correction and restoration of photoacoustic tomography," *IEEE Trans. Med. Imag.*, vol. 41, no. 9, pp. 2543-2555, Sep. 2022.

[89] D. M. Brandner, X. Cai, J. Foiret, K. W. Ferrara, and B. G. Zagar, "Estimation of tissue attenuation from ultrasonic B-mode images—Spectral-log-difference and method-of-moments algorithms compared," *Sensors*, vol. 21, no. 7, p. 2548, Apr. 2021.

[90] J. Fincke, X. Zhang, B. Shin, G. Ely, and B. W. Anthony, "Quantitative sound speed imaging of cortical bone and soft tissue: Results from observational data sets," *IEEE Trans. Med. Imag.*, vol. 41, no. 3, pp. 502-514, Mar. 2022.

[91] W. A. Simson, M. Paschali, V. Sideri-Lampretsa, N. Navab, and J. J. Dahl, "Investigating pulse-echo sound speed estimation in breast ultrasound with deep learning," 2023, arXiv:2302.03064.

[92] C. Holmes, B. W. Drinkwater, and P. D. Wilcox, "Post-processing of the full matrix of ultrasonic transmit-receive array data for non-destructive evaluation," *NDT & E Int.*, vol. 38, no. 8, pp. 701-711, Dec. 2005.

[93] M. Pérez-Liva, J. L. Herraiz, J. M. Udfás, E. Miller, B. T. Cox, and B. E. Treeby, "Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion," *J. Acoust. Soc. Amer.*, vol. 141, no. 3, pp. 1595-1604, Mar. 2017.

[94] C. Cai, X. Wang, K. Si, J. Qian, J. Luo, and C. Ma, "Feature coupling photoacoustic computed tomography for joint reconstruction of initial pressure and sound speed in vivo," *Biomed. Opt. Exp.*, vol. 10, no. 7, pp. 3447-3462, 2019.

[95] C. Huang, K. Wang, R. W. Schoonover, L. V. Wang, and M. A. Anastasio, "Joint reconstruction of absorbed optical energy density and sound speed distributions in photoacoustic computed tomography: A numerical investigation," *IEEE Trans. Comput. Imag.*, vol. 2, no. 2, pp. 136-149, Jun. 2016.

[96] T. P. Matthews and M. A. Anastasio, "Joint reconstruction of the initial pressure and speed of sound distributions from combined photoacoustic and ultrasound tomography measurements," *Inverse Problems*, vol. 33, no. 12, Dec. 2017, Art. no. 124002.

[97] T. Ding, K. Ren, and S. Vallélian, "A one-step reconstruction algorithm for quantitative photoacoustic imaging," *Inverse Problems*, vol. 31, no. 9, Sep. 2015, Art. no. 095005.

[98] Z. Yuan, Q. Zhang, and H. Jiang, "Simultaneous reconstruction of acoustic and optical properties of heterogeneous media by quantitative photoacoustic tomography," *Opt. Exp.*, vol. 14, no. 15, pp. 6749-6754, 2006.

[99] J. Zhang, K. Wang, Y. Yang, and M. A. Anastasio, "Simultaneous reconstruction of speed-of-sound and optical absorption properties in photoacoustic tomography via a time-domain iterative algorithm," in *Proc. SPIE*, vol. 6856, pp. 427-434, Feb. 2008.

[100] A. Javaherian and S. Holman, "Direct quantitative photoacoustic tomography for realistic acoustic media," *Inverse Problems*, vol. 35, no. 8, Aug. 2019, Art. no. 084004.

[101] C. M. Macdonald, S. Arridge, and S. Powell, "Efficient inversion strategies for estimating optical properties with Monte Carlo radiative transport models," *J. Biomed. Opt.*, vol. 25, no. 8, Aug. 2020, Art. no. 085002.

[102] R. Hochuli, S. Powell, S. Arridge, and B. Cox, "Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance," *J. Biomed. Opt.*, vol. 21, no. 12, Dec. 2016, Art. no. 126004.

[103] A. Capart, J. Wojak, R. Allais, M. Ghiss, O. Boiron, and A. Da Silva, "Quantitative photoacoustic reconstruction of the optical properties of intervertebral discs using a gradient descent scheme," *Photonics*, vol. 9, no. 2, p. 116, Feb. 2022.

[104] Z. Wang, W. Tao, and H. Zhao, "The optical inverse problem in quantitative photoacoustic tomography: A review," *Photonics*, vol. 10, no. 5, p. 487, Apr. 2023.

[105] G. Montaldo, M. Tanter, and M. Fink, "Real time inverse filter focusing through iterative time reversal," *J. Acoust. Soc. Amer.*, vol. 115, no. 2, pp. 768-775, Feb. 2004.

[106] B. Ursin and T. Toverud, "Comparison of seismic dispersion and attenuation models," *Stud. Geophys. Geodaetica*, vol. 46, no. 2, pp. 293-320, 2002.

[107] T. Zhu, J. M. Carcione, and J. M. Harris, "Approximating constant-Q seismic propagation in the time domain," *Geophys. Prospecting*, vol. 61, no. 5, pp. 931-940, 2013.