# Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study

Charlotte Constans, Ivana Falco, Bastien Arnal

► To cite this version:

Charlotte Constans, Ivana Falco, Bastien Arnal. Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study. 2025. hal-04728441v2

HAL Id: hal-04728441

https://hal.science/hal-04728441v2

Preprint submitted on 3 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Mitigating the unicity problem in quantitative photoacoustic imaging using ultrasound blood flow imaging, a simulation study

Charlotte Constans¹, Ivana Falco¹, and Bastien Arnal¹

¹ Univ. Grenoble Alpes, CNRS, LIPhy, Grenoble 38058, FRANCE

E-mail: bastien.arnal@univ-grenoble-alpes.fr

January 2025

**Abstract.** *Objective.* Quantitative photoacoustic (PA) imaging aims at recovering optical properties in optically turbid media. Here, we aim to recover uniqueness of the solution in absorption and scattering coefficient in a single illumination and wavelength setting. *Approach.* By considering a numerical tissue model composed of optically contrasted regions surrounded by blood vessels, we combine PA imaging with ultrasound power Doppler (UPD). We show that quantitative absorption coefficients inside the blood vessels can be derived from UPD, serving as a prior for a radiance Monte Carlo inverse model. *Main results.* We show that using this vascular imaging prior mitigates the non-uniqueness issue, providing quantitative reconstructions of absorption and scattering even in the absence of background or blood absorption prior information. We generalize to unresolved capillary beds by degrading the resolution of UPD and PA. *Significance.* Since dual PA-ultrasound systems are already widespread, this approach could become an essential brick of quantitative PA imaging.

*Keywords:* Inverse problems; Quantitative Photoacoustic imaging; Blood flow imaging; ultrasound Power Doppler; Monte Carlo radiance Submitted to: *Phys. Med. Biol.*

## 1. Introduction

Photoacoustic (PA) imaging is an emerging biomedical imaging technique that has proven able to reveal optical contrast at depths in soft tissues [1, 2, 3]. This technique has the potential to offer real-time imaging [4] at a high spatial resolution and penetration depth of endogenous, exogenous or genetically-encoded absorbers [5]. PA imaging theoretically offers quantitative chromophore concentration imaging [6]. This applies to ranges of applications going from oxygenation imaging [7, 8, 9] to molecular imaging [9, 10]. For this task, one needs to use spectral features from images obtained in a PA spectroscopic acquisition and operate spectral unmixing of the different chromophores
[6]. However, PA sources are proportional to the absorbed energy density, equal to the product between optical absorption coefficient and optical fluence $\mu_a(\lambda, \mathbf{r}) \times \phi(\lambda, \mathbf{r})$. As optical transport is wavelength-dependant, an effect referred to as spectral coloring [11] requires $\phi(\lambda, \mathbf{r})$ to be computed for each wavelength to apply a fluence correction before any demixing algorithm. When the transport or scattering properties of the background medium can be considered homogeneous, techniques have been proposed to account for the fluence distribution in the images [12, 13, 14, 15]. When unknown heterogeneous media are considered, a light propagation model has to be employed [16] within an inverse algorithm. From a PA image, the basic goal of the inverse model is to recreate a digital twin by finding the optimal optical properties throughout the medium. Ultimately, this allows to accurately recover decoupled absorption and scattering properties, providing further quantitative spectral analysis. Learning techniques based on simulated examples appear as simpler ways than the full inverse problem [8, 17, 18, 19, 20, 21] but they might lack generality since they rely on an ensemble of simulated examples and usually consider uniform background scattering properties.

As PA imaging is both optical and acoustical, complete inverse problems involve two distinct parts [22, 23]. In this work, we focus on the optical inversion part, considering the acoustic reconstruction ideal. Optical transport is governed by the Radiative transport equation (RTE) which can be in some cases simplified to the diffusion approximation [24, 12]. As diffusion approximation can be applied only a few mean free paths away from the illuminated surface, modelling of the RTE is preferable for more accuracy [16]. RTE can be modelled using Finite-element methods [25] but the advent of parallel GPU computing makes Monte Carlo-based methods more attractive for their computational speed [26, 27, 28]. A radiance Monte Carlo inverse method using the directional information of photons was introduced to inverse the RTE without the diffusion approximation [29, 30, 31].

In quantitative PA imaging, inverse problems aim at recovering the optical scattering $\mu_s(\mathbf{r})$ and absorption $\mu_a(\mathbf{r})$ coefficients, from measured data. It has been demonstrated that the solution to the photoacoustic inverse problem is not unique when using a single wavelength and illumination. The reason why is when the calculated absorbed energy $H$ differs from the measured one, it is impossible to know whether it comes from an error in the absorption or scattering distribution, which are both responsible for light attenuation. This means that multiple couples $(\mu_s(\mathbf{r}), \mu_a(\mathbf{r}))$ satisfy the inversion [6]. To solve the non-uniquity issue, researchers proposed to use multiple illumination patterns [32, 33, 25], but this requires dedicated equipment and illumination calibrations which might not be adequate for medical applications. Another approach would be to constrain the problem by having strong priors. Using multiple wavelengths, priors on the chromophore spectral properties can lead to a unique solution [6] in both $\mu_s(\mathbf{r})$ and $\mu_a(\mathbf{r})$. There are a few experimental results of multispectral inversions but some assumptions are made on $\mu_s$ and inversions only estimate $\mu_a$ [30, 31]. Clearly, the simultaneous inversion of $\mu_s(\mathbf{r})$ and $\mu_a(\mathbf{r})$ is considered as a bottleneck when the
medium is heterogeneous with no prior information. In this issue, 2D simulated data input is considered. Our aim is to propose a method to facilitate the recovery of both scattering and absorption coefficients with single illumination and wavelength. To this end, we will rely on blood vessels.

Blood vessels are ubiquitous in the body. Contained within these vessels, hemoglobin is also the main optical absorber present in tissues for near-infrared light. In this paper, we propose to use a multimodal approach with the coupling of PA imaging and another modality offering blood vessels imaging to gain prior information. This modality is supposed to be of a distinct physical nature and to be insensitive to the optical transport phenomenon. A good candidate as a coupled modality is ultrasound imaging since simultaneous PA and US imaging is common [34, 35, 36]. Ultrasound-informed PA reconstructions were proposed using ultrasound-driven anatomical information [37, 38]. While this approach is relevant, it does not exploit quantitative information that can be given by ultrasound. One interesting modality is Ultrasound power Doppler (UPD) [39] that enables to image specifically the vasculature. Recent developments improved significantly the sensitivity of the technique that is now applied for functional imaging in the brain [40, 41]. State-of-the-art UPD was recently combined with PA imaging [42, 43]. Different high-resolution instrumentation able to image blood vessels (ex: CT angiography, MR-angiography scanners or optical imaging) may also be coupled to PA imaging but it seems adequate to consider ultrasound-based techniques for their resolution, penetration and as the same equipment can be used for dual PA and ultrasound imaging. Considering blood vessels as a support of valuable information, we propose to use UPD to address the non-unicity problem with a single-illumination and wavelength PA+US system.

UPD consists in evaluating the fluctuations caused by blood flow. It consists in integrating the whole Doppler spectrum. Consequently, UPD measures the agitation of blood and is insensitive to the flow speed. More precisely, it quantifies blood volume (the fraction of blood within one imaging pixel/voxel) under the assumption of uniform hematocrit [44, 40]. For a given concentration of hemoglobin in blood vessels, it is theoretically possible to convert a UPD picture into a prior optical absorption coefficient due to blood, including in resolved vessels with blood volume of 100% and unresolved vessels with intermediate blood volume down to unresolved capillary beds. The conversion is done through a multiplication by a factor that is a priori unknown. We demonstrate it can be optimally determined within our proposed inversion procedure. The absorption due to blood is considered as dominant over other chromophores around 800 nm. In the general case, we then consider UPD gives the effective absorption coefficient locally within vascularized area. Furthermore, while PA imaging and UPD can both image blood, it is essential to account for resolution differences between both techniques. However, resolution degradation will lead to a blurring of the vasculature seen by UPD, which can be problematic if the blurring overlaps with some heterogeneities. A strategy is further defined to use the prior in the case of an optical heterogeneity surrounding or covering some blood vessels.
We based our inverse algorithm on a previously-reported Monte Carlo radiance inverse problem [29]. Thanks to the prior brought by UPD, we show we can simultaneously recover the absorption and scattering coefficients in regions surrounded or covered by blood capillaries with or without prior background information. The media in this paper are composed of blood capillaries (or vessels) among a background medium that includes optical absorption and scattering heterogeneities that are considered as markers of pathologies. For the sake of generality, we assume these pathologies result in local changes of optical absorption, scattering, or both. We consider that these pathologies are not visible from the vascular imaging perspective which implies that they are solely brought to light from their optical properties.

In this issue, we present the inversion model in the general case, where $\mu_s$ and $\mu_a$ are retrieved from a complete PA measurement. We then present the principles of the UPD prior information used in this study and show simulated examples. We demonstrate the UPD prior mitigates the non-uniqueness encountered in single-wavelength and illumination PA imaging. We show the inversion does not need prior information on background, inclusions and hemoglobin concentration and is robust to the presence of noise. Finally, we address the resolution aspects with various decimation of a high resolution medium.

## 2. Theory

### 2.1. Monte Carlo inverse radiance algorithm

The reader can refer to the studies that developed this specific theoretical framework for PA imaging inversion using a Monte Carlo inverse radiance scheme [25, 29, 30]. We briefly expose here the main equations and formulas without repeating their demonstrations. We define the ideal PA reconstruction as $H(\mathbf{r}) = \mu_a(\mathbf{r}).\Phi(\mathbf{r})$ where $\mu_a(\mathbf{r})$ is the optical absorption in the medium and $\Phi(\mathbf{r})$ is the fluence. We distinguish the "measured" reconstruction $H_{meas}$ and the reconstruction by the model at given iterations $H(\mu_a, \mu_s)$. The radiance $\phi(\mathbf{r}, \hat{s})$ is solution to the RTE, also referred to as forward problem:

$$ (\hat{s} \cdot \nabla + \mu_a(\mathbf{r}) + \mu_s(\mathbf{r}))\phi(\mathbf{r}, \hat{s}) - \mu_s(\mathbf{r}) \int_{S^{N-1}} P_{\theta}(\hat{s}, \hat{s}')\phi(\mathbf{r}, \hat{s}')d\hat{s}' = q(\mathbf{r}, \hat{s}) \quad (1) $$

where $\mathbf{r}$ is the position, $\hat{s}'$ and $\hat{s}$ are the original and scattered propagation directions, $P_{\theta}(\hat{s}, \hat{s}')$ is the scattering phase function, $q(\mathbf{r}, \hat{s})$ is a source term and $S^{N-1}$ indicates integration over angle in $N-1$ dimensions ($N=2$ in 2D). The forward model is run based on an initial guess of the model parameters (priors). The inversion is then performed using a gradient descent algorithm. A classical way is to derive functional gradients from the definition of the adjoint equation:

$$ (-\hat{s} \cdot \nabla + \mu_a(\mathbf{r}) + \mu_s(\mathbf{r}))\phi^*(\mathbf{r}, \hat{s}) - \mu_s(\mathbf{r}) \int_{S^{N-1}} P_{\theta}(\hat{s}, \hat{s}')\phi(\mathbf{r}, \hat{s}')^*d\hat{s}' = Q_{adj} \quad (2) $$
where $\phi^*$ is the solution to the adjoint equation, $Q_{adj} = \mu_a \cdot (H - H_{meas})$ is the source term of the adjoint problem. The error function is defined as: $\epsilon(\mu_a(\mathbf{r}), \mu_s(\mathbf{r})) = \frac{1}{2} \int (H(\mathbf{r}, \mu_a(\mathbf{r}), \mu_s(\mathbf{r})) - H_{meas}(\mathbf{r}))^2 d\mathbf{r}$.

To update the model at each iteration, one has to run the forward and the adjoint problem, allowing to calculate functional gradients as :

$$ \frac{d\epsilon}{d\mu_a} = -\Phi(\mu_a, \mu_s)(H_{meas} - H(\mu_a, \mu_s)) + \int_{S^{N-1}} \phi(\hat{s}, \mu_a, \mu_s)\phi^*(\hat{s}, \mu_a, \mu_s)d\hat{s} \\ \frac{d\epsilon}{d\mu_s} = \int_{S^{N-1}} \phi(\hat{s}, \mu_a, \mu_s)\phi^*(\hat{s}, \mu_a, \mu_s)d\hat{s} - \int_{S^{N-1}} \int_{S^1} \phi(\hat{s}, \mu_a, \mu_s)P_{\theta}(\hat{s}, \hat{s}')\phi(\hat{s}', \mu_a, \mu_s)^*d\hat{s}' \quad (3) $$

Monte Carlo light simulations usually output the angularly integrated fluence. In order, to compute the radiance scalar products, a decomposition in an angular harmonics basis is performed. In 2D, we use the circular harmonics basis (also called classical Fourier basis). Note that in 3D, one would use the spherical harmonics basis. During the simulation, each photon will impact local Fourier coefficients on their path, depending on their propagation direction resulting in $a_n(\mathbf{r})$ and $b_n(\mathbf{r})$ maps for the forward simulation and $a_n^*(\mathbf{r})$ and $b_n^*(\mathbf{r})$ [30].

During the Monte-Carlo transport simulation, each pixel/voxel accumulates the weights $\delta W$ deposited by the photons along their paths. The Fourier coefficients are updated following, for the 2D case, $a_0 = \frac{\delta W}{\pi}$, $a_n = \frac{\delta W}{\pi} \cos(n\theta')$ and $b_n = \frac{\delta W}{\pi} \sin(n\theta')$, where $\theta'$ is the incoming angle of the photon.

As demonstrated previously [30], using the Fourier decomposition, the gradients (3) become in the 2D case:

$$ \frac{d\epsilon}{d\mu_a} = -\Phi(\mu_a, \mu_s)(H_{meas} - H(\mu_a, \mu_s)) + \frac{\pi}{2}a_0a_0^* + \pi \sum_n (a_na_n^* + b_na_n^*) \\ \frac{d\epsilon}{d\mu_s} = \pi \sum_{n=1} [a_na_n^* + b_na_n^*](1 - g^n) \quad (4) $$

where $g$ is the anisotropy factor of the medium. The computation of this gradient at each iteration allows updating $\mu_a(\mathbf{r})$ and $\mu_s(\mathbf{r})$ towards the optimal solution.

## 2.2. Quantitative Power Doppler

After removing the tissue signals with a filtering step, the power Doppler image $\sigma_{US}^2(\mathbf{r})$ is calculated as the temporal variance in each pixel. In a previous publication, we have introduced a quantitative Power Doppler formula involving the imaging part [45]:

$$ \sigma_{US}^2(\mathbf{r}) = P_i^2(\mathbf{r}) \times \left(\frac{\Delta Z(\mathbf{r})}{Z_0}\right)^2 \times \eta(\mathbf{r})W[\eta(\mathbf{r})]V_{RBC}(\mathbf{r}) \times \int f(\mathbf{r}')|h_{US}|^2(\mathbf{r}', \mathbf{r})d\mathbf{r}' \quad (5) $$

where $P_i$ is the incident acoustic pressure, $\Delta Z(\mathbf{r})/Z_0$ is the relative mean impedance fluctuation induced by a single red blood cell, $\eta$ is the volume fraction of blood
(hematocrit), $W(\eta)$ is the packing factor, $h_{US}$ is the US imaging point spread function (PSF) and $f(\mathbf{r})$ defines the blood vessels. This formula is composed of a product between excitation amplitude $P_i$, medium parameter $\Delta Z(\mathbf{r})/Z_0$, statistical medium blood parameters $\eta W(\eta)$ and object-PSF integral (on the right side of the equation). Object-PSF integral involves the squared envelope of the US PSF which, thanks to its positivity, do not allow destructive interferences and enable a faithful representation of the blood vessels at the given resolution. The square factor enables a reduction of the effects of potential side lobes. The object-PSF integral is proportional to the blood volume $BV(\mathbf{r})$, partial volume of vessels in the pixel, or equivalently fractional vessel density.

It is important to note that modern Power Doppler, involving spatio-temporal SVD-filtering [46], efficiently reconstructs all the blood vessels at the given resolution and is insensitive to the flow speed and direction, unlike other modes such as Color Doppler. In the case the PSF is non-stationary, the blood volume image is derived from a normalization of the object PSF integral by the local squared L2 norm of the PSF: $||h(\mathbf{r})||_2^2 = \int |h_{US}|^2(\mathbf{r}', \mathbf{r})d\mathbf{r}'$. The normalized object-PSF integral is equal to 1 in vessel larger than the PSF all over the imaging domain. Assuming stationary blood properties (uniform hematocrit and red blood cell impedance and volume) and a correction of the image by a calibrated incident pressure $P_i$, one can assume the whole normalized image $\sigma_{US,n}^2$ is proportional to the blood volume:

$$ \sigma_{US,n}^2(\mathbf{r}) = (P_i^2(\mathbf{r}).||h(\mathbf{r})||_2^2)^{-1} \sigma_{US}^2(\mathbf{r}) \propto BV(\mathbf{r}) \quad (6) $$

This holds true in resolved vessels that are delineable on the images but also in unresolved small vessels. Using a wavelength range where hemoglobin absorption dominates over the other chromophores (at 800 nm), Power Doppler reflects the pixel effective absorption coefficient. In a previous study [45], we demonstrated experimentally that the ratio between Photoacoustic fluctuation imaging and UPD is proportional to the fluence in blood vessels, proving that Power Doppler detains the blood absorption quantity.

Blood vessels regions are defined, as a domain $M_{BV}$ satisfying $\mathbf{r} \in M_{BV}$ if $\tilde{\sigma}_{US}^2(\mathbf{r}) > \epsilon$, where $\epsilon$ is a threshold set on the UPD image amplitude. Thus, we consider that the absorption coefficient in blood vessels can be derived from the Power Doppler image through:

$$ \text{if } \mathbf{r} \in M_{BV}, \mu_{a,blood}(\mathbf{r}) = K_1.K_2.\sigma_{US,n}^2(\mathbf{r}) = \xi.\sigma_{US,n}^2(\mathbf{r}) \quad (7) $$

where $K_1 = 1/\sigma_{US}^2(BV(\mathbf{r}) = 1)$ is constant dependent on the imaging device, which can be eventually calibrated. It is the inverse of the value taken by the US Doppler picture in a large vessel that has a maximum $BV = 1$ value. Besides, $K_2$ depends on the hemoglobin and blood characteristics of the subject, that can vary among individuals but can be considered uniform within the same individual. Overall, $\xi = K_1.K_2$ is the resulting factor. To summarize, the Power Doppler image only gives a relative $\mu_a$ value
where blood vessels are present. $\xi$ is a priori unknown, we further show $\xi$ can be retrieved from a one-dimensional optimization procedure.

When a fluence model can be applied to the measurement, photoacoustic imaging can also probe blood volume through the total hemoglobin concentration in the pixel [47]. This does not hold true when spectral coloring affects the measurements or if the medium optical properties are unknown unless inverse problems are used.

In [45], quantitative fluence imaging is achieved experimentally after compensating the resolution differences between the UPD and PA fluctuation images. This was done in the most general case of non-stationary PSF in 3D, using a PSF filter defined to degrade the resolution of UPD to the one of PA fluctuation imaging. For a N-points spatial grid, the PSF filter is a NxN matrix where each column represents a kernel that applies to each position. In the present framework, this PSF filter can be obtained from model-based matrices $H_{PA}$ and $|H_{US}|^2$ containing all the reconstructed PSF for each point of the grid. Consequently, this formalism can be applied to degrade $\sigma_{US}^2$ to the PA resolution. However, in this study, we assume no resolution mismatch between Power Doppler and PA imaging. Nonetheless, such a PSF filter could be useful for an experimental validation of our findings.

# 3. Methods

## 3.1. Implementation

To compute the light radiance inside a domain of size $N_z \times N_x$ defined by optical parameters $\mu_a(\mathbf{r})$, $\mu_s(\mathbf{r})$ and $g(\mathbf{r})$, we use a Monte-Carlo implementation of the radiative transfer equation (RTE) called MCXLAB [48], with GPU acceleration. This open-source software offers light fluence computing in heterogeneous media. We modified the source code in order to project the radiance onto the harmonic basis following the implementation proposed by Hochuli et al. [29]. The modified version of MCXLAB was compiled on a UNIX-based operating software (OS) and then on a Microsoft Windows OS to achieve this study. MCXLAB is piloted from MATLAB 2021a (Mathworks, MA, USA) as is the inversion algorithm based on outputs from MCXLAB.

## 3.2. Inverse problem implementation

The forward part simulates a real, physical event where the photons propagate from the source in an absorbing and scattering medium. They lose part of their weight along their path until they reach a threshold (typically $10^{-6} W_0$ where $W_0$ is their initial weight) and die. The adjoint model, however, relates to a virtual situation where the photons propagate backward, due to the minus sign before the convection term (2). With a change of variable $\hat{s} \rightarrow -\hat{s}$, we see that the field $\phi(\mathbf{r}, -\hat{s})^*$ is solution to the original RTE equation. Thus, the same model is run for adjoint to get $\phi(\mathbf{r}, -\hat{s})^*$ and the angular variable is taken as its opposite (reversed directions) to obtain $\phi(\mathbf{r}, \hat{s})^*$. Before computing the scalar products, the photon directions are reversed, ie $\theta' => \theta' + \pi$
where $\theta'$ represents the incoming angles of arrival of the photon. Equivalently, the odd Fourier coefficients ($a_1^*, b_1^*, a_3^*, b_3^*$, etc) change sign. To maximize the reach of the photons emitted in the adjoint simulation, we consider a loss less adjoint propagation by conserving the weight of the photons along their paths. The photons die when they leave the domain.

A Total Variation (TV) regularization [49] is used to dampen the effects of noise inherent to MC simulations and potentially to PA measurements. The TV of $\mu_a$ or $\mu_s$ is defined as

$$T(\mu) = \int_{\Omega} \|\nabla \mu(\mathbf{r})\|_2 d\mathbf{r} \qquad (8)$$

Note that the Euclidean norm is not taken to the square. A prox-function smoothing version using the Huber function instead of the L2 norm is used in order to make it a derivable function and compute the gradients $\nabla T_{\tau(\mu_s)}$ and $\nabla T_{\tau(\mu_a)}$ where $\tau$ is the proximal smoothing parameter [49]. Compared to the Tikhonov regularization, the TV presents the advantage of preserving the edges and the high values. The cost function and gradients are changed to:

$$\epsilon(\mu_a(\mathbf{r}), \mu_s(\mathbf{r})) = \frac{1}{2} \int (H(\mathbf{r}) - H_{meas}(\mathbf{r}))^2 d\mathbf{r} + \alpha_a T(\mu_a(\mathbf{r})) + \alpha_s T(\mu_s(\mathbf{r})) \\ \left[ \frac{d\epsilon}{d\mu_a} \right]_{TV}(\mathbf{r}) = \frac{d\epsilon}{d\mu_a}(\mathbf{r}) + \alpha_a \nabla T_{\tau}(\mu_a(\mathbf{r})) \qquad (9) \\ \left[ \frac{d\epsilon}{d\mu_s} \right]_{TV}(\mathbf{r}) = \frac{d\epsilon}{d\mu_s}(\mathbf{r}) + \alpha_s \nabla T_{\tau}(\mu_s(\mathbf{r}))$$

The proximal smoothing $\tau$ was fixed to $5.10^{-4}$. The TV regularization parameters $\alpha$ were adjusted depending on the situation especially if noise is added ($\alpha_s = 1.10^{-5}$ or $1.10^{-4}$ and $\alpha_a = 3.10^{-5}$ or $3.10^{-4}$, resp. without or with noise).

An Adam scheme was used to optimize the convergence rate [50]. We consider a 8 × 8mm 2D medium illuminated by a single homogeneous source on the top of the sample with a 0.1mm step size. $5.10^6$ photons are launched at each iteration on a standard desktop PC with a NVIDIA GeForce RTX 3090 GPU, each iteration taking about 2.0 seconds.

### 3.3. Tissue model and image generation

Retinal capillary structures are used to represent blood vessels, taken from the DRIVE dataset (https://drive.grand-challenge.org/) [51]. We consider these images give a binary definition of the vessels at high resolution (see Fig.1a). Equation (5) states the Power Doppler is obtained by a convolution between the object and the squared envelope of the US PSF $|h_{US}(\mathbf{r})|^2$. $|h_{US}(\mathbf{r})|^2$ is a smooth function typically of the form of a 2D gaussian function. This convolution performs a low-pass filter on the object and the result can be downsampled without loss of information. Effectively, Power Doppler represents the blood volume, which is defined as a relative amount of blood in each given
pixel. Our Power Doppler image shown in Fig.1b is obtained by this decimation process,
resulting in the blood volume quantity (between 0 and 1). To simplify the test object,
resulting values lower than 0.18 are set to zero. Local blood absorption is also governed
by the amount of blood in each pixel. Therefore, to obtain the input $\mu_a$ map, we convert
the blood volume picture to blood absorption using the $\xi$ coefficient (whole blood:
$\mu_a = 4 \text{ cm}^{-1}$ at 800 nm) and we add the inclusions ($\mu_a = 0.8 \text{ cm}^{-1}$ and $\mu_a = 1 \text{ cm}^{-1}$)
and background values ($0.25 \text{ cm}^{-1}$) resulting in the Fig.1c. The anisotropy factor
and refraction index are supposed homogeneous and taken as $g = 0.9$ and $n = 1.37$.
The scattering coefficient of blood is assumed equal to the one of background, as a
simplification but it could be set to a more accurate representative value if experiments
require this in the future. The $\mu_s$ map is defined as shown in Fig.1d, with one of the
scattering inclusions co-aligned with one of the absorbing inclusions ($\mu_s = 30 \text{ cm}^{-1}$,
$\mu'_s = 3 \text{ cm}^{-1}$) and another on top left ($\mu_s = 40 \text{ cm}^{-1}$, $\mu'_s = 3 \text{ cm}^{-1}$) and a background
value of $\mu_s = 10 \text{ cm}^{-1}$ ($\mu'_s = 1 \text{ cm}^{-1}$). The addition of these scattering inclusions allows
demonstrating that the method is efficient with an arbitrary heterogeneous medium in
terms of scattering properties. All the optical coefficient values describe a realistic tissue
as previously reported in [52, 53].

The incident illumination is defined as a uniform plane wave coming from the top.
Following a forward Monte-Carlo simulation giving the fluence $\phi(\mathbf{r})$, the PA image is
obtained by the product $\mu_a(\mathbf{r}).\phi(\mathbf{r})$ (see Fig.1e).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Model and measurements. a) Vascular network high resolution binary map. -b) Power Doppler image. c-d) True $\mu_a$ and $\mu_s$ maps. e) PA imaging $H_{meas}$. f) Binary thresholded vascular map. The white color indicates the region where $\mu_a$ is updated through the iterations, $\mu_a$ is locked in the black vessel region. Units are $\text{mm}^{-1}$ for absorption and scattering</figcaption></figure>
### 3.4. Priors

To start the inversion, the numerical model needs to be initialized with $\mu_a$ and $\mu_s$ maps, called priors $\mu_{a,0}$ and $\mu_{s,0}$. In the presented examples, we always consider a uniform prior $\mu_{s,0}$ that can be equal to the true background value or not. The same principle is applied to $\mu_{a,0}$ except we use the $\mu_{a,blood}$ values in blood vessels (for $\mathbf{r} \in M_{BV}$). Note that the $\mu_{a,blood}$ and associated conversion factor $\xi$ is initially considered to be the exact one but we further show we can recover $\xi$ if unknown.

### 3.5. Inverse algorithm

The MC forward model is run using the defined medium properties to obtain the "measured" PA image $H_{meas}$ which would be measured in an experiment. In parallel, a numerical model is defined and the previously defined priors allow running a first forward model to obtain the model PA image $H_0$ (see Fig.2).

At each iteration, the model PA image $H_n$ is computed from $\mu_{a,n}$ and $\mu_{s,n}$ and compared to the measured value $H_{meas}$, leading to the adjoint source $Q_{adj} = \mu_{a,n} (H_{meas} - H_n)$. The adjoint propagation is computed from that source, resulting in radiance coefficients which, combined with the forward coefficients, give the gradients expressed in equations (4) (Fig.1). Since the source $Q_{adj}$ is signed (positive or negative), 2 simulations are run separately with respectively the positive and the negative part, then the results are subtracted. 4 Fourier harmonics were found sufficient in our model to construct correct gradients.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Summary of model inversion using a prior coming from Ultrasound Power Doppler.</figcaption></figure>
### 3.6. Optimization of the Power Doppler - $\mu_a$ coefficient

Equation (7) expresses the absorption coefficient in blood vessels. As previously stated, the factor $\xi$ presented in the equation depends on the hemoglobin and blood characteristics of the subject, which can vary among individuals but can be considered uniform within the same individual. Quantifying and fixing this factor is crucial for the convergence of the proposed inverse problem, as we will demonstrate in the results section. The absorption coefficient in blood vessels can potentially be quantified by taking a blood sample from the specific patient and measuring it through a spectrophotometer [54]. However, in an in-vivo application, this approach would make ineffective the non-invasive benefits of photoacoustic and ultrasound imaging. Hence, we propose a solution that relies on a one-dimensional optimization where the only unknown is $\xi$. Specifically, we implemented a simple optimization method, which is performed only once before the start of the entire inverse radiance algorithm. Starting from an initial guess, $\xi$ is changed using a fixed step size of $\pm 0.5 \cdot 10^{-3} \text{mm}^{-1}/\text{a.u.}$, in a direction determined by the initial sign of $d\epsilon_{\mu_a}/d\mu_a$. The error function defined in section 2.1 decreases and we let it re-increase for 10 iterations: its minimum is recovered. The routine requires around 60 iterations to find the optimal $\xi$. The validation of this method is discussed in section 4.3 and applied to all the results presented subsequently.

### 3.7. Quantitative aspects

In some of the presented results, some additive noise is added to the input PA image. This noise is generated using a normal distribution of given standard deviation $\sigma_{image}$. We define Signal-to-Noise Ratio (SNR) as the ratio between the maximum of the image and background noise standard deviation. To compare different situations, a quantitative analysis of results is presented in terms of normalized root mean square error:

$$ NRMSE = \left\langle \left( \frac{\mu(\mathbf{r}) - \mu_{\text{true}}(\mathbf{r})}{\mu_{\text{true}}(\mathbf{r})} \right)^2 \right\rangle \quad (10) $$

where $\mu$ can be either $\mu_s$ or $\mu_a$ and $\langle\rangle$ denotes spatial average over the whole domain.

### 3.8. Resolution aspects

We show the effect of a resolution degradation in results in Fig.6. To obtain the UPD picture, the high-resolution binary vessel map is filtered by a 2D Gaussian before getting decimated. In parallel, the true $\mu_a$ map is first defined at high resolution with binary vessels and inclusions. Fluence distribution is simulated in the medium at high resolution and the PA image $H_{0,HR}$ is computed. Resolution is then degraded the same way as UPD to obtain the input $H_0$ of the inversion. The width of the 2D Gaussian is 20, 40, and 60 pixels corresponding to 0.4, 0.8, and 1.2 mm. The decimation factor remained at the same value of 5 leading to the same image size for the 3 cases.

Resolution degradation blurs the contour of the vascular network which eventually creates pixels containing mixed inclusion and blood. As our algorithm is based on
keeping the prior blood absorption values in blood vessels, the inclusion can not be retrieved this way. To optimize a region containing a supposed inclusion, the prior blood absorption is set to the background $\mu_a$ value in a disk enclosing the inclusion. This means the prior blood absorption is only used around the inclusion. We will show our approach still has a beneficial effect on stabilizing a unique solution in the inversion. The way $\mu_s(\mathbf{r})$ is treated is not affected by this resolution degradation.

# 4. Results

## 4.1. Effect of blood prior

We first consider the case where background $\mu_a$ and $\mu_s$ are known. Fig.3 shows the inversion results with (c-e) and without knowledge (c'-e') of the vascular structure. Without any prior, the reconstruction of both $\mu_a$ and $\mu_s$ is impossible due to the ill-posed problem. With the vascular information, $\mu_a$ is known and fixed in the vessels. This constraint is sufficient to overcome the non-uniquity problem and allows the model to converge to an accurate solution (middle column). We remind that when the prior is considered, $\mu_a$ is not updated in the blood vessels throughout the inversion. However, when no prior is considered, $\mu_a$ is indeed updated everywhere along with $\mu_s$. The NRMSE for $\mu_a$ is 0.76% with prior and 12.9% without prior, and for $\mu_s$ the NRMSE is 8.71% with prior and 24.8% without prior. Interestingly, the cost function minimal value ends up higher with the prior $(1.2 \times 10^{-7})$ than without prior $(0.9 \times 10^{-7})$ which depicts how non-uniqueness can lead to false results. However, using the prior leads to a 16-fold decrease in the initial cost function value.

## 4.2. Effect of inaccurate background priors

We consider an error of +20% for both $\mu_a$ and $\mu_s$ priors in the background. We consider the blood prior as exact. The recovered $\mu_a$ in Fig.4a.2 is very close to the true $\mu_a$ in Fig.3a with an NRMSE of 4.7%. The effect is more visible for $\mu_s$ with a recovered background slightly perturbated (see Fig.4b.2) than in Fig.3e. The NRMSE value is increased to 23%. Most importantly, convergence is conserved when inaccurate background priors are considered.

## 4.3. Effect of inaccurate blood absorption value

In section 3.6, we presented a secondary algorithm that can be run to initially optimize the $\mu_a$ blood value through the $\xi$ coefficient 1D optimization. Here, while keeping the situation of Fig.4 with wrong background priors, we consider an error in the prior $\xi$ coefficient of 10%. Keeping this false blood absorption still leads to a solution, shown in Fig.5a1-b1. The recovered $\mu_a$ is not visually affected (NRMSE=9.5%) but the $\mu_s$ reconstruction is degraded (NRMSE = 31.9%).

The optimization routine of $\xi$ is performed to lead to improved results (see Fig.5a2-b2, with NRMSE=5.0% for $\mu_a$ and 26.5% for $\mu_s$). $\xi$ is only optimized at the first iteration
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Results of the inversion of $\mu_a$ and $\mu_s$ with a prior on $\mu_a$ in the vessels (c-e) and without (c'-e'). True values are displayed in a-b). Units are $mm^{-1}$</figcaption></figure>

of the inversion, based on the homogeneous wrong $\mu_s$ and $\mu_a$ background priors. The evolution of the cost function during this 1D-optimization is shown in Fig.5c. Through this 62 iterations process, the optimal $\xi$ has an absolute error of 1.95% compared to the true value. This result shows considering the blood prior map is effective even if the absolute blood absorption and background priors are erroneous.

### 4.4. Effect of resolution

The results presented in this last section account for resolution degradation by considering the effect of the PSF on both photoacoustic and ultrasound prior images. The three cases shown in columns 1, 2 and 3 of Figure 6 involve, as specified in section 3.8, a degradation of $H_0$ via a 2D Gaussian filtering, with a width of the Gaussian of 20, 40, and 60 pixels, respectively. To simulate a more realistic scenario, these examples include not only the PSF-object convolution but also all the conditions previously mentioned, as a 20% error in the background values of $\mu_a$ and $\mu_s$, and an unknown $\xi$. Additionally, some noise was added to the image before resolution degradation. As an initial step in the inversion, the $\xi$ parameter has been estimated using the optimization algorithm, yielding errors of 2.5%, 5.8%, and 6.6% for the three examples.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Inversion of $\mu_a$ and $\mu_s$ with wrong background priors. Units are $mm^{-1}$. a.1-b.1 Priors with 20% error in background values. a.2-b.2 Inversion results. True $\mu_a$ and $\mu_s$ maps are found in Fig.3a-b.</figcaption></figure>

Before discussing the findings, we computed the $SNR_{image}$ on the degraded $H_0$ images using the equation provided in section 3.7. The signal-to-noise ratios were evaluated for the first 15 pixels at the top and the last 15 pixels at the bottom of the image. Across the three degradation scenarios, the SNR at the top was 26.7, 23.9, and 19.9, while at the bottom, it was 9.2, 11.8, and 14.9. The decrease in fluence is responsible for a lower SNR at the deeper region.

From the inversion of $\mu_s$ presented in Fig.6f.1-f.3, we can notice that the recovered values are not significantly affected by the resolution degradation (resp. NRMSE = 33.9%, 43.9% and 42.9%). Considering the absorption, Fig.6e.1-e.3 shows that the retrieved $\mu_a$ closely matches the PSF-convolved $\mu_a$ in Fig.6c.1-c.3. The NRMSE values are equal to 76.5%, 44.1% and 31.0%.

It is worth noting that the inversion is performed using a PSF-convolved $H_0$, and the PSF-convolved $\mu_a$ is not provided to the algorithm. Nevertheless, as fluence is slowly varying, the PSF-convolved $\mu_a$ is a good target to compare with the inversed $\mu_a$.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Optimization of the UPD to blood absorption coefficient $\xi$ starting with a 10% error and with a 20% error in background properties. a1-b1 Result without optimization of $\xi$. a2-b2 Result with optimization of $\xi$. c Trend of the cost function as the optimization algorithm iterates. The marker indicates the optimal point, representing the value of $\xi$ resulting in the lowest cost.</figcaption></figure>

**5. Discussion**

In this issue, we propose a strategy to add valuable information for quantitative photoacoustic imaging. We show how Ultrasound Power Doppler can help to estimate a prior absorption map that reduces the number of unknowns and leads to accurate inversions of $\mu_a$ and $\mu_s$. Thanks to the prior $\mu_{a,blood}$ that covers the whole domain, the adjoint sources, which are strongly dependent on current $\mu_a$ values are more effective to minimize the cross-talk between $\mu_a$ and $\mu_s$. We show the algorithm is robust to unknown background as well as unknown blood properties. While degrading the imaging resolution, we present an effective strategy when there is no clear separation between optical heterogeneities and blood vessels, in the presence of noise.

The implementation of a radiance Monte Carlo inverse problem was done using MCXLAB. This was made possible by adding features to the source code of MCXLAB to retrieve harmonic coefficients describing the angular information of radiance from the simulation. This preliminary study presents only 2D results. The 3D nature of optical transport and PA imaging requires a 3D framework for working on real data. The exposed implementation and theory are ready to use for 3D data and will be part of a future study. Our algorithm uses an Adam optimizer which speeds up convergence but the computational cost could be reduced, for example using Machine learning to adjust the number of needed photons in each Monte Carlo simulation [55].
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Effect of resolution added to background errors, $\mu_{a,blood}$ error, and noise. a-b: True high-resolution $\mu_a$ and $\mu_s$. Columns 1, 2 and 3 refer to the degradation via a 2D Gaussian filtering, with a width of the Gaussian of 20, 40, and 60 pixels, respectively. Line c: degraded $\mu_a$, for comparison to the recovered $\mu_a$ in line e. Line d: $\mu_a$ priors including the disk set to the background value enclosing the inclusion. Line e: Recovered $\mu_s$.</figcaption></figure>

We applied the inverse model to a medium composed of absorption and scattering inclusions surrounded by some blood vessels. We assumed prior $\mu_a$ could be retrieved from a combined vascular imaging modality (based on ultrasound Power Doppler), allowing the model to start with more knowledge on the medium. The conversion factor $\xi$ can be optimized at different times of the inversion. Different strategies were tried such as optimization every 5 iterations but this turned out to be less efficient compared to a single optimization at the beginning, at least with our medium model.
The ξ optimization algorithm was chosen here very simple and there is room to improve
its efficiency, which was not our goal here.

Note that the adjoint equation is a modified RTE with a source term $Q_{adj}$ defined
as the residual error at each position **r**, multiplied by $\mu_a(\mathbf{r})$. Hence, $\mu_a(\mathbf{r})$ has a great
impact on this source term which implies that the inversion will be much faster if a
$\mu_a(\mathbf{r})$ prior covering most of the domain is used. In the presented situation where blood
is the strongest absorber, this means that the vessels are the main contributors to $Q_{adj}$.
We can see in the result of Fig.4b.2) that the algorithm is unable to update the borders
of the medium properly. This might be due to the low values taken by $Q_{adj}$ when $\mu_a$
values are small.

We show that the sensitivity to noise can be improved using a TV-regularization.
This type of regularization uses an empirically-tuned parameter. An automatic
algorithm could be useful for general applications. A strong assumption in this study is
that the PA reconstructions are "perfect" and so are UPD reconstructions. Behind this
assumption, lies the fact that PA imaging and UPD have the same PSF and thus, image
in the same manner which is not true in the general case. For example, ultrasound Power
Doppler lateral and axial resolution is different from the one obtained in PA imaging
[43] [45], since ultrasound compounding performs focus in emission, and also since the
round-trip nature of ultrasound improves the axial resolution two-fold. If we apply the
method proposed here on real data, it means that the prior $\mu_a$ from vascular imaging
would not have the same native resolution as $H_{meas}$. $Q_{adj}$ being the product of two
different quantities with different resolutions, it will have strong spatial fluctuations,
for example, on the border of a smooth object, which might prevent converging to
an acceptable solution. For application to real data, a resolution-equalization needs
to be performed and the easiest manner to proceed is to degrade the best resolution
towards the worst of the two techniques as performed in [45]. Applying model-based
reconstructions for both modalities might still result in different resolutions [56].

Blood vessels can have diameters ranging from centimeters down to 10µm. This
multiscale nature has been approached in this study with the resolution degradation. In
UPD, the evaluated quantity on the image is the blood volume [40], the same quantity
as MRI or CT can evaluate [57]. Blood volume images are proportional to the blood
volume fraction within each pixel. A maximum value is reached when the pixel is
saturated with blood. Particularly, when observing unresolved capillary-bed structures,
blood volume is still quantitatively determining the volume fraction of blood. Thus,
correct mean effective $\mu_a$ pixel values can still be imposed by this adequate vascular
imaging techniques in the case of unresolved vessels.

However, when resolution blurs the contour of the vasculature prior, it is necessary
to ignore the prior blood absorption in a region that needs to be updated in $\mu_a$. This
strategy supposes that there is a way to identify the presence of a $\mu_a$ heterogeneity
not due to flowing blood. This situation can also arise in the presence of additional
chromophores (lipid, melanin, contrast agent accumulation) inside or very close to some
vessels. The delineation of this region could be done manually or automatically based
on cost function values after a first run of the inversion or based on an Ultrasound
anatomical image showing a suspicious structure. Herein, we considered vessels in a
given range of scales. When applying the method *in vivo*, we could face the case where
there is a higher variety of vessel scales, with a capillary bed composing the background.
In this case, two possible strategies can be considered. First, if UPD sensitivity allows
probing capillary bed structures as an effective vessel density, also responsible for the
background optical absorption, this information could lead to a prior optical absorption
in the background. Hence, the PA contrast will enable recovering scattering properties
and optical absorption contributions other than blood, provided that optical absorption
coefficients are left free to some extent during the updates, using regularizations, for
example. The other strategy involves considering the UPD prior only in larger vessels,
simply by using a threshold on the UPD images. In this case, the same method as
proposed here, which locks the prior $\mu_a$ only in larger vessels, will let the inverse problem
recover background optical absorption and heterogeneities. The practical use of either
strategy depends on the blood vessel structure as well as the PA reconstruction ability
to recover given spatial-frequency contents. These considerations will be the focus of
further work.

Visibility limited issues were not studied here. The limitations due to the detectors
were neglected. It could be interesting to study the effect of the vascular imaging prior
(which is not affected by visibility artefacts) in the inversion of PA visibility-limited data.
This implies considering a model of detectors and acoustic reconstruction in addition.

Another assumption of our study is that optical absorption in pixels containing
blood vessels is only due to hemoglobin. This is not the general case, since
other chromophores such as lipids, water and melanin can participate in absorption
and spectral unmixing needs to be performed before quantitative blood analysis.
However, our approach uses the isobestic point at 800 nm where the effects of other
chromophores are minimized [58]. Additionally, oxygenation modulations do not affect
blood absorption coefficients at this wavelength. In the advent of a highly-absorbing
heterogeneity mixed with blood, the strategy presented in section (3.8) is effective to
estimate the total $\mu_a$ in the region. Eventually, from the difference $\mu_a(\mathbf{r}) - \mu_{a,blood}(\mathbf{r})$,
one could estimate the absorption due to the materials other than blood.

Over such a small propagation distance, we consider acoustic attenuation can
be neglected. However, depending on the frequency band of the signals, acoustic
attenuation effects could be significant and would have to be corrected within the
employed reconstruction.

Our tissue model was chosen with a finer resolution in $\mu_a$ than in $\mu_s$. When testing
the inversion on $\mu_a$ and $\mu_s$ maps using smaller inclusions (results not shown here),
the resolution of the result was significantly higher for $\mu_a$, while such small inclusions
were not recovered for $\mu_s$. This happens because the PA signal is significantly more
sensitive to $\mu_a$ than to $\mu_s$. Additionally, the gradient analytical formulae contribute to
this difference, since $d\epsilon/d\mu_a$ contains a local term proportional to the measurement itself
$\Delta H$ whereas $d\epsilon/d\mu_s$ involves forward and adjoint radiance which plays a secondary role.
## 6. Conclusion

As part of an optical inverse problem, we show the combination of two modalities that measure blood vessels in different manners, Ultrasound Power Doppler and PA imaging, can bring the additional information needed for convergence. The fact that the support of information constituted by blood vessels is omnipresent in vivo and holds strong signals is a great opportunity. This methodology brings a promising strategy for quantitative imaging in photoacoustics.

## 7. Acknowledgments

This work was supported by the Agence National de la Recherche (ANR) under the FULBOX project (ANR-21-CE19-0003-01). The authors wish to thank Qianqian Fang for technical help with MCXLAB.

## 8. Data availability

Data will be made available upon reasonable request.

## References

[1] Beard P 2011 *Interface focus* **1** 602–631

[2] Taruttis A and Ntziachristos V 2015 *Nature photonics* **9** 219–227

[3] Wang L V and Hu S 2012 *science* **335** 1458–1462

[4] Deán-Ben X L, Ford S J and Razansky D 2015 *Scientific reports* **5** 10133

[5] Jathoul A P, Laufer J, Ogunlade O, Treeby B, Cox B, Zhang E, Johnson P, Pizzey A R, Philip B, Marafioti T et al. 2015 *Nature Photonics* **9** 239–246

[6] Cox B T, Arridge S R and Beard P C 2009 *Journal of the Optical Society of America A* **26** 443 ISSN 1084-7529, 1520-8532

[7] Laufer J, Delpy D, Elwell C and Beard P 2006 *Physics in Medicine & Biology* **52** 141

[8] Tzoumas S, Nunes A, Olefir I, Stangl S, Symvoulidis P, Glasl S, Bayer C, Multhoff G and Ntziachristos V 2016 *Nature communications* **7** 12121

[9] Lutzweiler C, Meier R, Rummeny E, Ntziachristos V and Razansky D 2014 *Optics letters* **39** 4061–4064

[10] Li P C, Wei C W, Liao C K, Chen C D, Pao K C, Wang C R C, Wu Y N and Shieh D B 2007 *IEEE transactions on ultrasonics, ferroelectrics, and frequency control* **54** 1642–1647

[11] Hochuli R, An L, Beard P C and Cox B T 2019 *Journal of biomedical optics* **24** 121914–121914

[12] Jeng G S, Li M L, Kim M, Yoon S J, Pitre Jr J J, Li D S, Pelivanov I and O'Donnell M 2021 *Nature communications* **12** 716

[13] Brochu F M, Brunker J, Joseph J, Tomaszewski M R, Morscher S and Bohndiek S E 2016 *IEEE transactions on medical imaging* **36** 322–331

[14] Ranasinghesagara J C and Zemp R J 2010 *Journal of biomedical optics* **15** 046016–046016

[15] Bu S, Liu Z, Shiina T, Kondo K, Yamakawa M, Fukutani K, Someda Y and Asao Y 2012 *IEEE transactions on biomedical engineering* **59** 1354–1363

[16] Cox B, Laufer J G, Arridge S R and Beard P C 2012 *Journal of biomedical optics* **17** 061202–061202

[17] Kirchner T and Frenz M 2021 *Journal of biomedical optics* **26** 085001–085001

[18] Bench C, Hauptmann A and Cox B 2020 *Journal of Biomedical Optics* **25** 085003–085003
[19] Shi M, Vercauteren T and Xia W 2024 Learning to compensate spectral coloring in a led-based photoacoustic/ultrasound imaging system *Photons Plus Ultrasound: Imaging and Sensing 2024* vol 12842 (SPIE) pp 69-72

[20] Gröhl J, Kirchner T, Adler T J, Hacker L, Holzwarth N, Hernández-Aguilera A, Herrera M A, Santos E, Bohndiek S E and Maier-Hein L 2021 *Scientific reports* **11** 6565

[21] Liang Z, Zhang S, Liang Z, Mo Z, Zhang X, Zhong Y, Chen W and Qi L 2024 *Photoacoustics* **37** 100601

[22] Poudel J, Lou Y and Anastasio M A 2019 *Physics in Medicine & Biology* **64** 14TR01

[23] Huang C, Wang K, Nie L, Wang L V and Anastasio M A 2013 *IEEE transactions on medical imaging* **32** 1097-1110

[24] Arridge S R 1999 *Inverse problems* **15** R41

[25] Saratoon T, Tarvainen T, Cox B T and Arridge S R 2013 *Inverse Problems* **29** 075006 ISSN 0266-5611, 1361-6420

[26] Leino A A, Lunttila T, Mozumder M, Pulkkinen A and Tarvainen T 2020 *IEEE Transactions on Medical Imaging* **39** 2985-2995 ISSN 0278-0062, 1558-254X

[27] Hänninen N, Pulkkinen A, Arridge S and Tarvainen T 2022 *Journal of Biomedical Optics* **27** ISSN 1083-3668

[28] Macdonald C M, Arridge S and Powell S 2020 *Journal of Biomedical Optics* **25** ISSN 1083-3668

[29] Hochuli R, Powell S, Arridge S and Cox B 2016 *Journal of Biomedical Optics* **21** 126004 ISSN 1083-3668

[30] Buchmann J, Kaplan B A, Powell S, Prohaska S and Laufer J 2019 *Journal of Biomedical Optics* **24** 1 ISSN 1083-3668

[31] Zheng S, Yingsa H, Meichen S and Qi M 2023 *Physics in Medicine & Biology* **68** 065009

[32] Zemp R J 2010 *Applied optics* **49** 3566-3572

[33] Shao P, Cox B and Zemp R J 2011 *Applied optics* **50** 3145-3154

[34] Niederhauser J J, Jaeger M, Lemor R, Weber P and Frenz M 2005 *IEEE transactions on medical imaging* **24** 436-440

[35] Kolkman R G, Brands P J, Steenbergen W and van Leeuwen T G 2008 *Journal of biomedical optics* **13** 050510-050510

[36] Merčep E, Deán-Ben X L and Razansky D 2017 *IEEE transactions on medical imaging* **36** 2129-2137

[37] Pattyn A, Mumm Z, Alijabbari N, Duric N, Anastasio M A and Mehrmohammadi M 2021 *Photoacoustics* **23** 100275

[38] Yang H, Jüstel D, Prakash J, Karlas A, Helfen A, Masthoff M, Wildgruber M and Ntziachristos V 2020 *Photoacoustics* **19** 100172

[39] Murphy K J and Rubin J M 1997 Power doppler: it's a good thing *Seminars in Ultrasound, CT and MRI* vol 18 (Elsevier) pp 13-21

[40] Macé E, Montaldo G, Cohen I, Baulac M, Fink M and Tanter M 2011 *Nature methods* **8** 662-664

[41] Deffieux T, Demené C and Tanter M 2021 *Neuroscience* **474** 110-121

[42] Robin J, Özbek A, Reiss M, Dean-Ben X L and Razansky D 2022 *IEEE Transactions on Medical Imaging* **41** 846-856 ISSN 1558-254X conference Name: IEEE Transactions on Medical Imaging

[43] Godefroy G, Arnal B and Bossy E 2023 *Scientific Reports* **13** 2961 ISSN 2045-2322 number: 1 Publisher: Nature Publishing Group

[44] Rubin J M, Adler R S, Fowlkes J B, Spratt S, Pallister J E, Chen J F and Carson P L 1995 *Radiology* **197** 183-190

[45] Falco I, Bossy E and Arnal B 2024 *Physics in Medicine Biology*

[46] Demené C, Deffieux T, Pernot M, Osmanski B F, Biran V, Gennisson J L, Sieu L A, Bergel A, Franqui S, Correas J M et al. 2015 *IEEE transactions on medical imaging* **34** 2271-2285

[47] Wang X, Xie X, Ku G, Wang L V and Stoica G 2006 *Journal of biomedical optics* **11** 024015-024015

[48] Yu L, Nina-Paravecino F, Kaeli D and Fang Q 2018 *Journal of biomedical optics* **23** 010504-010504

[49] Jensen T L, Jørgensen J H, Hansen P C and Jensen S H 2012 *BIT Numerical Mathematics* **52**
329–356

[50] Kingma D P and Ba J 2014 *arXiv preprint arXiv:1412.6980*

[51] Khanal A and Estrada R 2019 Dynamic deep networks for retinal vessel segmentation (Preprint 1903.07803)

[52] Sandell J L and Zhu T C 2011 *Journal of biophotonics* **4** 773–787

[53] Jacques S L 2013 *Physics in Medicine & Biology* **58** R37

[54] Thomas A, Rietberg M, Akkus M, van Soest G and Francis K J 2024 *Photoacoustics* 100673

[55] Macdonald C M, Arridge S and Powell S 2020 *Journal of Biomedical Optics* **25** 085002–085002

[56] Vilov S, Arnal B, Hojman E, Eldar Y C, Katz O and Bossy E 2020 *Scientific reports* **10** 4637

[57] Zaharchuk G 2007 *American Journal of Neuroradiology* **28** 1850–1858

[58] Cox B T, Arridge S R and Beard P C 2009 *JOSA A* **26** 443–455