{'fraction_non_alphanumeric': 0.032594589198362785, 'fraction_numerical': 0.011954677049016671, 'mean_word_length': 4.685965659145736, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Correction of wavelength-dependent laser fluence in swept-beam spectroscopic photoacoustic imaging with a hand-held probe

Recently, we demonstrated an integrated photoacoustic (PA) and ultrasound (PAUS) system using a kHz-rate wavelength-tunable laser and a swept-beam delivery approach. It irradiates a medium using a narrow laser beam swept at high repetition rate (∼1 kHz) over the desired imaging area, in contrast to the conventional PA approach using broad-beam illumination at a low repetition rate (10−50 Hz). Here, we present a method to correct the wavelength-dependent fluence distribution and demonstrate its performance in phantom studies using a conventional limited view/bandwidth hand-held US probe. We adopted analytic fluence models, extending diffusion theory for the case of a pencil beam obliquely incident on an optically homogenous turbid medium, and developed a robust method to estimate fluence attenuation in the medium using PA measurements acquired from multiple fiber-irradiation positions swept at a kHz rate. We conducted comprehensive simulation tests and phantom studies using well-known contrast-agents to validate the reliability of the fluence model and its spectral corrections.


Introduction
Quantitative spectroscopic photoacoustic (PA) imaging relies on the known optical absorption spectrum of target chromophores to estimate their concentrations. The PA signal, however, depends not only on the optical absorption spectrum of a target but also on the wavelengthdependent optical fluence at the same site. A critical challenge for in vivo, PA spectroscopic imaging is predicting light transport noninvasively in a turbid medium [1,2]. The spatial distribution of fluence depends on the optical absorption and scattering properties of that medium, but they are usually unknown in advance. In addition, since light attenuation in a medium is a function of wavelength, the medium induces spectral distortion (sometimes called 'coloring') between the nominal PA-measured optical absorption spectrum and the true target's spectrum [1]. Therefore, for true quantitative spectroscopic (or multiwavelength) PA imaging that can identify molecular constituents in the medium under study, the optical fluence must be simultaneously estimated and compensated. This requires not only an appropriate fluence model, but also an effective method to estimate the fluence distribution at each wavelength from PA measurements.
Fluence correction methods have been proposed for homogenous absorption and scattering media such that the light distribution and PA measurements can be closely represented by simple formulas or mathematical expressions [3][4][5][6][7]. Unfortunately, very few of them translate into clinical tools because most explicitly depend on prior knowledge of the medium's optical properties. Indeed, optical constants reported in the literature may vary a few orders of magnitude depending on the measurement technique, tissue condition, and geometry [8]. Background tissue properties may also change dynamically based on tissue blood content and oxygenation level. Optical constants measured even with the same device, the same experimental conditions, and for the same person may vary. Furthermore, these properties may change during medical procedures and interventions.
Many recent studies have shown the potential of deep learning schemes to access more complex tissue structures and heterogeneous optical properties within a short computation time [13][14][15][16][17][18][19]. Their primary challenge is the generation of realistic training data/images. Despite success in simulations, they have not been shown yet to work broadly for in-vitro or in-vivo conditions. Unfortunately, full view and bandwidth systems used for small animal PA imaging cannot be easily adapted to most clinical applications in humans. Due to limited light penetration, PA systems intended for most clinical applications must use hand-held, limited view array transducers [21], leading to strong artifacts in structural PA images [22]. Diagnostic ultrasound (US) co-registered with PA (e.g. PAUS methods) has many advantages for clinical translation of PA technology, but the limited view/bandwidth of US arrays degrades PA reconstruction. Even with these constraints, PA imaging, particularly spectroscopic imaging, brings a molecular dimension to US. However, there have been no successful demonstrations to date of wavelengthdependent laser fluence correction for these systems that is appropriate for in vivo, let alone ex vivo, applications.
This study focuses on wavelength-dependent laser fluence correction for PA imaging using a limited view/bandwidth diagnostic US probe [21] appropriate for real-time monitoring of interventional procedures. The main target is identifying small and sparse absorbers such as labelled drugs and contrast agents or quantifying functional changes in microvessels in quasi-homogenous optically turbid tissue. In this geometry, fluence correction should be constrained and physics-based because limited view and narrow bandwidth lose some low-frequency components during reconstruction (beamforming) of light-absorbing structures [23]. Note that laser fluence estimation without the help of other tools and at near real-time rates is highly desirable.
In Ref. [24], point source tissue illumination from different positions along the tissue surface was proposed to calculate the laser fluence distribution within the medium and then apply the estimated fluence for PA spectral decomposition. However, the authors did not show how to integrate this method into real-time PA scanners since broad beam illumination is most commonly used in PA imaging systems [24]. Leveraging the results of this study, we recently introduced a real-time interleaved photoacoustic-ultrasound (PAUS) fast-sweep scanner [25], where unlike previous delivery systems coupling laser pulses into all fibers in a bundle simultaneously, light is coupled into individual fibers sequentially (see Fig. 1), but at a very high rate.
A unique diode-pumped wavelength tunable (700 nm -900 nm) laser emitting about 1 mJ pulses at 1000 Hz, with wavelength switching in less than 1 ms for any arbitrary wavelength order, was designed especially for the fast-sweep PAUS scanner. To maximize exposure, we illuminate with a narrow (∼ 1 mm in diameter) beam and switch it from fiber-to-fiber at 1000 Hz, resulting in one loop around the US probe per single-wavelength frame in only 20 ms (Fig. 1). The next loop uses another wavelength without delay; the procedure repeats over all wavelengths. That is, instead of illuminating with a broad beam, we use fast-scanning (or fast-sweep) over the same illumination area. In our PAUS system, ten fibers are uniformly spaced along each elevational edge of the US array (e.g. 20 fibers in total, Figs. 1,2). Every laser shot from a single fiber is followed by a sub-image reconstruction, i.e. 20 sub-images in all, which are then coherently summed to form the full PA frame. The kHz rate enables 50 full 'loops' of the laser beam around the probe per second, resulting in a 50 Hz PA frame rate. For stable spectral decomposition, 10 wavelengths (i.e., 700, 715−875 nm every 20 nm) form the spectroscopic sequence. More details on the system can be found in Ref [25].
Here we describe in detail how to use partial PA images from every fiber to estimate laser fluence. Indeed, when light emerges from different fibers, it propagates different distances to a target and, therefore, the amplitude of a partial PA image reconstructed from a single fiber illumination will depend on that fiber's position around the US probe (see Fig. 2). For every fiber light source, we adopted analytical fluence models extending diffusion theory for narrow beam (pencil beam) illumination obliquely incident on a semi-infinite homogeneous scattering medium [26]. Using this model, we explore methods to extract optical parameters to assess the fluence distribution within the medium using data acquired at different fiber positions and wavelengths. In particular, estimated parameters include the effective light attenuation coefficient µ eff , and reduced light scattering µ s coefficients of a turbid medium.
Building on the basic analytic model, we then explore robust methods to estimate these parameters from noisy data. Simulations were conducted to clarify quantitative errors arising from adopted optical fluence models, variations in medium properties, and measurement noise levels. Finally, we demonstrate the reliability of the correction methods via phantom studies.

Photoacoustic signal
Acoustic pressure obtained for the k th optical fiber, optical wavelength j at discrete position r i can be denoted as where … k {1, 2, , 20} and … j {0,1, , 9}, is the Grü nesien parameter, is the optical absorption coefficient of a target and is the light fluence. We assume here that the parameter r ( ) governing PA efficiency is constant over space = x y z r r r { | ( , , ), } and wavelength. This assumption is not critical for relative spectroscopic measurements, but is important for absolute concentration estimates of specific Real-time integrated photoacoustic and ultrasound (PAUS) system used in these studies. The ultrasound system programmably controls the laser, motor controller, and US transducer. The motor controller synchronizes emission with the centers of 20 fibers in the bundle, delivering a trigger to the US system when properly aligned for each fiber. The US system then externally triggers the compact laser, transmitting a pulse at about a 1 kHz repetition rate with a wavelength switchable from pulse-to-pulse over the range from 700 nm to 900 nm. With absolute position control, a precise rate is not needed for external laser triggering, ensuring maximal light delivery to each fiber. Motor speed variations only slightly alter the overall frame rate of 50 Hz. A total of 20 fibers are arranged on two sides of the linear array US transducer, as shown in the zoomed front view in the bottom right corner. The system mixes laser beam and focused US transmissions for interleaved PA spectroscopic imaging and US B-mode imaging. chromophores producing the PA signal. The fluence distribution depends on the optical properties of the medium. Absorption arising from imaged targets μ a also disturbs the fluence, but the effect is insignificant if the target is small [27]. We thus assume that the acoustic pressure (PA signal) is linearly proportional to the absorption coefficient of the target.
The absorption spectrum is given as where L is the number of chromophore (absorber) types and C l and l j ( ) are the concentration and absorption coefficient of the l th chromophore type, respectively. The ultimate goal of quantitative spectroscopic PA imaging is to estimate the relative, or absolute, concentration of a particular chromophore of interest at each position within the imaging field from measurements given the known absorption spectrum l j ( ) of that chromophore. Thus, correct modeling and computation of the wavelength-dependent optical fluence distribution r ( ) at every source position is required to guarantee accurate quantification of the target chromophore.

Optical fluence model
Light illumination in a medium can be simply expressed using the diffusion approximation, where Q denotes the source, µ a denotes the absorption coefficient within a medium, and denotes the fluence [26]. The parameter D is the diffusion coefficient given as where µ s is the scattering coefficient of the medium, µ s is the reduced scattering coefficient and g is the anisotropy factor. This expression is derived from the radiative transfer equation (RTE) for high scattering (µ µ ) s a and a nearly isotropic medium [26]. We assume that biological tissue can be regarded as 'macro-homogeneous'. If the scale of optical micro-heterogeneities is smaller than a few photon transport mean free paths, light fluence will be smooth and can be characterized by a macro-homogeneous medium model [28]. For a homogeneous a a s s and a point source , the solution to the diffusion equation using the Green's function is given as [29][30][31] is the effective attenuation coefficient and r is the spatial position of the source. The solution can be modified for our acquisition environment in which a pencil beam is obliquely incident on a semi-infinite scattering medium. As illustrated in Fig. 2, the beam can be converted into two isotropic point sources mirror-symmetric about the extrapolated boundary = z z b due to the refractive-index mismatch between media [26]. The source below the boundary helps describe the transition from ballistic to scattering regimes, and its symmetrical source is added to satisfy a boundary condition described in Appendix A1. Note that every fiber is mounted on the transducer and covered by a thick acrylic holder and protection glass (BK-7 optical glass). Thus, the main contributor to the mismatch with the scattering medium is the solid transducer rather than the ambient medium (air). Appendix A1 describes the detailed computation of z b using the refractive indexes n of the media.
The relative positions of the point sources with respect to the boundary are determined by the transport mean free path = l t µ 1 s [26].
For the k th fiber whose tip position is = x y z r ( , , ) k k k k and incident angle , the fluence solution can be expressed as . Note that the unknown parameters are 1 , µ s ' and µ eff . Other parameters can be derived from them or known initially.
If the positions of two imaginary isotropic sources in Fig.2 are close to each other due to a very small l t (e.g. very high µ s ), the fluence can be simplified to an asymptotic expression where 2 is a scalar. The derivation is presented in Appendix (B). Here, unknown parameters are only 2 and µ eff , i.e. µ eff is the only parameter defining the relative distribution of laser fluence in the medium. Held et al. employed this model [22].
In this paper, analytic fluence expressions in Eqs. (6) and (7) are called Model I and II, respectively. Thus, when Model II is quite accurate, it makes fluence estimation from experimental data very simple and stable. We will explore below how well both models compare to Monte-Carlo simulations, and under what conditions these models can be used for laser fluence assessment from swept-beam PA measurements.

Optical fluence estimation
PA measurement can be expressed as where j and k are wavelength and optical fiber indexes, respectively. n k is the system noise containing a DC component b k and zero-mean Gaussian noise n k . Control data were recorded by assigning zero laser power at the first wavelength ( j = 0) to estimate the noise bias as , where | | is the cardinality of the entire domain . The unbiased measurement can then be obtained as Note that the unknown parameters needed to estimate p r ( ) provided Model I is employed. To enhance estimation efficiency, we selected position indexes as where denotes the threshold value. Also, we normalized measurements to reduce the degrees of freedom as where i , and y r ( ) ranges from 0 to 1. The optimal parameters for the j th wavelength can be estimated as is the weight such that a position with higher SNR contributes more to the estimate. Likewise, Model II can be used for estimation, but the search parameter is only μ eff j ( ) for every wave- can be obtained by substituting optical parameter estimates into either Model I or II.

Fluence correction of light absorption spectrum
Assume that one type of chromophore is located at a local position of interest and its absorption spectrum is known as j . The optical absorption spectrum obtained from PA measurements = d y j k j k , is distorted due to the wavelength-dependent fluence of the surrounding turbid medium (biological tissue, for example), as shown in Eq. (1). The spectrum can be corrected using the fluence estimateˆj k , obtained by solving the least-squares problem, , for every wavelength j. By comparing normalized d j and c j with a j , fluence estimation accuracy can be computed.

Simulation parameters
The primary purpose of simulations is to verify the proposed fluence models and associated estimation methods described in Sections 2.2 and 2.3. First, Models I and II were compared to ground truth Monte Carlo simulations. We adopted the medium geometry of Fig. 2 for all simulations. All media were approximately represented as cubical shapes. Their sizes, optical parameters, and refractive indexes are summarized in Table 1. We assumed that one of the fibers transmits light into the scattering medium at the interface. The location and tilt angle of the fiber tip are (0 mm, 5.68 mm, 0 mm) and 35 , respectively. The total number of photon packets for Monte Carlo tests is 20 million. The simulation used the open source program, MCX Studio [32].
We also compared errors in parameter estimation for Model I compared to Model II at a fixed optical wavelength. We generated synthetic data using Model I (Eq. (6)) assuming that Model I is very close to ground truth. The noise power is set to zero to find the contribution of model difference to the error. We assumed a point target located at position = x z r { ( mm, 0 mm, mm)} in a scattering medium whose optical parameters are µ a , µ s and µ eff . For every (z, µ s ), we conducted 100 simulation tests to estimate µ eff by randomly changing values of other parameters in the range. A fractional error sample for every test is given as = × µ µ µ (ˆ)/ 100 eff eff eff (%), where μ eff denotes the estimate from Model I. We averaged 100 error samples for every (z, µ s ) to monitor the mean error pattern z ( , µ ) s over penetration depth and scattering coefficient.
In the last simulation, we investigated the accuracy of parameter estimation for different wavelengths in a turbid medium where µ a and µ s are wavelength dependent. We generated synthetic data using Model I as in the previous simulation, and added white Gaussian noise to mimic experimental conditions and explore reconstruction algorithm stability to noise. The total number of SNR levels was 15, ranging from 20 dB to 50 dB, and the number of test data sets for each SNR was 100. SNR is defined as where n 2 is the noise variance. For every test, we assumed SNR is constant over wavelength. A point target is located at position providing the most complicated conditions for light absorption spectrum reconstruction using PA signals. We used this scattering function for the test medium in the third simulation but assumed the absorption coefficient = µ ( 0.03 cm ) a j ( ) 1 is constant for all wavelengths. This does not reduce the generality of the results because (as we show below) reconstruction error is not highly sensitive to variations in light absorption.
Optical parameters were estimated from noise contaminated data . We examined the mean correlation coefficient using 100 samples for every SNR level. Fig. 3(a) shows slices of a 3D fluence distribution x y z ( , , ) from the Monte Carlo simulation. We used the distribution of optical fluence z (0,0, ) along the z axis as a reference to compare with corresponding distributions calculated using Model I and Model II (given by Eqs. (6) and (7) respectively). Fig. 3(c-e) show the distributions obtained for different µ s (2 cm 1 , 5 cm 1 and 10 cm 1 , respectively). The elevational position y of the source is fixed at 5.70 mm, which corresponds to the fiber positions in our experimental transducer array (see Fig. 2). For extremely small = µ ( 2 cm ) s 1 , Monte Carlo produces a spike since ballistic transport crosses the image plan. As expected, Model I is very close to Monte Carlo results for diffusive conditions because scattering dominates in the medium. Model II approximates Model I when the measurement point is located at distances from the source and the interface much larger than the transport mean free path l t . The higher the scattering, the closer the agreement between Model II and Monte Carlo, as shown in Fig. 3(d) and (e). We also compared fluence models at different lateral source positions (y is 2.85 mm, 5.70 mm and 11.40 mm, see Fig. 3(f-h)). Scattering µ s was set to 10 cm 1 . The position z max of maximum optical fluence along the z axis increases with increasing distance y between the image plane and light source, and Model II converges quickly to Model I soon after z max . Fig. 3(b) summarizes the shift in z max as a function of source position y for different µ s typical of biological tissue [8]. For the source position in our system ( = y 5.7 mm), z max barely changes for µ s varying from 5-20 cm 1 . This means that simplified Model II converges to MC and Model I at very similar depths over a wide range of tissue scattering. For our transducer configuration, z max is around 3 mm. This fact is very important for practical fluence assessment because µ s is not known a priori.

Simulation results
Based on these initial simulations, it is clear that Model I can be used as ground truth if the light source is located more than a distance l t from the light source. We now focus on the difference in parameter estimation error (bias) between Models I and II in more detail. The  Fig. 4(a) shows the error over z l / t for several µ s over the tissue range. The error is under 10 % when z l / t is larger than ∼10. Fig. 4(b) shows the error for several assumed µ a values. Note that the µ a variation rarely changes the error. Fig. 4(c) illustrates the fractional error z ( , µ ) s for < < µ 5 35 cm s 1 and < < z 0 40 mm when Model II is used to estimate µ eff . Fig. 4(d) displays pixel values of (c) over the z l / t (= zµ s ) axis. Fig. 5 presents optical parameter estimates using Models I and II when data are contaminated by noise. We used Eq. (13) to set groundtruth scattering coefficient variations over wavelength (715−875 nm). Fig. 5(a-c) show estimates when the SNR is the highest (50 dB). Fig. 5(d) shows the correlation between the estimated optical fluence spectrum and ground truth over a wide range of SNR. Model I estimation error (bias and deviation) is only a function of measurement noise. Estimates of µ eff are unbiased at high SNR while those of µ s are biased and their errors are high even at this high SNR level. However, despite large error in µ s estimated by Model I, the correlation between resultant fluence estimates and ground truth is close to 1. As shown in Fig. 5(d), the lower SNR level causes lower correlation because both estimation deviation and bias are higher. In Model II, estimation error is caused by model discrepancy as well as measurement noise. Note that Model II underestimates µ eff at high SNR, but the ratio of estimate μ eff to ground truth µ eff remains almost constant over the wavelength range. Thus, it provides competitive correlation between fluence estimate and ground truth despite the bias. The lower SNR, of course, leads to lower estimation performance (see Fig. 5(d)).

Experiments
Phantom experiments produce realistic PA measurements to help validate the proposed methods. We conducted two studies.

Phantom study I
In the first study, we explored the accuracy of optical fluence estimates in a turbid medium using a human hair as the absorbing target. Each hair was aligned parallel to the transducer elevational direction (y-axis) so that it appears as a point target in the z-x image plane (see Fig. 6(a)). The target was positioned in a cubic tank (open on top) filled with an optical scattering medium, and a transducer was positioned at the medium surface. The tank was sufficiently large so that boundary effects from all faces except the top interface between media can be neglected. The optical medium was an intralipid solution (20 % IV fat emulsion, Fresenius Kabi, Deerfield, USA) with homogeneous scattering. We diluted the original 20 % emulsion to 0.5 %-4 % ones to control the scattering coefficient.
Radio frequency ultrasound data were recorded using the Verasonics system and a 128 element linear-array transducer (LA 15/ 128-1633, Vermon S.A. Tours, France). The transducer center frequency is 15 MHz and the 3-dB bandwidth is 11−19 MHz. The array element pitch is 0.1 mm. The laser pulse energy was between 0.4 mJ and 0.5 mJ at the tissue surface, depending on the wavelength. A little (4% portion of it) was taken from the main path by a beam splitter and recorded for every laser pulse with a photodetector located before the fiber bundle. Thus, the energy of every laser pulse was measured to compensate its dependence on wavelength and take into account pulse- to-pulse energy variations (∼ 8%). The number of axial samples and scan-lines were adjusted to the field of view (25.2 mm × 12.7 mm). One PA data set consists of 2048 samples × 128 channels × 20 fibers × 10 wavelengths × 40 frames acquired at a 1 kHz repetition rate, representing a total acquisition period of 8 s. We averaged every data set over frames to enhance SNR and processed PA signals using  conventional delay-and-sum beamforming to obtain an image (512 × 128) for every fiber and wavelength. A Hilbert transform was then applied to obtain the smooth envelope of axial image samples. The enveloped image was used as y j k , in Eq. (8) for the k th fiber and j th wavelength. Fig. 6(b) illustrates optical parameter estimates over the range of laser wavelengths (715−875 nm) for the 0.5 % intralipid solution as the scattering medium. Each marker and error bar indicate mean and standard deviation of 8 samples, respectively. The reference spectrum employed here was obtained from Ref. [24]. Since this spectrum was measured for a 1.2 % intralipid concentration in [24], it was scaled here to compare with the 0.5 % solution and plotted as dotted lines (Fig. 6(c)). The expected reduced scattering coefficient shown on the right is presented as a shaded region since reference papers [24,[34][35][36] give different measurement values. The scattering coefficient µ s is under 5 cm 1 over the measured wavelength range, close to the lower edge of medium scattering used in numerical simulations. As shown in Fig. 6(b,c), estimated µ s , and µ eff match very closely with those reported in the literature under similar measurement conditions. Due to the small µ s , Model II provides a less accurate spectral estimate. Estimates of effective attenuation coefficient for three higher intralipid concentrations (0.5 %, 1% and 2%) are presented in Fig. 6(b), (d) and (f) respectively. As expected, both Models I and II provide very accurate assessment of µ eff , whereas µ s estimates from Model I are clearly not accurate for these higher scattering media (Fig. 6(e,g)). However, as shown in the numerical simulations, one parameter, µ eff , is enough when Model II is valid. Thus, inaccuracies in the reconstructed µ s do not significantly affect the accuracy of laser fluence estimates.

Phantom study II
The purpose of the second study is to (i) perform optical fluence estimation in the medium, and (ii) apply optical fluence corrections to spectroscopic PA measurements. For this test, the container held three cylindrical tubes aligned parallel to the y-axis (see Fig. 7(a)). We injected a nanoparticle solution (gold nanorods (GNR), width (11.4 nm), length (44.8 nm), mass concentration (2.2 mg/mL), longitudinal peak (776 nm), NanoHybrids Inc. Austin, USA) [25] and a black ink solution (Higgines Black Magic Ink, Chartpark Inc., Leeds, USA) as absorbers in Tubes I and III, respectively, where the absorption spectra of the solutions are well-known. Also, we injected water into Tube II as a control. The maximum absorption coefficients of the nanoparticle solution and black ink solution are roughly 35 cm 1 and 82 cm 1 , respectively, over the wavelength range (715 nm to 875 nm). The container was filled with a 1% intralipid solution as the scattering medium. We additionally added customized Prussian blue ink [25] to the solution to increase spectral distortion at the expense of high attenuation (low data SNR).
We used a total of 8 datasets and averaged them to improve SNR. Fig. 7(a) presents PA images at three particular wavelengths, where every pixel value is proportional to the photoacoustic signal magnitude presented on a log scale. Due to the limited view and bandwidth of the transducer, the signal only appears at the top and bottom of the tube. The signal beneath the tubes is from a reverberant wave in the tube. The region near Tube II has a weak signal because the tube material itself weakly absorbs light over this wavelength range. The ranges of SNR in Tube I and Tube III are 27.8-39.9 dB and 41.1-47.5 dB, respectively, over the wavelength range. Fig. 7(b) and (c) show estimation results using signals in Tube I and Tube III for laser fluence assessment, respectively. In other words, we demonstrate here that any target, independent of its absorption properties (for example, for a solution GNR or black ink), effectively absorbing light over the spectral range probed by the system can be used for laser fluence estimation in the swept-beam concept. In practice, such targets can be blood vessels, injected contrast agents or labelled drugs, or other absorbing targets.
The first and second row in columns (b) and (c) of Fig. 7 show effective light attenuation and reduced scattering coefficient estimates, respectively. We smoothed all estimates over wavelength to reduce bias and applied them for fluence compensation. Note that the signals in Tube III provide more stable estimates than those in Tube I due to higher SNR, as predicted in numerical simulations above. The third and fourth row compare raw PA spectra and spectra after fluence correction for nanoparticles and black ink, respectively. Corrected spectra more closely reflect true absorption spectra and approach the ground truth. As predicted in simulations, high SNR increases correction accuracy. Note that laser fluence estimation was performed here using only PA signals from a single target, yet reconstruction accuracy is already reasonable. If a larger number of points within the PA image are used simultaneously for fluence correction, then more accurate corrections can be obtained without the need for signal averaging.

Discussion and conclusions
In simulations, we investigated optical parameter estimation using measurements on chromophores acting as targets in a turbid medium. Although we did not fully validate over a wide range of absorption target types, simulations and phantom studies presented here yield consistent estimation results. We can clearly infer that estimation using more sparse targets can lead to lower bias and variance at any noise level if their signal strengths are similar. Also, when multiple targets are used for estimation, targets located near image edges will contribute more to estimation results than those near the center because they provide a wider range of distances r r | | k between target and source positions.
We approximated the narrow beam as an ideal pencil beam. Even though every fiber produces a diverging beam corresponding to their NA = 0.22, we found through MC simulations that the fractional error in the fluence with respect to the true fluence (e.g. for the diverging beam of the same diameter at the medium surface) is less than 3% in the image domain. This minor discrepancy does not affect estimation performance.
The primary advantage of Model I is that it more closely matches the results of MC simulations when the source-target distance exceeds the photon transport mean free path. Thus, estimates of µ eff based on this model are unbiased and accurate at zero noise levels, as shown in Fig. 6. A disadvantage of this model is the additional parameter, µ s , that must be simultaneously estimated. A two-parameter search often produces higher estimation errors under real, not ideal, conditions, especially when experimental measurements are contaminated by noise. Thus, using constraints for narrowing the search range of µ s based on prior knowledge would increase the performance under low SNR conditions.
Model II is an approximate version of Model I where µ s can be considered high. In reality, this approximation works when a target is located far (i.e. at distances much larger than l t ) from all light sources and interfaces. As shown in Figs. 3 and 4, within the range of typical light scattering of biological tissue, Model II can be used for targets located at depths deeper than 5−7 mm from the medium surface, i.e. at depths where the wavelength dependence of laser fluence starts affecting spectroscopic measurements. As depth z increases, Model II converges to Model I very quickly. Note, that the ratio of μ eff to µ eff remains mostly constant as µ eff varies over wavelength, as shown in Fig. 5, even though we used the highest variation of µ s (for brain tissue) from all biological tissue reported in the literature. Model II is a single parameter model that provides fast and stable laser fluence estimates for the range where it is valid. In our signal model, we assumed that the only factor creating PA intensity variation over fibers is the fluence. However, during data acquisition using our system, every fiber delivered slightly different laser pulse energy due to positioning deviations and pulse-to-pulse laser emission fluctuations. Thus, we took into account pulse energy variations with simultaneous measurements of laser pulse energy for every laser pulse to prevent additional estimation biases. Phantom studies validated the overall performance of spectral corrections based on optical parameter estimation. Fluence-corrected spectra closely match ground-truth spectra, as shown in Figs. 6 and 7. Estimation variance depends on SNR at each wavelength. Thus, estimates based on the GNR target are less accurate because of relatively small signal amplitude at the edges (715 nm and 875 nm wavelengths) of the GNA absorption spectrum.
In our studies, we neglected shadowing (laser fluence at one absorber shadowed by another absorber) because absorbers are located sparsely over the imaging field. However, in clinical applications, the presence of large-scale absorbers such as arterial vessels can induce unavoidable fluence attenuation in the entire image field. We conducted simple MC simulations to explore this effect ( Fig. 8(a, b)) assuming a cylindrical-shaped vessel close to a point target in a medium (µ eff = 0.9 cm 1 ). When the vessel diameter is 2.4 mm, our estimation method using noise-free measurements from the target position yields µ eff = 1.57 cm 1 . Note that our primary purpose is estimating true fluence rather than estimating true optical parameters of the medium. As shown in Fig. 8(c), the actual in-depth fluence attenuation, calculated through the target (as indicated in the upper left panel of Fig. 8(a)), is much closer to our estimate than to the fluence without the vessel. Indeed, the best fit to the ground-truth using the analytical fluence model (blue dotted line) yields µ eff = 1.41 cm 1 (using the depth range 5−10 mm). Thus, our estimates are not perfect, but they account for large absorbers and provide accurate compensation of fluence variations within the medium. Our future studies will develop more detailed algorithms to further improve reconstruction accuracy.
Our system emits around 1 mJ pulses at 1 kHz. This complies with the maximum permissible exposure (MPE) for skin [37] over the Fig. 7. Spectroscopic PA imaging of absorbing targets in a turbid medium. (a) Examples of PA images displayed on a log (dB) scale at 3 different wavelengths (715 nm, 795 nm and 875 nm) for Phantom Study II. The first, second and third tubes were filled with GNR solution, water and black ink, respectively. The scattering medium was a solution of intralipid and Prussian blue ink. The pixel value of the image is associated with the PA pressure. Dotted circles in the image indicate tube cross-sections. (b) and (c) columns show parameter estimates using nanoparticle signals in Tube I and black signals in Tube III, respectively, and spectrum correction results using the estimates. First and second rows represent estimates of effective light attenuation and reduced scattering coefficients in a turbid medium, respectively. Points marked (o) and ( ) denote estimates using Model I and Model II, respectively. The dotted line is obtained by smoothing estimates over wavelength. Third and fourth rows show absorption spectra of the nanoparticle and black ink, respectively. Magenta curves with points marked (◇) correspond to measured PA absorption spectra, i.e. without fluence correction. Orange and blue curves with points marked (o) and ( ) respectively correspond to ink and GNR absorption spectra obtained after applying fluence correction using Model I and Model II, respectively. Dotted lines indicate the reference target spectra measured independently with optical spectrophotometry.
wavelength range (700-900 nm) given a small standoff between fiber and skin. During data acquisition, the maximum laser pulse energy at the fiber output was less than 1 mJ due to ∼ 60 % efficiency of light coupling into optical fibers. Coupling efficiency can be improved and will be a focus of future studies.
In Phantom Study II, we used signal averaging to increase SNR because of the high optical attenuation of the medium. Since most human soft tissues have lower attenuation [8,33], our compensation method should require much lower averaging to obtain comparable results for in vivo studies. Also, as we have shown, any absorber independent of its absorption spectrum can be used for fluence reconstruction. Therefore, the larger the number of image points participating in fluence estimation, the better the accuracy and stability, and fewer averages will be required. In addition, our ongoing research on noise reduction using deep-learning may potentially reduce the need for averaging [23].
In conclusion, our study shows that one of the fundamental problems of PA imaging, i.e. decoupling the local light absorption spectrum from wavelength-dependent optical fluence, may be solved automatically using a swept-beam imaging concept at least when the medium under study can be considered macro-homogeneous. We adopted an analytic laser fluence model for a homogenous scattering medium and estimated wavelength-dependent fluence variations within the image plane from PA measurements without knowledge of optical properties and without large computational costs. In phantom studies, we have shown that fluence correction considerably improves the accuracy of measured chromophore spectra. Real-time PAUS systems providing fluence-compensated spectroscopic PA imaging have the potential to enable molecular imaging for many clinical applications, such as interventional procedure guidance using molecularly labelled therapeutic agents and procedure validation based on spectroscopic confirmation of modifications in microvascular networks.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

A.1 Refractive-index-mismatched boundary
Assume that the boundary between medium I and II is refractive-index-mismatched and one source to medium I is located at the interface. Outgoing light toward II from I is partially reflected at the boundary due to the mismatch. This can be mathematically expressed as where L and R denote the radiance and Fresnel reflection, respectively [26]. The vectors r, s and n denote the position on the boundary, unit direction vector, and unit normal vector pointing toward medium I, respectively. d denotes a differential solid angle element. Using a spherical harmonics expansion, the radiance can be approximated as = + L r s r r s ( , ) where J is the current density, or the energy flow per unit area. The Fresnel reflection is given as