Abstract
In quantitative photoacoustic tomography (qPAT), raw photoacoustic (PA) images offer only an indirect representation of the structural and physiological details of biological tissues. To address this issue, precise absorption coefficient (AC) extraction algorithms are essential for converting PA images into accurate AC maps, which involves mitigating the effects of non-uniform light fluence (LF). This study employs a dual-modality approach using photoacoustic tomography and ultrasound (PAUS), where ultrasonic segmentation identifies key structural boundaries, guiding the construction of Monte Carlo (MC) optical transport models. Leveraging sparse signal representation theory, we introduce a novel quantitative reconstruction algorithm that efficiently separates LF components from PA signals, refining AC imaging accuracy in complex tissues. Through numerical simulations and experiments with tissue-mimicking phantoms and in vivo mouse models using a PAUS system, we demonstrate our algorithm's capability. Results indicate significant improvements in feature visibility, boundary definition, and overall image quality, along with enhanced structural and functional information. This work aims to advance the detection and quantitative imaging of blood vessels and organs.
© 2025 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Photoacoustic (PA) imaging is a hybrid biomedical modality that offers highly contrasting and sensitive insights into the structure and function of biological tissues [1–3]. It relies on the photoacoustic effect, wherein photoacoustic signals generated by short-pulse laser irradiation reflect the tissue's light absorption characteristics. By utilizing optical absorption contrast and ultrasonic spatial resolution, PA imaging can visualize both endogenous and exogenous contrast agents [3–6]. A key objective of quantitative photoacoustic tomography (qPAT) is to eliminate the influence of light fluence (LF) on images, allowing for precise quantitative estimation of the absorption coefficient (AC) of chromophores in biological tissues [7–9].
The contrast in photoacoustic imaging arises from light absorption, which is closely tied to tissue composition. As light travels through tissues, its heterogeneous absorption by various components and reflections at different boundaries result in uneven energy deposition. The detected photoacoustic response is influenced by the local photon density, causing the resulting PA image to reflect both the spatial variations in absorption coefficient and the corresponding LF. However, the accurate representation of the absorption coefficient is also affected by other factors, such as limited bandwidth and limited view effects. These factors can complicate the interpretation of images, leading to inaccuracies in the estimation of the local AC [8]. This leads to a situation where the local AC is not accurately represented in the PA image, blurring the contrast between target tissues and biological chromophores. Consequently, while photoacoustic imaging leverages the photoacoustic effect, the original images only indirectly depict the biological tissue structure and physiological information. To address this, it is essential to eliminate the depth-dependent decay of LF and accurately reconstruct the AC distribution, which directly reveals the spatial distribution of absorbers [9–12]. This approach ultimately enables a precise and quantitative estimation of chromophore distribution in tissues [10–25]
Many studies have focused on improving the accuracy of tissue absorption coefficient imaging by quantitatively estimating the absorption distribution [26–35]. Cox et al. developed a simple iterative method for quantitative reconstruction, fitting the AC by measuring the initial pressure map. Their algorithm converges to a precise AC estimate [11,12]. Additionally, Cox et al. employed a gradient descent scheme to recover either absorption or scattering coefficients when the other is known, using a two-dimensional model with two chromophores. They used a diffusion approximation model and finite difference method to reconstruct chromophore concentrations [13]. Banerjee et al. proposed a non-iterative method to restore optical absorption coefficients by achieving a more uniform photon density with multisource irradiation, calculating the AC distribution using recovered LF, with boundary measurements obtained via the Monte Carlo method [14]. Yang et al. introduced a deep residual recurrent neural network to generate quantitative blood oxygen images, extracting LF information from optical absorption images through Monte Carlo simulations [25]. Bu et al. integrated fluence compensation into the reconstruction process, iteratively updating the AC distribution by minimizing residuals until convergence to the true absorption distribution [26]. Yuan et al. developed a quantitative reconstruction method based on diffusion equations and the regularized Newton method, incorporating prior structural information from PAT image segmentation [27,28]. Zhang et al. proposed a pixel-level AC map using a non-segmented iterative algorithm for alternating optimization, showing improvements in image quality and feature visibility compared to segmentation-based methods [29]. Liang et al. developed an automatic 3D segmentation method for small animal PAT images, using an optimal graph search algorithm to simulate the volumetric LF distribution across the body and generate corrected images [30]. Hochuli et al. introduced a model-based reconstruction framework that uses gradient-based minimization and a finite element model to recover both optical absorption and scattering coefficients, particularly effective in heterogeneous tissue environments [36]. Ranjbaran et al. developed an iteratively refined wave field reconstruction inversion method, improving the precision of AC estimates by accounting for complex tissue interactions. Their simulation study demonstrated reduced reconstruction errors and enhanced image clarity [37].
Recent advances in qPAT have introduced various approaches, including diffusion approximation, iterative reconstruction, and deep learning-based methods. While diffusion approximation offers computational efficiency, it often fails to accurately model light propagation in heterogeneous tissues due to its assumption of homogeneous optical properties [38]. Iterative reconstruction methods, on the other hand, provide higher accuracy but are computationally expensive and prone to convergence issues [39]. Deep learning-based techniques, such as convolutional neural networks (CNNs) and physics-informed neural networks, have shown promise in improving reconstruction speed and accuracy; however, they require large datasets for training and lack interpretability [40]. Comparison of different methods in qPAT as shown in Table 1.
The challenge of quantitative photoacoustic imaging is to eliminate the effects of LF [10,31–42]. In order to obtain more accurate tissue absorption images, correction or elimination of the fluence distribution is necessary. The main challenge of model-based approaches comes from the difficulty in modeling the propagation of light in organizations. Such modeling requires prior knowledge of the optical properties of complex geometries related to biological tissues. This knowledge often includes incomplete estimates of absorption and scattering coefficients of spatial variation in the medium, which can lead to inaccurate photon prediction models and correction results.
In this work, we propose a novel qPAT reconstruction algorithm that leverages sparse decomposition and ultrasound-guided Monte Carlo (MC) modeling to address these limitations. Unlike conventional methods, our approach explicitly accounts for depth-dependent scattering and absorption variations, enabling more accurate fluence correction. Furthermore, compared to deep learning-based methods, our framework offers interpretability, data efficiency, and generalizability, making it suitable for a wide range of clinical and preclinical applications.
This paper proposes a novel method for extracting the AC from photoacoustic images by leveraging the general behavior of LF rather than relying on traditional light propagation models. The method decomposes the photoacoustic image into two components using sparse signal representation: a low spatial frequency global component, caused by photon diffusion, and a high spatial frequency local component, which reflects changes in the AC. This approach enables more accurate quantification of the AC in complex tissues by effectively separating the influence of LF from the photoacoustic signal. In this study, we employed a dual-modality photoacoustic and ultrasound (PAUS) imaging system, where ultrasound segmentation helps preserve key structural boundaries, offering prior structural knowledge to guide the construction of the Monte Carlo optical transport model. Based on sparse signal representation theory, a new quantitative reconstruction algorithm was developed to effectively separate LF components from photoacoustic signals. By correcting for non-uniform fluence and absorption, the accuracy of AC imaging in complex tissues is enhanced. We constructed a dictionary for absorption and fluence distribution using the kernel singular value decomposition (K-SVD) algorithm [34,35,41], and employed orthogonal matching pursuit (OMP), a modified version of matching pursuit, to improve decomposition accuracy [35,41,43–45].
2. Materials and methods
2.1. Experimental setup
In order to demonstrate the proposed method, a PAUS dual-mode imaging system was developed for experimental study. The system consists of 5 main parts: optical illumination, motion control platform, PA and US signal acquisition, data processing and display module. In the optical section, a Q-switched Nd:YAG laser (Vibrant 355 II HE, Opotek, USA) is utilized, emitting 5 ns pulses within a tunable near-infrared range. The laser operates at a repetition rate of 10 Hz and a wavelength of 532 nm. A custom-made optical fiber bundle ensures uniform target illumination. The laser fluence on the sample surface is capped at 13 mJ/cm2, adhering to ANSI safety standards. For PA and US imaging, a linear array transducer with 128 elements captures the signals. An external laser trigger activates the receive-only mode on the US platform. The transducer has a 0.30 mm pitch and each element operates at a 7.5 MHz center frequency with fractional bandwidth of 75%. Optical fiber bundles on both sides are integrated into a custom housing, enabling the design of a handheld PA probe. The linear array transducer is scanned by a motion control platform. Figure 1 shows a schematic of the PAUS dual-modality imaging system integrated setup.
Fig. 1. Schematic diagram of the PAUS dual-modality imaging system. FB: fiber bundle; TH: transducer holder. The laser trigger synchronizes all events involved (laser firings, PA receiving). Tx: transmit; Rcv: Receive.
2.2. Reconstruction algorithm
In photoacoustic imaging, when the pulsed laser irradiates the tissue surface, the light energy is deposited in the irradiating area, and part of the energy is absorbed by the tissue, which generates photoacoustic effect to excite ultrasound. The resulting initial pressure of PA can be described as Eq. (1):
Here, $\Gamma $ is the thermo-elastic Grüneisen parameter, ${\mu _a}(r)$ and ${\mu _{s^{\prime}}}(r)$ are absorption and scattering coefficients, and $F$ is the LF at position r. Initially, we assume $p(r)$ is accurately reconstructed from acoustic data with negligible distortion. In soft tissue, $\Gamma $ is nearly constant [18].
The absorbed energy map is reconstructed by measuring the photoacoustic signal emitted after absorbing the laser energy in the tissue [10,11]. Thus, the absorbed energy $p^{\prime}(r)$ is found by removing constant terms from $p(r)$, the measured absorption energy map as a photoacoustic image, which is the product of AC and LF, which can be written as Eq. (2):
Where ${\mu _a}(r)$ and $F(r,{\mu _{s^{\prime}}}(r))$ are the AC and LF respectively. Here, we assume that the temperature, scattering, and other parameters are known, and the absorption energy map and the raw PA map are proportional.Here, we will introduce a new quantitative photoacoustic imaging algorithm based on signal sparsity decomposition. This method relies on differences in the spatial properties of the AC and optical energy density in the medium. The method performs signal decomposition based on sparse representation theory, using these different features to extract the AC and LF of the imaged object from the photoacoustic image. The algorithm schematic is shown in Fig. 2. For quantitative PA imaging, we aim to recover AC by disentangling it from LF. Traditional methods assume homogeneous fluence, but in practice, LF exhibits depth-dependent attenuation. Our approach leverages ultrasound-guided sparse decomposition to address this challenge.
Prior knowledge of the structure was incorporated into the Monte Carlo simulation process to obtain more accurate absorption and LF distributions. A dictionary learning algorithm, K-SVD [35], is introduced that generalizes the clustering of k-means and efficiently learns an over-complete dictionary from a set of training signals. Based on AC and LF to train the dictionary separately, the K-SVD algorithm is applied to the process of finding the best dictionary, which can sparsely represent each signal.
Structural priors from co-registered US images constrain MC simulations to generate anatomically realistic optical property maps. US-defined tissue boundaries (e.g., vessels, muscle layers) are segmented to create spatially resolved AC and LF maps. Trained on AC patches from MC simulations, capturing sharp, localized absorption features. Trained on LF patches, encoding smooth, depth-dependent attenuation patterns.
K-SVD generalizes k-means clustering by learning an overcomplete dictionary D from training signals. For AC and LF dictionaries, the optimization is formulated as:
Then, the dictionaries with absorption distribution characteristics and light propagation characteristics are formed into a joint dictionary as input to the sparse decomposition process. The other input is the raw photoacoustic signal, which represents the measured data. The output of the orthogonal matching pursuit (OMP) algorithm is a matrix composed of sparse vectors of coefficients, and after extracting the sparse vector corresponding to the absorption coefficient part, the result calculated in the dictionary of the corresponding AC characteristics is the absorption distribution of the final quantitative extraction.
Sparse representation theory is used to decompose a photoacoustic image into two components. The sparse representation in this article can be defined as that the photoacoustic image can be sparsely expressed based on the appropriate dictionary or library, and more obviously, the input image or the photoacoustic signal can be represented as the sum of a small number of basic functions [16,17,34,35,41,43,44]. Moreover, the sparse representation is adaptive to select the appropriate super-complete dictionary according to the characteristics of the photoacoustic signal. The purpose of sparse representations of photoacoustic images is to find an adaptive dictionary to make the input signal the most sparsely expressed. To match the use of sparse representation theory, we convert the signal from the product formula to the summation formula., which can be written as Eq. (4):
The two component functions $\log [{\mu _\textrm{a}}(r)]$ and $\log [F(r)]$ have very distinguishable characteristics. The feature changes corresponding to the absorption distribution are sharp and locally changing; in contrast, the feature changes corresponding to LF are slow and globally attenuated. Extraction of AC and LF extraction based on different features.
Two important steps to achieve signal sparse representation of photoacoustic images are the generation of dictionaries and the sparse decomposition of signals. For the choice of dictionaries, we use a learning dictionary with strong adaptive ability, compared with the use of fixed dictionaries, the advantage is that it can better adapt to data with different characteristics. Our decomposition of $\log [p^{\prime}(r)]$ is based on the assumption that we use the K-SVD algorithm to design two libraries ${\phi _n}$ and ${\varphi _m}$ which have rapid change characteristics and global attenuation characteristics, respectively. We will refer to these conditions as the sparsity conditions. When the sparsity conditions are met, the photoacoustic image can be sparsely represented in the joint library $D = {\{{{\phi_n},{\varphi_m}} \}_{n,m}}$, which can be written as Eq. (5):
In our work, we use the OMP algorithm to improve the accuracy of the decomposition process. In OMP, the optimal coefficients for all elements are recalculated with each iteration, allowing errors in earlier iterations to be compensated for. OMP converges faster than traditional MP algorithms. Here, a Monte Carlo model of light transport is utilized, which simulates initialization data, including the measured PA map and the LF distribution of the pre-set model. The corresponding dictionary is trained according to the characteristics of AC and LF, and the distribution of AC is obtained based on the sparse decomposition algorithm. The sparse decomposition based on the OMP algorithm is shown in Table 2.
Table 2. Sparse decomposition algorithm.
In the process of sparse decomposition using Orthogonal Matching Pursuit (OMP), we focus on efficiently extracting the sparse representation of the photoacoustic signal based on the previously constructed joint dictionary. When executing the OMP algorithm, the input will be sparsely represented by the photoacoustic vector as X, the joint dictionary D and the maximum number of coefficients for each signal. Initialize the residual matrix and the sparse coefficient vector $\alpha$. The corresponding dictionary elements are extracted by calculating the maximum inner product of the residual matrix and the joint dictionary. Finally, the sparse coefficient vector is calculated by the Moore-Penrose pseudoinverse equation, which can be written as Eq. (6):
Finally, the signal decomposition is performed, and the data corresponding to the absorption distribution are extracted for the sparse coefficient vector calculated in the Eq. (5), the reconstruction of LF distribution ${F^e}$ is carried out based on dictionary${\varphi _m}$ with LF distribution characteristics, and the absorption distribution $\mu _a^e$ is reconstructed based on ${F^e}$ with LF compensation. The parameter δ is a regularization factor that ensures numerical stability and suppresses noise during the sparse coefficient recovery. It is typically a constant [11,12,46,47]. AC calculation process can be written as Eq. (7):
3. Experiments and results
3.1. Simulation results of spatial characteristic distribution
Monte Carlo simulation is used here to simulate tissue light propagation models, including vascular and tumor imaging targets. The whole geometry is a cube of the same width and depth, where the width and depth are 0.512 cm, and the grid matrix is 512 × 512. The z-x plane configuration for the simulation model is shown in Fig. 3. Figure 3 shows the design of three different simulation models with increasing structural complexity, as seen in Fig. 3. Light transport simulations were conducted using the Monte Carlo method, resulting in the corresponding LF distribution for each model, as shown in Fig. 3. Since LF attenuates with increasing depth, it affects the quality of the PA imaging, as illustrated in Fig. 3. Table 3 shows the simulation model parameters.
Fig. 3. Spatial distribution of AC, LF, and PA in three simulation models. (a-c) AC, LF, and PA distributions in a simple vascular structure model. (d-f) AC, LF, and PA distributions in a complex vascular structure model. (g-i) AC, LF, and PA distributions in a more complex vascular structure model.
Table 3. Simulation model parameters
The quantitative photoacoustic imaging algorithm is based on the theory of signal sparsity decomposition. This method relies on differences in the spatial properties of the absorption coefficient and optical energy density in the medium. The K-SVD algorithm is used to find adaptive dictionaries so that the signal is expressed in its most sparse form. The OMP algorithm is used to perform sparse decomposition and finally accurately extract the AC and LF distributions. To decompose the photoacoustic image, two libraries should be found that meet the sparsity criteria. These functions depend on the absorption coefficient and the spatial characteristics of the fluence. Figure 4 shows that the spatial properties of the absorption coefficient and fluence are distinguishable.
Fig. 4. Spatial intensity variation of AC, LF, and PA in three simulation models. (a-c) Spatial intensity distribution of AC, LF, and PA in the simple vascular structure model. (d-f) Spatial intensity distribution of AC, LF, and PA in the medium-complexity vascular structure model. (g-i) Spatial intensity distribution of AC, LF, and PA in the complex vascular structure model.
Figure 4 presents the spatial intensity variation corresponding to the results shown in Fig. 3. In the three models, the spatial distribution of the absorption coefficient exhibits sharp changes, while the LF shows a global attenuation trend, gradually decreasing with increasing depth. The photoacoustic images are influenced by the combined effects of both AC and LF.
Figure 5 illustrates the spectral intensity variation of the AC and LF corresponding to the results shown in Fig. 3. The figure reveals the unique spectral characteristics of AC and LF:The AC exhibits multiple frequency components, reflecting its complex behavior influenced by heterogeneous tissue structures (e.g., blood vessels, tumors). The LF is dominated by low-frequency components, consistent with the smooth photon diffusion in turbid media.These findings are crucial for understanding the optical properties of biological tissues under various conditions. The spectral distribution changes across the three models further confirm the distinct frequency-domain characteristics of AC and LF, providing a foundation for quantitative photoacoustic imaging technology based on sparse representation.
Fig. 5. Spectral intensity variation of AC, LF, and PA in three simulation models. (a-c) Spectral intensity distribution of AC, LF, and PA in the simple vascular structure model. (d-f) Spectral intensity distribution of AC, LF, and PA in the medium-complexity vascular structure model. (g-i) Spectral intensity distribution of AC, LF, and PA in the complex vascular structure model.
3.2. Simulation experiments and results
For simulation experiments, the Monte Carlo model is used to model tissue light propagation, simulating targets such as blood vessels and tumor imaging. In the simulation experiment, the mouse torso model and its corresponding PA image are shown in Fig. 6. The raw PA image exhibits structural and detail loss, as demonstrated in Fig. 6. The simulation model shown in Fig. 6(a) was designed to simulate the mouse trunk. The entire geometry is a cube of the same width and depth, with a width and depth of 0.512 cm and a grid matrix of 512 × 512. The Z-X plane configuration of the simulation model is shown in Fig. 6. The absorption coefficient of the absorption model is 100 cm-1, the attenuation scattering coefficient is 100 cm-1, and the anisotropy factor g is 0.9. A uniform laser beam of 0.50 cm in diameter is used to illuminate the central imaging region. The laser wavelength is 532 nm and the simulation time is 1.0 minutes.
Fig. 6. Mouse torso simulation model and corresponding PA image. (a) Mouse torso simulation model, including Mouse torso, Spine, Liver, Intestine, Abdominal Aorta (AA), Portal Vein (PV), and Inferior Vena Cava (IVC). (b) PA image.
Our simulation experiment only considers the optical inversion problem which needs to be solved in this paper, but does not consider the acoustic inversion problem. The original PAT image for the simulation experiment is obtained by adding about 40 dB noise to the product of the ideal AC image and the estimated LF distribution. In the simulation experiment, structural information from the tissue model is used instead of US imaging modality for structural priors to initialize the Monte Carlo model. By exploiting the distinct spatial variation characteristics of LF and AC, a sparse decomposition method is applied to extract the LF distribution from the raw PA data. This process ultimately compensates the LF to achieve the reconstructed AC distribution.
Figure 7 presents a detailed comparison between ideal AC images, reconstructed AC distributions, raw PA images, and LF distributions. Figure 7 (a) shows the ground truth AC image, serving as the reference for assessing the accuracy of both the raw PA image and the reconstructed AC distribution. By comparing this ideal image with others, we can evaluate the effectiveness of the reconstruction algorithm. Figure 7 (b) provides structural information, likely derived from simulations or other imaging modalities, which helps guide the reconstruction process. This structural map may influence the initialization of models, such as Monte Carlo simulations and sparse decomposition, which aid in LF extraction and AC reconstruction. Figure 7(c) shows the raw PA image, which contains both AC and LF information. Without proper separation techniques, this image doesn't directly reveal the AC distribution. For better analysis, the raw PA data is typically normalized and filtered. Figure 7 (d) displays the reconstructed AC distribution. Using methods like sparse decomposition, the LF effects are separated from the raw PA image, providing a more accurate AC representation, which is comparable to the ideal AC image. Figure 7(e) and (f) show LF distributions obtained from Monte Carlo simulations and sparse decomposition, respectively. The successful separation of LF is critical for accurate AC reconstruction. Figure 7(g) and (h) present error maps that compare the raw PA and reconstructed AC images to the ideal AC image. These maps quantify the deviation, with lower errors in Fig. 7 (h) demonstrating the success of the reconstruction method. Overall, Fig. 7 demonstrates the effectiveness of sparse decomposition in separating LF and accurately reconstructing the AC distribution, resulting in significantly reduced errors compared to raw PA data.
Fig. 7. Simulation results. (a) Ideal AC image. (b) Structural information map. (c) Raw PA image. (d) Reconstructed AC distribution. (e) Initial LF distribution from Monte Carlo model. (f) LF distribution extracted using sparse decomposition. (g) Error map between raw PA and ideal AC. (h) Error map between reconstructed AC and ideal AC.
Figure 8 presents the visualization results of the dictionaries, including the dictionaries for LF and AC, as well as their joint visualization. Figure 8 (a) displays the dictionary associated with LF distribution, where the dictionary elements represent variations in light intensity under different spatial characteristics, aiding in the accurate extraction of LF information in complex tissues. Figure 8 (b) shows the dictionary related to AC distribution, with its elements capturing the absorption properties of tissues, thereby enhancing the understanding of tissue structures. Finally, Fig. 8 (c) provides a joint visualization of the LF and AC dictionaries, offering a view of the relationship between the two. This joint representation not only strengthens the awareness of the relationship between light and absorption features in the model but also provides important structural priors for subsequent reconstruction processes, facilitating more accurate photoacoustic imaging results.
Fig. 8. Dictionary visualization results. (a) Dictionary for LF. (b) Dictionary for AC. (c) Joint Dictionary Visualization.
Figure 9 presents a quantitative comparison of reconstruction results, focusing on spatial intensity variations and LF distribution models. In Fig. 9(a), the horizontal intensity profiles of the ideal AC, reconstructed AC, and raw PA images are compared. The reconstructed AC should closely follow the ideal AC, indicating successful reconstruction, while the raw PA deviates due to unseparated AC and LF components. Figure 9(b) shows the vertical intensity variations. Depth-wise accuracy is critical for maintaining the internal structure, and the reconstructed AC should match the ideal AC at different depths. Any discrepancies between the raw PA and reconstructed AC show LF interference before separation. Figure 9(c) illustrates the initial LF distribution model, which serves as a baseline for separating LF from the raw PA image. The success of the separation algorithm depends on the accuracy of this model. Figure 9(d) displays the LF distribution extracted by the separation algorithm. Comparing this to the initial model shows how well the algorithm isolates LF, improving the reconstruction of the AC.
Fig. 9. Quantitative reconstruction results comparison. (a) Comparison of horizontal spatial intensity variations for ideal AC, reconstructed AC, and raw PA. (b) Comparison of vertical depth-wise spatial intensity variations for ideal AC, reconstructed AC, and raw PA. (c) Initial LF distribution model. (d) LF distribution from separation algorithm output.
3.3. Phantom experiments and results
To further validate the feasibility of the proposed method, we tested its performance using a tissue-mimicking phantom with embedded absorbers. The phantom's background consisted of 2% Intralipid and 1% agar powder. Three tubes, each filled with a 2% mass fraction of Indian ink as absorbers, were embedded at 5 mm depth intervals within the phantom. Each tube had a maximum diameter of 3 mm, as depicted in Fig. 10.
In Fig. 11, we present the experimental results from the tissue-mimicking phantom using the dual-modal PAUS technique. Figure 11 (a) illustrates the quantitative reconstruction of fused imaging results from US and PA data prior to any reconstruction process. This initial fusion provides insight into the spatial distribution of both imaging modalities. Figure 11 (b) further compares the reconstructed imaging results, demonstrating the effectiveness of the proposed reconstruction method in integrating information from both US and PA modalities, highlighting improvements in image quality and accuracy following reconstruction.
Fig. 11. Tissue-mimicking phantom PAUS dual-modal experimental results. (a) Quantitative reconstruction of US and PA fused imaging results before reconstruction. (b) Comparison of quantitative reconstruction results of US and PA fused imaging.
Table 4 presents a comparison of size calculations for three imaging targets based on US, raw PA, and reconstructed AC data. The size calculated from the US data exhibited a deviation of 0.6% larger than the actual diameter of 3 mm. Conversely, the size calculated from the raw PA data showed a deviation of 12.5% smaller than the true size. The reconstructed AC data yielded a size calculation error of 3%. These results indicate that US provides reliable structural information, while the raw PA data is affected by noise and artifacts, leading to distorted structural information. However, the results from the quantitative reconstruction demonstrate notable improvements.
Table 4. Phantom diameter parameters
The PA signals within the tube structure were compensated through co-registration based on the structural and boundary information provided by US image. This aims to reflect the true absorption distribution. The US image provides spatial constraints to exclude artifact-contaminated regions during training patch extraction. Training samples, derived from a Monte Carlo simulation model, are selectively extracted from US-defined homogeneous regions, ensuring that the learned dictionary atoms encode true absorption features rather than boundary artifacts. The US image delineates the approximate boundaries of the tube structure. The spatial alignment of US and PA ensures that the simulated LF distribution (via Monte Carlo) matches the true geometry, thereby reducing depth-dependent fluence errors. The trained dictionaries are shown in Fig. 12, where Fig. 12(a-b) present dictionaries with different features for the LF and the AC, respectively. Based on the theory of sparse decomposition, further separation of the absorption coefficient and LF is performed.
Fig. 12. Dictionary visualization. (a) Dictionary for LF. (b) Dictionary for AC. (c) Joint Dictionary.
In Fig. 13, we present the quantitative experimental results from the tissue-mimicking phantom using PA imaging. Figure 13 (a-c) showcase the PA imaging results for targets 1 to 3 prior to any reconstruction, highlighting the inherent features captured by the raw data. Figure 13 (d-f) display the PA imaging results after quantitative reconstruction for the same targets, indicating improvements in clarity and detail achieved through the reconstruction process. Finally, Fig. 13(g-i) illustrate the corresponding LF results for targets 1 to 3, Fig. 13(j-k) illustrate the corresponding log-LF results for targets 1 to 3, allowing for a comprehensive comparison of the different imaging techniques applied to the phantom.
Fig. 13. Tissue-mimicking phantom PA quantitative experimental results. (a-c) Raw PA imaging results for targets 1-3. (d-f) Quantitative reconstruction AC results for targets 1-3. (g-i) Corresponding LF results for targets 1-3. (j-k) Corresponding log-LF results for targets 1-3.
In Fig. 14, we display the results of spatial variations in both horizontal and vertical dimensions for targets 1 to 3. Figure 14 (a-c) compare the horizontal spatial intensity variations, illustrating the differences between raw PA data and reconstructed AC results. This comparison emphasizes the improvements in spatial resolution and intensity rendition achieved through the reconstruction process. Figure 14 (d-f) present a similar analysis for vertical depth-wise spatial intensity variations, enabling a thorough examination of the performance of the reconstruction technique in different spatial orientations.
Fig. 14. Horizontal and vertical spatial variation results. (a-c) Comparison of horizontal spatial intensity variations for targets 1-3, showing raw PA and reconstructed AC results. (d-f) Comparison of vertical depth-wise spatial intensity variations for targets 1-3, illustrating raw PA and reconstructed AC results.
3.4. Animal experiments and results
In Fig. 15, we present the PAUS dual-modality quantitative imaging results for a mouse hind-limb. The US map primarily provides information about the muscles and bones of the mouse hind-limb, while the PA map showcases the functional information of a major blood vessel and two branching vessels extending from the mouse's hindquarters down through the thigh. Mouse hind-limb experimental images and a schematic of PAUS imaging are shown in Fig. 15(a) . The reference length of the mouse hind limb (from hip to ankle) was approximately 25 mm, and we imaged the region about 22 mm in length. Figure 15 (b) shows the unprocessed PA image, capturing initial signals from hind limb tissues. Vascular structures like the Popliteal Artery (PAr), Femoral Vein (FV), and Femoral Artery (FA) are visible, but clarity and contrast are low, making smaller structures hard to delineate. Figure 15 (c) presents the reconstructed image, enhancing anatomical features. Vascular structures such as PAr, PV, FA, and FV are clearer, with better resolution and contrast. This highlights the AC method's effectiveness in isolating and enhancing blood vessels. Figure 15 (d) and (e) display fusion maps of PA and US images. Based on the US image in Fig. 15(d), Fig. 15(e) enhances the visualization of the hind-limb muscles (HM) and bones (HB) in the mouse, providing anatomical boundaries (e.g., muscle-bone interfaces) that serve as structural priors to guide sparse decomposition. Furthermore, by separating absorption features from fluence artifacts, the reconstructed absorption coefficients (Fig. 15(c)) align more precisely with the US-defined anatomy. The fused PA-US image (Fig. 15(e)) overlays the processed PA absorption map (artifact-reduced) onto the US image, highlighting co-localized features. The fusion integrates anatomical context from US data with functional imaging from PA, improving clarity and localization of features like Hind-limb Muscles (HM) and Bones (HB). Figure 15 (f) likely represents a structural model of the hind limb's key anatomical elements, providing a reference for comparing raw PA and reconstructed AC images and interpreting spatial relationships. Figure 15 (g) shows the LF distribution map, offering insight into LF within tissues. Understanding this is crucial for accurate quantitative PA imaging, affecting signal strength and tissue property reconstruction. To effectively represent the signal intensity at each position and better illustrate spatial variations, we extracted the maximum intensity values along the profile lines. Figure 15 (h) and (i) plot spatial intensity variations horizontally and vertically, comparing raw PA and reconstructed AC results. Both directions show improved signal uniformity and intensity post-reconstruction, enhancing clarity and accuracy, vital for assessing tissue depth and blood vessel distribution. In summary, Fig. 13 shows that PAUS dual-modal imaging enhances vascular signal detection, spatial resolution, and quantitative accuracy, with an average signal intensity increase of 38% to 75%.
Fig. 15. Mouse hind-limb PAUS dual-modality quantitative imaging results. (a) Experimental image of mouse hind-limb. (b) Raw PA image. (c) Quantitative reconstructed AC image. (d) Fusion map of raw PA and US. (e) Fusion map of reconstructed AC and US. (f) Structural information model. (g) LF distribution. (h) Horizontal spatial variation curves for raw PA and reconstructed AC. (i) Vertical depth-wise spatial variation curves for raw PA and reconstructed AC. Description of markers: Popliteal Artery (PAr),Popliteal Vein(PV),Femoral Vein(FV),Femoral Artery(FA),Hind-limb Muscles (HM), Hind-limb bones(HB)
We calculated the spatial overlap between the reconstructed photoacoustic (PA) image and the ultrasound (US) vessel mask. The Dice coefficient improved from 0.62 for the raw PA image to 0.85 for the reconstructed image, indicating a significant enhancement in the alignment of vessel contours. The diameter of the reconstructed femoral artery ranged from 0.35 to 0.78 mm, which is consistent with the literature-reported range of 0.3–0.8 mm [42,48], further confirming the physiological plausibility of the reconstruction results. Additionally, the peak signal-to-noise ratio (PSNR) of the reconstructed image improved from 17.46 dB to 30.18 dB, demonstrating a significant enhancement in signal quality relative to noise. The structural similarity index (SSIM) increased from 0.67 to 0.92, indicating a substantial improvement in the structural similarity between the reconstructed and original images. The gradient magnitude at vessel boundaries improved by 22%, suggesting sharper and more distinct edges. PAr is smaller, about 0.2 to 0.5 mm in diameter. FV is generally slightly larger than the artery and is about 0.5 to 1.0 mm in diameter. The PV diameter is usually between 0.3 and 0.7 mm. By using the proposed PAUS dual-mode quantitative imaging method, the improvement of the inner structure and edge information of the blood vessel is achieved by 11.4% to 33.3% . The vascular structural parameters of mice are shown in Table 5.
Table 5. Mouse hind-limb blood vessels structural parameters
Figure 16 presents the visualization results of the dictionaries related to Fig. 15, further revealing the features extracted during the photoacoustic and ultrasound imaging processes. By analyzing these dictionaries, a better understanding of the specific manifestations and interrelationships of light fluence and absorption coefficient in the imaging of the mouse hind-limb can be achieved.
Fig. 16. Dictionary visualization. (a) Dictionary for LF. (b) Dictionary for AC. (c) Joint Dictionary.
Figure 17 presents a quantitative comparison of PA results across three regions of interest (ROIs), focusing on different vascular structures. It includes both raw PA images and reconstructed AC. Figure 17 (a, d, g) show visual comparisons of the reconstructed vascular structures in the 3 ROIs, highlighting areas such as the FA, FV, PAr, and PV. Figure 17 (b, e, h) depict the horizontal spatial intensity variations, comparing raw PA data (red) with AC-reconstructed data (blue). The AC data demonstrate generally enhanced intensity, indicating improved contrast and resolution of vascular features. Figure 17 (c, f, i) illustrate vertical spatial intensity variations, revealing depth-wise intensity profiles. The AC-reconstructed curves show better definition and higher peak intensities, suggesting a more accurate representation of depth-related structures within the ROIs. Overall, the comparison emphasizes that the AC-reconstructed data provide superior spatial resolution and contrast in both horizontal and vertical profiles, enhancing the clarity of vascular structures compared to raw PA images.
Fig. 17. Comparison of quantitative PA results for 3 ROIs. (a, d, g) Comparison of quantitative reconstruction results for vascular structures for ROIs 1-3. (b, e, h) Comparison of horizontal spatial intensity variation curves between raw PA and reconstructed AC. (c, f, i) Comparison of vertical depth-wise spatial intensity variation curves between raw PA and reconstructed AC.
Table 6 presents a quantitative comparison of absorption coefficients between raw PA and reconstructed AC images across three regions of interest (ROIs) in the mouse hind-limb. The results demonstrate significant improvements in absorption coefficient values post-reconstruction, with enhancements ranging from 31.34% to 64.86%. Specifically, in ROI 1, the absorption coefficient increased from 0.74 ± 0.16 mm−1 (raw PA) to 1.22 ± 0.15 mm−1 (reconstructed AC), representing the most substantial improvement. Similarly, ROI 2 and ROI 3 showed increases from 0.41 ± 0.17 mm−1 to 0.59 ± 0.13 mm−1 and from 0.67 ± 0.08 mm−1 to 0.88 ± 0.06 mm−1, respectively. Notably, the standard deviations of the reconstructed absorption coefficients are consistently lower than those of the raw PA data, indicating improved precision and reduced variability. These findings underscore the effectiveness of the proposed reconstruction method in enhancing both the accuracy and reliability of absorption coefficient measurements, which is critical for quantitative PA imaging applications.
Table 6. Mouse hind-limb blood vessels absorption parameters comparison
4. Discussion and conclusions
In this study, we proposed an innovative method for quantitatively extracting absorption distributions in tissue imaging, particularly in cases where the light fluence is non-uniform in deeper tissues. By applying a sparse representation framework, we successfully extracted the light fluence from photoacoustic images. Through LF compensation, we were able to accurately and quantitatively extract the absorption distribution. Our approach integrates photoacoustic tomography with ultrasound dual-modality imaging, leveraging ultrasonic segmentation for structural guidance and enhancement. This dual-modality setup provides essential structural prior knowledge for the construction of Monte Carlo optical transport models, which further refines the reconstruction process by correcting for non-uniform LF and absorption effects.
Our experimental results demonstrate significant improvements in both feature visibility and quantitative accuracy compared to conventional qPAT techniques. For instance, the Dice coefficient increased from 0.62 to 0.85, indicating better alignment with ultrasound-defined anatomy. Additionally, the absorption coefficient precision improved by up to 31.34%, highlighting the effectiveness of our sparse decomposition framework. Through numerical simulations and experiments with tissue-mimicking phantoms and in vivo mouse models, our results demonstrate significant improvements in both structural and functional information. Furthermore, our method enhances feature visibility, boundary definition, and overall image quality. The advancements presented in this paper hold significant promise for quantitative photoacoustic imaging applications, such as vascular system and tumor detection, where precise absorption distribution estimation is crucial. When compared to deep learning-based methods, our approach offers several unique advantages. First, the interpretability of our method allows for transparent and physically meaningful reconstructions, which is critical for clinical applications. Second, our framework achieves robust performance with minimal training data by leveraging MC simulations and ultrasound priors, addressing the data efficiency challenges faced by deep learning models. Finally, the generalizability of our method across different tissue types and imaging geometries makes it a versatile tool for quantitative PA imaging.
Despite its advantages, our method has some limitations that need to be acknowledged. The sparse decomposition process, particularly the orthogonal matching pursuit (OMP) algorithm, can be computationally intensive, especially for high-resolution images. Additionally, MC simulations, while accurate, are computationally expensive for complex tissue geometries. To mitigate this, we plan to optimize the algorithm implementation using parallel computing or GPU acceleration in future work. Another limitation is the sensitivity to noise, particularly in low signal-to-noise ratio (SNR) conditions, which can degrade the accuracy of sparse decomposition. We intend to incorporate noise suppression techniques, such as wavelet-based denoising or regularization strategies, to improve robustness. Furthermore, our method relies on ultrasound images for structural guidance, which may limit its applicability in cases where ultrasound data is unavailable or of poor quality. To address this, we aim to develop more robust registration algorithms and explore alternative approaches that do not require ultrasound priors.
In future research, we plan to synergize physics-informed neural networks (PINNs) with multi-wavelength photoacoustic datasets to optimize optical parameters, mitigate the impact of noise and experimental parameter variations, and achieve deeper tissue imaging [37–40,42,48,49]. While the current method focuses on addressing the optical inverse problem under heterogeneous absorption conditions, future work will integrate full-wave acoustic inversion to account for realistic acoustic effects, such as heterogeneous sound speed and detector geometry [46,47,50–52]. This will involve coupling our sparse decomposition framework with advanced acoustic solvers, such as k-Wave or time-reversal algorithms, to enable end-to-end quantitative photoacoustic tomography.
Funding
National Key Research and Development Program of China (Grant No. 2022YFC2402400); Shenzhen Medical Research Fund (Grant No. D2404002); National Natural Science Foundation of China (Grant No. 62275062); Project of Shandong Innovation and Startup Community of High-end Medical Apparatus and Instruments (2023-SGTTXM-002 and 2024-SGTTXM-005); the Shandong Province Technology Innovation Guidance Plan (Central Leading Local Science and Technology Development Fund) (Grant No. YDZX2023115); the Taishan Scholar Special Funding Project of Shandong Province, and the Shandong Laboratory of Advanced Biomaterials and Medical Devices in Weihai (Grant No. ZL202402).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data may be obtained from the authors upon reasonable request.
References
1. L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). [CrossRef]
2. M. L. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 305–598 (2006). [CrossRef]
3. P. Beard, “Biomedical photoacoustic imaging,” Interface Focus. 1(4), 602–631 (2011). [CrossRef]
4. X. Wang, Y. Pang, L. V. Ku, et al., “Noninvasive laser-induced photoacoustic tomography for structural and functional in vivo imaging of the brain,” Nat. Biotechnol. 21(7), 803–806 (2003). [CrossRef]
5. L. V. Wang, H. Wu, and B. R. Masters, “Biomedical optics: Principles and imaging,” J. Biomed. Opt. 13(4), 049902 (2008). [CrossRef]
6. S. L. Jacques, “Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation,” Photoacoustics 2(4), 137–142 (2014). [CrossRef]
7. S. Tzoumas, A. Nunes, I. Olefir, et al., “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nat. Commun. 7(1), 12121–10 (2016). [CrossRef]
8. W. Choi, D. Oh, and C. Kim, “Practical photoacoustic tomography: realistic limitations and technical solutions,” J. Appl. Phys. 127(23), 230903 (2020). [CrossRef]
9. B A Kaplan, J Buchmann, S Prohaska, et al., “Monte-Carlo-based inversion scheme for 3D quantitative photoacoustic tomography,” Photons Plus Ultrasound: Imaging and Sensing 2017, 10064, 802–814 (SPIE, 2017).
10. B. Cox, J. G. Laufer, S. R. Beard, et al., “Quantitative spectroscopic photoacoustic imaging: A review,” J. Biomed. Opt. 17(6), 061202 (2012). [CrossRef]
11. B T Cox, S R Arridge, P C. Kostli, et al., “Quantitative photoacoustic imaging: Fitting a model of light transport to the initial pressure distribution,” Photons Plus Ultrasound: Imaging & Sensing 2005: The Sixth Conference on Biomedical Thermoacoustics, Optoacoustics, and Acousto-Optics5697, 49–56 (2005).
12. B. T. Cox, S. R. Arridge, P. C. Köstli, et al., “Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,” Appl. Opt. 45(8), 1866–1875 (2006). [CrossRef]
13. B. T. Cox, S. R. Arridge, and P. C. Beard, “Gradient-based quantitative photoacoustic image reconstruction for molecular imaging,” Photons Plus Ultrasound: Imaging and Sensing 2007: The Eighth Conference on Biomedical Thermoacoustics, Optoacoustics, and Acousto-optics, 6437: 445–454 (2007).
14. B. Banerjee, S. Bagchi, D. Vasu, et al., “Quantitative photoacoustic tomography from boundary pressure measurements: Noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,” J. Opt. Soc. Am. A 25(9), 2347–2356 (2008). [CrossRef]
15. R. J. Zemp, “Quantitative photoacoustic tomography with multiple optical sources,” Appl. Opt. 49(18), 3566–3572 (2010). [CrossRef]
16. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Quantitative photoacoustic signal extraction using sparse signal representation,” IEEE Trans. Med. Imaging 28(12), 1997–2006 (2009). [CrossRef]
17. A Rosenthal, D Razansky, and V. Ntziachristos, “Sparse signal representation at the service of quantitative photoacoustic tomography,” In Photons Plus Ultrasound: Imaging & Sensing. International Society for Optics and Photonics; International Society for Optics and Photonics, 7564: 75640S (2010).
18. P. Shao, B. Cox, and R. J. Zemp, “Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,” Appl. Opt. 50(19), 3145–3154 (2011). [CrossRef]
19. A. Q. Bauer, R. E. Nothdurft, J. P. Culver, et al., “Quantitative photoacoustic imaging: correcting for heterogeneous light fluence distributions using diffuse optical tomography,” J. Biomed. Opt. 16(9), 096016 (2011). [CrossRef]
20. S. Na and L. V. Wang, “Photoacoustic computed tomography for functional human brain imaging [Invited],” Biomed. Opt. Express 12(7), 4056–4083 (2021). [CrossRef]
21. M. A. Mastanduno and S. S. Gambhir, “Quantitative photoacoustic image reconstruction improves accuracy in deep tissue structures,” Biomed. Opt. Express 7(10), 3811–3825 (2016). [CrossRef]
22. F. M. Brochu, J. Brunker, J. Joseph, et al., “Towards quantitative evaluation of tissue absorption coefficients using light fluence correction in optoacoustic tomography,” IEEE Trans. Med. Imaging 36(1), 322–331 (2017). [CrossRef]
23. F. Asllanaj and A. Addoum, “Simultaneous reconstruction of absorption, scattering and anisotropy factor distributions in quantitative photoacoustic tomography,” Biomed. Phys. Eng. Express 6(4), 045010 (2020). [CrossRef]
24. J. Buchmann, B. Kaplan, S. Powell, et al., “Quantitative PA tomography of high resolution 3-D images: experimental validation in a tissue phantom,” Photoacoustics 17, 100157 (2020). [CrossRef]
25. C. Yang, H. Lan, H. Zhong, et al., “Quantitative photoacoustic blood oxygenation imaging using deep residual and recurrent neural network,” 2019 IEEE 16th International Symposium on Biomedical Imaging (IEEE, 2019), pp. 741–744.
26. S. Bu, Z. Liu, Y. Shiina, et al., “Model-based reconstruction integrated with fluence compensation for photoacoustic tomography,” IEEE Trans. Biomed. Eng. 59(5), 1354–1363 (2012). [CrossRef]
27. Z. Yuan, Q. Wang, and H. Jiang, “Reconstruction of optical absorption coefficient maps of heterogeneous media by photoacoustic tomography coupled with diffusion equation based regularized Newton Method,” Opt. Express 15(26), 18076–18081 (2007). [CrossRef]
28. Z. Yuan and H. Jiang, “Three-dimensional finite-element-based photoacoustic tomography: Reconstruction algorithm and simulations,” Med. Phys. 34(2), 538–546 (2007). [CrossRef]
29. S. Zhang, J. Liu, Z. Liang, et al., “Pixel-wise reconstruction of tissue absorption coefficients in photoacoustic tomography using a non-segmentation iterative method,” Photoacoustics 28, 100390 (2022). [CrossRef]
30. Z. Liang, S. Zhang, J. Wu, et al., “Automatic 3-D segmentation and volumetric light fluence correction for photoacoustic tomography based on optimal 3-D graph search,” Med. Image Anal. 75, 102275 (2022). [CrossRef]
31. Y. Liu, H. Jiang, and Z. Yuan, “Two schemes for quantitative photoacoustic tomography based on Monte Carlo simulation,” Med. Phys. 43(7), 3987–3997 (2016). [CrossRef]
32. M. Kirillin, V. Perekatova, P. Turchin, et al., “Fluence compensation in raster-scan photoacoustic angiography,” Photoacoustics 8, 59–67 (2017). [CrossRef]
33. A. Dolet, R. Ammanouil, V. Petrilli, et al., “In vitro and in vivo multispectral photoacoustic imaging for the evaluation of chromophore concentration,” Sensors 21(10), 3366 (2021). [CrossRef]
34. M. Elad and A. M. Bruckstein, “A generalized uncertainty principle and sparse representation in pairs of bases,” IEEE Trans. Inf. Theory 48(9), 2558–2567 (2002). [CrossRef]
35. M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” IEEE Trans. Signal Process. 54(11), 4311–4322 (2006). [CrossRef]
36. R. Hochuli, S. Powell, S. Arridge, et al., “Quantitative photoacoustic tomography using forward and adjoint Monte Carlo models of radiance,” J. Biomed. Opt. 21(12), 126004 (2016). [CrossRef]
37. S. M. Ranjbaran, H. S. Aghamiry, A. Gholami, et al., “Quantitative photoacoustic tomography using iteratively refined wavefield reconstruction inversion: a simulation study,” IEEE Trans. Med. Imaging 43(2), 874–885 (2024). [CrossRef]
38. R. Zhang, R. A’dawiah, T. W. J. Choo, et al., “Navigating challenges and solutions in quantitative photoacoustic imaging,” Appl. Phys. Rev. 11(3), 031308 (2024). [CrossRef]
39. J Gröhl, T R Else, L Hacker, et al., “Moving beyond simulation: data-driven quantitative photoacoustic imaging using tissue-mimicking phantoms,” IEEE Transactions on Medical Imaging 43, 1214 (2023). [CrossRef]
40. T. Tarvainen and B. Cox, “Quantitative photoacoustic tomography: modeling and inverse problems,” J. Biomed. Opt. 29(S1), S11509 (2023). [CrossRef]
41. Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition,” Proceedings of 27th Asilomar Conference on Signals, Systems and Computers (IEEE, 1993), 40–44.
42. M. Li, M. Beamont, C. Ma, et al., “. Three-dimensional deep-tissue functional and molecular imaging by integrated photoacoustic, ultrasound, and angiographic tomography (PAUSAT),” IEEE Transactions on Medical Imaging 41(10), 2704–2714 (2022). [CrossRef]
43. J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Trans. Inf. Theory 53(12), 4655–4666 (2007). [CrossRef]
44. T. T. Cai and L. Wang, “Orthogonal matching pursuit for sparse signal recovery with noise,” IEEE Trans. Inf. Theory 57(7), 4680–4688 (2011). [CrossRef]
45. F. Liu, X. Gong, L. V. Wang, et al., “Dictionary learning sparse-sampling reconstruction method for in-vivo 3D photoacoustic computed tomography,” Biomed. Opt. Express 10(4), 1660–1677 (2019). [CrossRef]
46. J. Prakash, D. Sanny, S. K. Kalva, et al., “Fractional regularization to improve photoacoustic tomographic image reconstruction,” IEEE Trans. Med. Imaging 38(8), 1935–1947 (2019). [CrossRef]
47. B. E. Treeby and B. T. Cox, “Kwave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15(2), 021314 (2010). [CrossRef]
48. X. Zhu, Q. Huang, A. DiSpirito, et al., “Real-time whole-brain imaging of hemodynamics and oxygenation at micro-vessel resolution with ultrafast wide-field photoacoustic microscopy,” Light Sci. Appl. 11, 138 (2022). [CrossRef]
49. A. Hauptmann, F. Lucka, M. Betcke, et al., “Model-based learning for accelerated, limited-view 3-D photoacoustic tomography,” IEEE Trans. Med. Imaging 37(6), 1382–1393 (2018). [CrossRef]
50. B. E. Treeby, E. Z. Zhang, and B. T. Cox, “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Problems 26(11), 115003 (2010). [CrossRef]
51. W. F. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]
52. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]
