Open Access
Add to BibSonomy
Abstract
Currently, there are no non-invasive experimental methods available for measuring optical fluence distributions in tissue. We present photoacoustic tomography (PAT) as a method to approximate the relative optical fluence distribution in a homogeneous optically scattering medium. Three-dimensional photoacoustic images were captured with a near-full view PAT scanner in phantoms with known optical absorption and scatter properties. Resultant 3D PAT images were compared to the expected optical fluence distributions from Monte Carlo simulations and diffusion theory using volumetric and shape analysis. Volumetric analysis of PAT images compared well with the optical fluence distributions from simulation. Dice similarity coefficients ranged from 51 to 82%. The reduced scattering coefficient estimated from PAT images compared well to estimates from simulations for values below 0.5 mm-1. Near full-view PAT has been found to be useful for estimating the optical fluence distribution in an optically scattering medium. Further development is needed to extend the measurement range.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Optical techniques such as hyperspectral imaging, near-infrared spectroscopy (NIRS), diffuse optical tomography (DOT), and photoacoustic tomography (PAT) are becoming increasingly important as clinical diagnostic tools. The success of these methods depends on accurate estimates of the optical fluence distribution (i.e., light dosimetry) used to interrogate tissue at depths exceeding the ballistic regime [1]. Understanding photon propagation within optically absorbing and scattering media is crucial for designing and validating imaging systems, improving image reconstruction, and obtaining accurate results in clinical applications such as photodynamic therapy. The radiative transfer equation and its approximation, the diffusion equation, can be used to model light propagation [2]. In 1983, Wilson and Adam adapted the Monte Carlo method for light transport in tissue [3], which was then applied by Wang et al. [4]. However, non-invasive methods such as DOT that estimate optical fluence have limited utility in estimating the detailed fluence within scattering media due to its inherently low resolution. Invasive placement of a point detector inside a material to sample the fluence remains the only method to make direct optical fluence measurements [5,6]. Moreover, the invasive presence of a detector at the point of measurement influences the outcome.
In PAT, the medium absorbs laser light and converts optical energy to sound waves. Detection and reconstruction of the sound waves results in images that are proportional to the localized absorption coefficient (µa), optical fluence, and thermodynamic properties of the medium (Gruneisen coefficient). Estimation of the absorption coefficient has been the target of almost all PAT studies to date since it is directly related to chromophore concentration (e.g., hemoglobin) and reflective of the molecular constituents in the medium. Most studies have assumed uniform light distribution and homogeneous thermodynamic properties when reconstructing PAT images. For example, the 3D distribution of blood vessels in a female human breast was clearly visible even without prior knowledge of the variations in optical fluence [7]. Conversely, in a PAT experiment where there is an optically and thermodynamically homogenous medium, the localized absorption of light will be proportional to the localized optical fluence. Based on this assumption, we hypothesized that 3D PAT images will be surrogates of optical fluence when the medium has both uniform optical absorption and uniform thermodynamic properties. Previously, methods have been developed for fluence estimation during 2D photoacoustic imaging by Jeng et. al. (2021) and for photoacoustic microscopy by Zhy et. al. (2021) [8,9]. Here, we intend to estimate 3D macroscopic optical fluence distributions using PAT.
However, the use of PAT to estimate localized optical fluence is not straightforward. Detector coverage and detector bandwidth are critical for accurate image reconstruction. If coverage is too low or bandwidth too narrow, then acoustic signals emitted from photoacoustic sources go undetected and introduce negativity artifacts into the images [10]. Recently, L. Yip et al. (2022) developed a PAT system approaching near full-view with 3.8π steradian coverage and showed substantially increased performance compared to PAT systems with lesser coverage [11]. In this study, we take advantage of the near-full view PAT system to estimate the optical fluence distribution from a directional pulsed light source. We validated the measured fluence estimates against Monte Carlo simulations and diffusion theory.
2. Materials and methods
2.1. Phantom construction
Phantoms were fabricated from 2% w/v agarose (VWRVN605, CAS #9012-36-6, VWR International, Radnor, PA), India ink (Speedball 3378 Super Black India Ink, Speedball Art, Statesville, NC), and Intralipid (Intralipid-20%, Fresenius Kabi, Toronto, ON, CA). Each phantom was spherical in shape (54-mm diameter) and had two layers. The first layer consisted of agarose alone and had a maximum thickness of 17 mm. Each phantom had a unique combination of India Ink and Intralipid (IL) concentration in the second layer. In total, 17 phantoms were produced.
Phantoms were split into three sets (Table 1, A-C). Set A had a constant IL concentration of 0.5% with India ink concentrations ranging from 7 × 10−5 to 9 × 10−4 ml/ml (ml of India ink per ml of water). Set B had a constant India ink concentration of 1 × 10−4 ml/ml with IL concentrations ranging from 0.5% to 2.5%. Finally, set C had a constant India ink concentration of 1 × 10−4 ml/ml with IL concentrations ranging from 2.5% to 4.5%. Table 1 also reports the resulting optical properties of the phantoms. Spectrophotometric measurements described in the supplementary data were used to estimate absorption and scattering coefficients. The reduced scattering coefficient was estimated from the anisotropy coefficient reported in literature in combination with our estimates of the scattering coefficient [12].
Table 1. Characteristics of phantoms. Absorption coefficients (µa) and scattering coefficients (µs) were measured using a spectrophotometer, while anisotropy (g) was taken from literature. The reduced scattering coefficients (µ’s) were calculated from the other variables. II: India ink, IL: Intralipid
Each phantom was constructed using a two-piece mold (see supplementary data), which was 3D printed out of photopolymer resin (Clear Resin, Form 2, Formlabs, Somerville, MA). A 3D printed mount was integrated into each phantom. For each phantom, agarose powder was mixed with deionized water and heated in a microwave oven until fully dissolved and bubbling. The liquid was weighed before and after heating to account for evaporative losses. The agarose solution was allowed to cool to 56°C in a water bath before India ink was added. A sample was then preserved for spectrophotometric analysis, followed by the addition of IL. The solution was then poured through the top hole of the mold until the liquid reached the level-setting holes part way up the shell. After cooling, a clear agarose layer was poured on top to fill the mold.
2.2. PAT imaging system
A near-full view 3D PAT system was used for all PAT experiments [11]. Briefly, a 64-channel transducer array was composed of 16 identical transducer modules mounted to an aluminum ring (Ø ≈ 280 mm) (Fig. 1(a)). Each transducer module contained 4 transducer elements. Two of the four elements had a center frequency of 0.4 MHz with a bandwidth of 160% (one-way bandwidth at -6 dB). The remaining two elements had a center frequency of 0.9 MHz and a bandwidth of 185% (one-way bandwidth at -6 dB). Since the imaging system was designed for bulk tissue imaging, the center frequencies of the transducers were not strictly optimized for fluence distribution imaging [11]. The transducer array was mounted to a 6-axis robot (C3, Epson America Inc., Los Alamitos, CA) (Fig. 1(b)). A tunable Nd:YAG laser (Phocus, Opotek Inc., Carlsbad, CA, 680-950 nm, 10-Hz pulse repetition rate, 5-ns pulse width) was used for illumination. The transducer array was rotated to achieve 3.8π steradian coverage and translated to increase the effective field of view.
Fig. 1. PAT system design illustrating the detection and illumination setups. (a) CAD drawing of the circular ring array with 16 detector modules (two transducer types, A and B) surrounding the spherical absorber, mounting frame, and fiber optic cable. (b) CAD drawing of the array mounted to the 6-axis robot with arrows in green denoting azimuthal angle phi ($\mathit{\Phi}$) and elevational angle theta ($\theta$) rotational directions. (c) CAD drawing close-up of the spherical absorber, mounting frame, and optical fiber. Note the axes.
2.3. Experimental and simulated optical fluence distribution
Phantoms were imaged with the PAT system described above. Two scans were performed for each phantom, one using directional illumination and the other, diffuse. For directional illumination, the laser was coupled to a high-power 1 to 2 fiber optic cable (Excelitas Canada Inc., Mississauga, ON, CA, NA 0.17 in water, fused input). One output of the fiber optic cable was directed at the phantom and the second output was terminated at a beam dump (Fig. 1(c)). Imaging was performed with 800-nm laser light (to maximize energy output). For diffuse illumination, the two-legged fiber optic cable was replaced with a high-power 1 to 8 fiber optic cable (Excelitas Canada Inc., Mississauga, ON, CA, NA 0.17 in water, fused input). The eight fiber outputs were evenly spaced along the perimeter of the transducer ring and aimed toward the geometric center of the ring. Imaging was performed with 680-nm laser light (to enhance surface absorption). The photodiode (DET10A, Thorlabs Inc., Newton, NJ) signal was used to correct the transducer signals for variations in shot-to-shot laser fluence. Photodiode-corrected PA signals were reconstructed using the universal back projection algorithm [13] with directivity weighting [14]. Reconstructions were performed in MATLAB using the PATLAB toolkit developed in our laboratory [15].
Monte Carlo simulations were performed with TracePro software (Lambda Research Corporation, Littleton, MA). For each simulation, 2 × 106 rays were used. A CAD model of the holder and fiber optic cable, as well as the two-layer phantom were imported into the software package. The fiber optic cable had a NA of 0.17 (in water) and an output diameter of 3 mm. Agarose was assumed to have the same optical properties as water, which were taken from Buiteveld et al. (1994) [16]. The optical properties of the phantoms were set to those in Table 1. The surrounding medium was modelled with the optical properties of water.
Photoacoustic images collected with directional and diffuse illumination were inherently co-registered as the location of the phantom was not changed between scans. The 3D PA volumes were imported in 3D Slicer (v.4.13.0) [17]. PAT images obtained with the directional illumination were cropped to remove signal from the non-optically absorbing agarose section. For visualization purposes, threshold and color scale were adjusted. Simulation results were exported as 3D data sets of irradiance and imported into 3D Slicer for processing and visualization. Simulation results were co-registered to the experimental data using the phantom holder and curved edges of the phantom as landmarks. An intensity threshold was applied to the simulation results.
2.4. Volumetric and statistical comparison between experimental and simulated optical fluence distribution
For each of the three sets of phantoms, the experimental and simulated optical fluence volumes were segmented in 3D Slicer and compared. Dice similarity coefficient (DSC) was then used to compare between each segmented pair of experimental and simulated data. The DSC measured the similarity of segmented volumes, where 0 represented no overlap and 1 represented a perfect match [18].
2.5. Depth-resolved fluence estimation and reduced scattering coefficient extraction
We computed the reduced scattering coefficient (µ′s) from the acquired data by modelling against the depth-resolve fluence model of optically turbid media (e.g., diffusion approximation). The model is applicable when µ′s ≫ µa. Typically µ′s/µa must be > 10, but in practice it has been shown that µ′s/µa > 5 or even µ′s /µa > 3 provides accurate results [19,20]. Moreover, the depth needs to exceed the transport mean free path (TMFP). Phantoms 9-17 (phantom sets B and C in Table 1) were included in the results as they satisfied this requirement. Intensity profiles along the axis of the directional light beam (i.e., depth-resolved fluence profiles along the X-axis (see Fig. 1)) were analyzed to compare the optical property estimates from PAT, Monte Carlo, and diffusion theory. Intensity values deeper than the TMFP at each cross-sectional plane were summed over the length of the image. For PAT images, intensity values < 8% were discarded to reduce the effect of noise on the processed results. The µ′s values for each phantom were extracted both from the PAT and Monte Carlo data using Eq. (1), where $\delta$ was the decay constant of the exponentially decaying fluence profile and µa was estimated from spectroscopic measurements [2].
3. Results
3.1. Qualitative comparison of experimental and simulated optical fluence distribution
Two-dimensional PAT images extracted from 3D PAT data are presented in Fig. 2(a)-(c), where PAT images from directional illumination (red color scale) are superimposed on PAT images from diffuse illumination (grayscale). In comparison, the simulated optical fluence distribution superimposed upon a CAD model of the phantom is presented in Fig. 2(d)-(f). Additionally, examples of raw photoacoustic signals with indications of FWHM are shown in Fig. 2(g)-(i). In Fig. 2, the first column (a, d, g) represents Phantom 7, the second (b, e, h) represents Phantom 3, and the third (c, f, i) represents Phantom 12 (see Table 1).
Fig. 2. Experimental PAT reconstructions compared to optical fluence profiles estimated with Monte Carlo simulations and examples of raw photoacoustic signals. Experimentally obtained 2D PAT images extracted from the 3D PAT images are shown in panels (a-c), where the image from directional illumination (red color scale) is superimposed upon the image obtained with diffuse illumination (grayscale). Simulated optical fluence profiles (d-f) are shown superimposed upon a profile of the CAD model. The directional light source in each case was directed in negative x-direction. Raw photoacoustic signal examples from a detector located approximately perpendicular to the xz plane in the middle of the directional illumination reconstructions (indicated by a star symbol in a-c) are shown in (g-i). Each signal measurement panel includes an indication of the measured FWHM. Data in panels (a), (d) and (g) were obtained with Phantom 7. Data in panels (b), (e) and (h) were obtained with Phantom 3. Data in panels (c), (f) and (i) were obtained with Phantom 12. Color, scale bars and axes are common to all images. The grayscale applies to each PAT image individually (i.e., images cannot be compared with each other).
When comparing the PAT images across phantoms (Fig. 2(a)-(c)), differences in size of the optical fluence profile were noted. Phantom 3 had the lowest total attenuation (including the contributions from scatter and absorption) and the longest fluence profile along the x-axis (Fig. 2(b)). Phantom 7 had the second lowest total attenuation and the second longest fluence profile in the x-direction (Fig. 2(a)). Phantom 12 had the highest total attenuation and the shortest fluence profile in the x-direction (Fig. 2(c)). The scattering coefficient was the same for Phantoms 3 (Fig. 2(b)) and 7 (Fig. 2(a)), but Phantom 7 had a higher absorption coefficient compared to Phantom 3. The width of the optical fluence profile along the z-axis appeared to be similar for Phantoms 3 and 7. The absorption coefficient was identical for phantoms 3 (Fig. 2(b)) and 12 (Fig. 2(c)), but Phantom 12 had a higher scattering coefficient and a wider fluence profile in the z-direction.
The simulated results were qualitatively similar to the experimental findings for each phantom in terms of the optical fluence profile length and width. Some of the noticeable differences with the selected threshold include shorter optical fluence profiles along the x-axis and wider profiles along the z-axis (Fig. 2(d)-(f)) compared to the experimental results (Fig. 2(a)-(c)).
The raw photoacoustic signals shown in Fig. 2(g)-(i) were measured along the y-axis of the acquired 3D volume (perpendicular to the reader’s perspective of Fig. 2(a)-(c)). Phantom 7 had the highest absorption and the highest signal intensity (Fig. 2 g), while the signal from Phantom 12 with the highest attenuation and scattering had a noticeably low signal-to-noise ratio (Fig. 2(i)). Assuming the fluence distribution was symmetrical about the x-axis, the measured FWHM of the raw photoacoustic signal was consistent with the observed width of the experimental fluence in the z-direction for phantoms 3 and 12 (Fig. 2(b) and (c) visually compared to FWHM measurements Fig. 2(h)-(i)).
3.2. Volumetric and statistical comparison between experimental and simulated optical fluence distribution
Photoacoustic images and simulated optical fluence were segmented using simple threshold processing. The resultant segmentations were compared using volume and shape metrics (Fig. 3). Phantom set A had constant µs but varied in µa. The segmented simulated volumes for phantom set A decreased exponentially as µa increased (R2 = 0.98, Fig. 3(a)). The segmented experimental volumes also decreased as µa increased until reaching µa values larger than ∼0.2 mm-1 (R2 = 0.53, Fig. 3(a)). The decay rate of the exponential fits for both simulated and experimental data was similar and was only offset by 71 mm-3. Phantom sets B and C had constant µa but varied in µs. In a similar manner to phantom set A, the simulated segmented volumes followed an exponential decay as µs increased in value (R2 = 0.99, Fig. 3(b)). The segmented photoacoustic volumes followed a decreasing trend with increasing µs, had good overall numerical correspondence and a similar decay rate (deviated by 0.1), however the exponential fit was not as good (R2 = 0.48, Fig. 3(b)).
Fig. 3. Volumetric and shape similarity comparison between segmented volumes obtained from simulations of optical fluence and PAT. (a) Dependence of volume estimated from segmented PAT image data on absorption coefficient (Phantom Set A). (b) Similar to (a), except in relation to scattering coefficient (Phantom Sets B and C). (c & d) Dice similarity coefficient for comparison between simulated and experimental data for Phantom Sets A and B&C, respectively.
Shape similarity of the segmentations was compared using the Dice similarity coefficient (DSC). In the case of phantom set A, values of DSC were in the range of 0.51 to 0.82 for all µa (Fig. 3(c)), with an average of 0.73 ± 0.10 (SD). The highest values of DSC were observed for the phantom with µa ≈ 0.03 mm-1. Values of DSC appeared to decrease with µa before and after the peak. In the case of phantom sets B and C, values of DSC ranged from 0.66 to 0.84 for all µs (Fig. 3(d)), with an average of 0.79 ± 0.05 (SD). No apparent dependance of DSC on µs was observed.
3.3. Depth- resolved fluence estimation and reduced scattering coefficient extraction
The depth-resolved experimental PAT intensity profile for Phantom 12 was compared to depth-resolved optical fluence profiles obtained from Monte Carlo simulations and diffusion theory (Eq. (1)) (Fig. 4(a)). In all cases, the profiles decreased monotonically with depth into the medium. At depths beyond a few millimeters (i.e., beyond one TMFP), each profile had a linear response on a log-linear plot with almost identical exponential decay constants and R2 values approaching 1 (Fig. 4(a)).
Fig. 4. Estimation of the depth-resolved PAT signal intensity and reduced scattering coefficient compared to the depth-resolved optical fluence, and reduced scattering coefficient estimated from Monte Carlo simulations and diffusion theory. (a) Optical fluence profile (experimental (PAT), simulated and theoretical using Eq. (1) and Table 1 optical properties) for phantom 12. (b) Reduced scattering coefficient recovered from experimental and simulated depth-resolved profiles for phantoms 9-17 and compared to expected values (Table 1).
Figure 4(b) shows the µ’s values estimated from experimental PAT images of phantoms 9-17 (Table 1), Monte Carlo simulations, and estimates from diffusion theory. The measured values from the PAT intensity profiles were compared to the expected (Table 1) and simulated values of µ′s. The µ′s estimates derived from the simulations compared well to the expected values. Simulated results had a slight slope deviation from the expected values (Fig. 4(b)). Estimates of µ′s from the PAT images matched well with theoretical µ′s estimates up to µ′s ≈ 0.5 mm-1 (R2= 0.96), but did not match well for larger values of µ′s.
4. Discussion
In this study, we demonstrated that near-full view PAT can be used to estimate optical fluence within a homogeneous tissue-mimicking phantom. Photoacoustic images of each phantom were obtained with a near-full view PAT system equipped with directional and diffuse pulsed light sources. The images captured with directional illumination were then compared to Monte Carlo simulations of the expected optical fluence distributions. Qualitative, volumetric, and shape similarity comparisons demonstrated that near-full view PAT provided an accurate representation of the expected optical fluence distribution computed from Monte Carlo simulation for a range of phantoms each with a unique combination of µa and µs. In addition, analysis of near-full view PAT images resulted in accurate estimations of µ′s for values of µ′s ≤ 0.5 mm−1. Estimates deviated from expected values for µ′s > 0.5 mm−1 (see Fig. 4(b)). We hypothesized that the reduced measurement range of the PAT method was related to the sensitivity and bandwidth of the PAT system. The transducers had low sensitivity at high frequencies and were relatively far from the phantom (≈ 14 cm). High-frequency emissions due to confined distributions of light at higher scattering may have been undetectable in the PAT images.
The segmented image volumes from PAT followed a similar dependance on attenuation coefficient compared to the segmented fluence distributions from Monte Carlo simulations (Fig. 3). Near-full view PAT images were dependent on the changes in the optical properties between phantoms. We also found that the segmented PAT images matched with the TracePro simulated optical fluence maps with an average of 76% (Fig. 3(c) and (d)). We observed that 3D PAT intensity distributions more closely matched the optical fluence distributions from simulations in some image planes compared to others (see supplementary data, profile analysis). This was likely due to the design of the 3D PAT system, which was missing ≈ 0.2$\pi$ of coverage, concentrated along a gap in one side [11]. It was expected that the gap in coverage would lead to artifacts in the direction where view-angle was limited. In addition, the 3D-printed phantom holder was located centrally with respect to the coverage gap and likely interfered with the normal propagation of PA signals. These artifacts can be seen in Fig. 5 and were major contributors to the mismatch in shape between simulated and experimental data. Additionally, when looking at Fig. 3, it can be seen that at higher absorption values, the experimentally determined volumes plateau, while simulated volumes exponentially decrease. This result is likely another effect of limited detector bandwidth, which resulted in our inability to accurately resolve the dimensions of smaller objects.
Fig. 5. Example PAT image results with thresholded segmentations (green) demonstrating artifacts (see green arrow) due to limited coverage and acoustic interference of the phantom holder. The red outline shows the boundaries of the agarose phantom. (a) XZ slice without artifacts. (b) XY slice with artifacts.
The study had several limitations. First, we did not account for acoustic effects when comparing the PAT images to the simulated optical fluence distributions. For example, the study did not factor in acoustic effects such as attenuation of agarose and water and potential reflections from the holder and fiber bundle. These acoustic effects may have caused our PAT results to differ from the simulated results. Future work in this area could potentially use the simulated optical fluence maps as initial pressure profiles for k-Wave [21], which could then be used to simulate PA signals and, in turn, reconstructed into a PAT image for comparison to the experimentally obtained PAT image. Second, for phantoms comprised of India ink and Intralipid, both µa and µs were measured experimentally (see supplementary data). However, we were unable to measure the anisotropy value (g) of Intralipid, which is normally recommended due to batch-to-batch variability [19]. This could have been a source of error leading to the observed small discrepancies between the PAT and Monte Carlo results.
In conclusion, we have demonstrated a method to experimentally and non-invasively measure fluence distribution within a tissue-mimicking phantom that approximates the expected fluence distribution based on diffusion theory and Monte Carlo simulations. While there are many areas for potential improvement, we propose that this method may be extended to imaging phantoms and clinical specimens with both homogeneous and heterogeneous optical properties, and lead to a better understanding of fluence distribution due to photon propagation in optically scattering media.
Funding
Canadian Institutes of Health Research (356794); Natural Sciences and Engineering Research Council (RGPIN-2019-06914).
Acknowledgments
PO was supported through a studentship award provided by the Lawson Health Research Institute, while LCMY was supported through studentships funded by the Breast Cancer Society of Canada and an Ontario Graduate Scholarship. We would also like to thank Drs. Androu Abdalmalak and Daniel Milej for their valuable insight. This study was conducted with the support of FACIT through funding provided by the Government of Ontario. This work was supported by research funding from Multimagnetics Inc.
Disclosures
Jeffrey Carson reports a relationship with Multimagnetics Inc. that includes equity or stocks and funding grants. Jeffrey Carson has patent #US9128032 assigned to Multimagnetics Inc.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. B. C. Wilson, “Modelling and measurements of light propagation in tissue for diagnostic and therapeutic applications,”. Laser Systems for Photobiology and Photomedicine, vol 252 (Springer, 1991).
2. L.V. Wang and H. I. Wu, Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).
3. B. C. Wilson and A. G. Adam, “A Monte Carlo model for the absorption and flux distributions of light in tissue,” Med. Phys. 10(6), 824–830 (1983). [CrossRef]
4. L. V. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Methods Programs Biomed. 47(2), 131–146 (1995). [CrossRef]
5. T. C. Zhu, J. C. Finlay, and S. M. Hahn, “Determination of the distribution of light, optical properties, drug concentration, and tissue oxygenation in-vivo in human prostate during motexafin lutetium-mediated photodynamic therapy,” J. Photochem. Photobiol., B 79(3), 231–241 (2005). [CrossRef]
6. W. M. Star, “Light dosimetry in vivo,” Phys. Med. Biol. 42(5), 763–787 (1997). [CrossRef]
7. R. A. Kruger, R. B. Lam, D. R. Reinecke, S. P. Del Rio, and R. P. Doyle, “Photoacoustic angiography of the breast,” Med. Phys. 37(11), 6096–6100 (2010). [CrossRef]
8. G. S. Jeng, M. L. Li, M. Kim, S. J. Yoon, J. J. Pitre Jr, D. S. Li, I. Pelivanon, and M. O’Donnell, “Real-time interleaved spectroscopic photoacoustic and ultrasound (PAUS) scanning with simultaneous fluence compensation and motion correction,” Nat. Commun. 12(1), 716 (2021). [CrossRef]
9. J. Zhu, C. Liu, Y. Liu, J. Chen, Y. Zhang, K. Yao, and L. Wang, “Self-Fluence-Compensated Functional Photoacoustic Microscopy,” IEEE Trans. Med. Imaging 40(12), 3856–3866 (2021). [CrossRef]
10. K. Shen, S. Liu, T. Feng, J. Yuan, B. Zhu, and C. Tian, “Negativity artifacts in back-projection based photoacoustic tomography,” J. Phys. D: Appl. Phys. 54(7), 074001 (2021). [CrossRef]
11. L. C. M. Yip, P. Omidi, E. Rascevska, and J. J. L. Carson, “Approaching closed spherical, full-view detection for photoacoustic tomography,” J. Biomed. Opt. 27(08), 1–19 (2022). [CrossRef]
12. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907–5925 (2008). [CrossRef]
13. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E 71(1), 016706 (2005). [CrossRef]
14. B. T. Cox and B. E. Treeby, “Effect of sensor directionality on photoacoustic imaging: a study using the k-wave toolbox,” SPIE Proceed. 7564, 75640I (2010). [CrossRef]
15. P. Omidi, L. C. M. Yip, E. Rascevska, M. Diop, and J. J. L. Carson, “PATLAB: A graphical computational software package for photoacoustic computed tomography research,” Photoacoustics 28, 100404 (2022). [CrossRef]
16. H. Buiteveld, J. H. M. Hakvoort, and M. Donze, “Optical properties of pure water,” Ocean. Opt. XII 2258 174–183 (1994).
17. A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J. C. Fillion-Robin, S. Pujol, C. Bauer, D. Jennings, F. Fennessy, M. Sonka, J. Buatti, S. Aylward, J. V. Miller, S. Pieper, and R. Kikinis, “3D Slicer as an image computing platform for the Quantitative Imaging Network,” Magn. Reson. Imaging 30(9), 1323–1341 (2012). [CrossRef]
18. A. A. Taha and A. Hanbury, “Metrics for evaluating 3D medical image segmentation: Analysis, selection, and tool,” BMC Med. Imaging 15(1), 29 (2015). [CrossRef]
19. S. L. Jacques and S. A. Prahl, “Diffusion Theory,” https://omlc.org/classroom/ece532/class5/limits.html, (1998).
20. H. Xu, T. J. Farrell, and M. S. Patterson, “Investigation of light propagation models to determine the optical properties of tissue from interstitial frequency domain fluence measurements,” J. Biomed. Opt. 11(4), 041104 (2006). [CrossRef]
21. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields,” J. Biomed. Opt. 15(2), 021314 (2010). [CrossRef]
