1.IntroductionConventional spectrophotometers are probably the most commonly used type of instruments for material analysis by determination of absorbance A or attenuation coefficient as a function of wavelength for a given sample. For turbid materials including biological tissues, additional optical parameters are needed to characterize molecular composition and particle sizes on scales close to by their ability to absorb and scatter light. One paradigm is to retrieve these parameters by the radiative transfer (RT) theory,1 and its spectroscopic implementation is termed multiparameter spectrophotometry (MPS). The RT theory quantifies light–matter interaction by the parameters of absorption coefficient and scattering coefficient in addition to a single-scattering phase function with and as the polar and azimuthal scattering angles. MPS thus requires measurement of multiple light scattering signals to solve multiple inverse scattering problems (ISPs) at selected values of that remains challenging for highly turbid samples. To reduce complexity, of an unknown sample is often modeled by an analytical function of pHG(cosθ) proposed by Heyney and Greenstein under the assumption of axially symmetric scattering.1–3 Alternative model functions have been studied for but is the preferred one with a form fully specified by an anisotropy factor g as the mean value of cos . By choosing pHG as the phase function, MPS is to measure multiple light scattering signals and determine inversely the RT parameters as function of at . These parameters can be expressed as a vector of . A widely used approach of MPS is to measure three signals of collimated transmittance , hemispherically integrated diffuse reflectance and transmittance with one or two integrating sphere.3–11 If only one integrating sphere is used, signals need to be measured in two steps with one for and another one for with acquired in either step. To solve the ISPs, one first derives () from by the Beer-Lambert law followed with retrieval of and from and using different methods of forward modeling guided by gradient descent based inverse solver. The modeling methods include numerically solving the RT boundary-value problem, adding-doubling algorithm or Monte Carlo (MC) simulations. Despite its popularity as a research tool, the need for two integrating spheres or two steps of signal measurement makes it difficult to translate this approach into an easy-to-use instrument like a conventional spectrophotometer. In addition, the use of integrating sphere limits severely the accessibility of the approach to non-specialists due to time-consuming sample assembly and system maintenance. Other approaches without integrating sphere have been investigated such as goniometric measurement and detection of and non-hemispherical and at fixed angles with MC based inverse algorithms.12,13 We have previously shown that the RT parameter vector can be uniquely determined from three simultaneously measured signals of non-hemispherical diffuse reflectance , diffuse transmittance and forward transmittance without integrating sphere.14,15 MC simulations are performed to accurately calculate signals as , and and repeated to match the measured ones. A gradient descent algorithm has been developed to guide iteration in the RT parameter space to minimize an objective function defined as the sum of squared percentage differences between the measured and calculated signals. An ISP at is deemed as solved with when in which represents a threshold based on the experimental errors in signal measurement. The inherent statistical variance in the calculated signals, however, makes it difficult to accurately solve ISPs for highly turbid and optically thick samples by gradient decent despite the existence of a unique solution.16 Specifically, sizable regions of small values exist in the RT parameter space of P for these ISPs of large scattering albedo and optical thickness with as sample thickness. In such regions, values of fluctuate considerably due to the variance of MC simulations that often lead to errors in the solution given by . To solve these challenging ISPs for highly turbid samples, one needs to either significantly increase the number of photons in MC simulations for variance reduction or search manually by contour analysis of distributions in the RT parameter space. Either way gives rise to high computational cost and prevents rapidly solving ISPs for MPS.16 In this report, we present a rapid inverse solver for MPS based on a particle swarm optimization (PSO) algorithm with a novel objective function through local averaging in the space of to determine .17 The function significantly reduces the effect of MC simulation variance in calculated signals on inverse calculation by PSO and computation time to solve ISPs by reducing the number of photons in MC simulations. Our validation results with 20% intralipid samples demonstrate that can be retrieved rapidly between 520 and 1000 nm by executing MC simulations on one GPU board. 2.Materials and Methods2.1.Signal Measurement by MPSAn experimental system has been constructed to measure , , and for validation of the new inverse solver with the signal detection configuration shown in Fig. 1(a). The details of the experimental system were previously reported.18 Briefly, a xenon light source (XL1-175-A, WavMed Technologies Corp.) and a monochromator (CM110, CVI Corp.) are employed to produce a monochromatic beam with λ adjustable between 520 and 1000 nm in steps of 20 nm and bandwidths around 5 nm. The beam is modulated at a frequency of by a mechanical chopper (SR540, Stanford Research Systems) and incident on an assembly consisting of a turbid sample confined in a spacer ring between two glass slides. The intensity of the incident beam is monitored by a photodiode of (FDS1010, Thorlabs, Inc.). Three photodiodes (FDS100, Thorlabs, Inc.) of to are used to measure respectively for diffusely reflected, for diffusely transmitted and for forwardly transmitted light intensity. The current signals of photodiodes were amplified by an in-house built four-channel lock-in amplifier to obtain the measured signals of and . Figure 1(a) shows the detection configuration for acquisition of measured signals from the sample assembly. 2.2.Signal Calculation by iMCAn in-house developed individual photon tracking MC (iMC) code was employed to calculate signals as and from given for a phantom of the same shape as the sample inside a spacer between two glass slides by tracking N0 photons, which imports the parameters of sample size and detection configuration as input data.16,19,20 The code injects each of the N0 photons incident on the sample assembly and then tracks the photon once it transports inside a glass slide or sample until it is either absorbed inside the sample or escapes into air. The exit location and propagation direction of an escaping photon on a glass slide surface are used to determine if it hits a detector for detection. A counter associated with each detector records the number of detected photons as , and by the detector , and , respectively. The above process repeats until the total number of injected photons reaches and the calculated signals are given by , and . Because of independence in trajectory among the photons, the numbers of detected photons follow Poisson distributions.21 To estimate the variance in these photon numbers, one can draw a random number of Poisson distribution with as probability mass function and as the mean. In the case of calculated by an iMC simulation, equals to if distribution of follows and yields the variance-free value of by tracking “infinite” number of photons. One thus can use to quantify the effect of variance on calculated signals and optimize the objective function for variance reduction. 2.3.PSO AlgorithmA stochastic algorithm based on PSO has been developed to guide iMC simulations and solve for in the RT parameter space from measured signals by optimizing an objective function as shown in Fig. 1(b). An objective function quantifies the difference between measured signals and calculated ones obtained by given , and optimization of this function reduces the difference to solve for from an initial choice of . The PSO algorithm was chosen for its high efficiency and ability to perform global search or avoid local traps in the RT parameter space. A search proceeds through multiple threads in PSO, and the threads are represented by a swarm of “particles” with index and as the number of threads or particles. Here, a particle symbolizes a thread of positions in the RT parameter space along which it moves from to as22 where ; is the number of iterations; is the particle’s velocity; and ’ are random numbers of uniform distribution from 0 to 1; is the particle-best position; is the swarm-best position; , and are weights of terms contributing to ; and is a damping constant to increase stability. A particle transits toward and through completed iterations by Eq. (1) that respectively yields the best position of that particle and all particles for optimized objective function to reduce difference between calculated and measured signals. By setting values of and equal and large in comparison to in Eq. (1), for example, enables fast convergence of all particles from their best positions towards on average as j increases. A region in the RT parameter space is designated as the search space and in Eq. (1) is set to if the former moves out of . The search stops once the objective function at reaches a preset threshold or exceeds and current is saved as for output.3.Results and Discussion3.1.Effect of Variance in Calculated Signals on Solving ISPsMPS requires solving one ISP for each value of from the measured signals. The first step is to perform forward calculations of signals from given by iMC simulations of light–matter interaction in the sample. The calculated signals can be normalized by respective measured signals and expressed as a vector of . An inverse algorithm is to iterate iMC simulations toward the objective of making as close to as possible. A function of squared error sum is often selected as the objective function to guide iteration and decide when stops. Ideally, vanishes for a perfect ISP solution given by . In practice, one deems an ISP as solved if with representing a threshold determined by the error level in signal measurement. To quantify the effect of variance in calculated signals by a stochastic iMC simulation, we define a fluctuation vector for calculated signals as where , , and are the variance-free signals calculated by tracking “infinite” numbers of photons. It should be noted that is independent of and one can express by with and “o” denoting the vector form of component-wise product. Combining these results, can now be written as where is the variance-free value of . It is clear that if represents a “true” solution of ISP.3.2.Comparison of Two Objective Functions for Variance SuppressionBased on above analysis, local averaging on in the RT parameter space was first employed to suppress variance and hence increase the stability of inverse solutions. We then compared two objective functions defined below on their abilities to further reduce the effect of variance where is a set of vectors consisting of and its neighbors on a 3D grid of . The components of as defined in Eq. (2) can be obtained as by drawing a random number q independently from and setting the values, related to , and , equal to focus on the effect of variance for simplicity. We further set with and its nearest neighbors as the members of and assumed that all of them have the same value of . The values of and are plotted in Fig. 2(a) against ranging from 0.1% to 1% with a solid line representing or . These results show clearly that fluctuates significantly less than for between 0.6% and 1%, which is the typical range of set by the error level of measured signals.16To gain further insight on fluctuation in calculated signals, the relative difference Δ between and or and are calculated and averaged as and where is uniformly sampled between 0.1% and 1.1% in each sum and is the total number of samples in the sums. The sampled values of correspond to the mean values of each signal term in ranging from 1.8% to 6.1% since is defined as the sum of squared relative differences of calculated and measured values on three signals. Figure 2(b) illustrates the dependence of the two Δ functions on values with set to 1000. As discussed before, represents the mean value of which is the probability mass function of Poisson distribution for modeling detected photon numbers, such as and related calculated signal Therefore, the ratio of equals to as the variance-free value of in the above case and larger corresponds to larger . The results in Fig. 2(b) demonstrate that use of instead of or as the objective function can significantly reduce the variance of calculated signals for the same value of . Alternatively, can be considerably reduced for the same variance to speed up forward calculations if is adopted. For example, computational time of iMC simulations can be cut in half with ), instead of , as the objective function by tracking half of photons with similar variance.3.3.Measurement of 20% Intralipid SampleTo validate the new approach, RT parameters of 20% intralipid (I141-100ML, Sigma-Aldrich) have been determined from three measured signals of , and for from 520 to 1000 nm in steps of 20 nm. Figure 3 presents the RT parameters of a sample with thickness obtained by the PSO based algorithm with selected as the objective function. Signal measurement was repeated three times to obtain the mean values and standard deviations which are plotted in Fig. S1 in the Supplementary Material. Previously determined values of real refractive index of the 20% intralipid were used to obtain interpolated values in iMC simulations at the wavelengths of measurement for the calculated signals and is presented in Fig. S2 in the Supplementary Material.23 The same measurements were repeated on another sample of and the inversely determined RT parameters agree well with those in Fig. 3 within the experimental errors. The search region for PSO based inverse algorithm at each wavelength was set between and for , 10.0 and for , 0.10 and 1.00 for . In additional to choosing for the particle number and , we set , , and for the PSO parameters defined in Eq. (1) for optimized performance of inverse calculations. The initial positions of for all particles were randomly distributed in each of nine equal partitions of to ensure globally optimized solutions of ISPs. We found that the convergence toward the final solution was not affected by the choice of under the above condition. After calculation of with , iterations as defined in Eq. (1) were performed in . Specifically, was set as for particle and as that has the largest value among all particles for to obtain and by Eq. (1). The search then continued as j increases to 2, 3,… and stopped when either exceeds for or with saved as and output together with and . The threshold corresponds to the value of in Eq. (5) with set to 0.8% and to 7. It took about eight iterations on average to solve an ISP of from the measured signals of 20% intralipid sample. The total numbers of iMC simulations were found to be about 800 for solving each ISP if 27 particles were employed in PSO based search. Signal were calculated by iMC simulations in which was gradually increased from to as became 1% or less. On average, over different values of RT parameters and , each iMC simulation took about 0.2 s and solving per wavelength took about 160 s by executing the code on one GPU board (Nvidia, GeForce RTX 2080 Ti). The results of in Fig. 3 compare similarly to those obtained by methods using integrating sphere.8,24 For example, both of the range and wavelength dependence of and g agree well with those in Ref. 8 for the 20% intralipid sample while the range of is between those of Ref. 24 for 10% intralipid and Ref. 8. Since the 20% intralipid is high turbid with scattering albedo over the measured wavelengths, the large differences in the values of can be attributed to the very large errors in determining with values smaller than by three orders of magnitude as previously discussed.8 We note further that the values of optical thickness and albedo of the intralipid sample reach the maximum values of 7.1 and 99.98%, respectively for . Results of our study on samples of 20% intralipid with sufficiently large values showed that current MPS method with the inverse solver reported here can yield unstable solutions of when the value of τ becomes larger than 7.1 for very large values of a above 99.9990%. A detailed comparison of the measured and calculated signals and analysis of background noise in signal measurement suggests that the instability of inverse solution is mainly caused by the variance in calculated signals, which may be mitigated without using much larger than those used for ISPs of for iMC simulations. A study is underway to further improve the PSO based inverse solver by taking into account the distribution obtained from completed iterations, which is expected to enable accurate measurement of RT parameters for highly turbid samples with up to 15 without significant increase of computational cost. Such an improvement can make the approach of MPS concerned here capable of characterizing optically thick samples by their RT parameters for which the conventional approach with integrating sphere fails for inability to accurately measure collimated transmittance . 4.ConclusionsWe have analyzed the effect of variance by MC simulations on solving ISPs and proposed a novel objective function to reduce variance in calculated signals and speed up simulations. Combination of the objective function with the PSO algorithm leads to a rapid inverse solver for determination of RT parameters of turbid samples from three measured signals without integrating sphere. The inverse solver has been validated by the results of RT parameter retrieval on samples of 20% intralipid in a wavelength range of 520 to 1000 nm with ISPs solved rapidly using one GPU board. Taken together, we have demonstrated a new approach of MPS, which allows its translation into a powerful and easy-to-use instrument for clear separation and characterization of molecular composition and turbidity. Code and Data AvailabilityAll data in support of the findings of this paper are available within the article or as supplementary material. Code and other materials associated with this article are available upon request sent to the corresponding author. AcknowledgmentsAuthors thank Dr. Kenneth M. Jacobs for help on development of the experimental system. J. Jing acknowledges grant support by Education Department of Hunan Province of China (Grant No. 22B0680). ReferencesH. C. van de Hulst, Multiple Light Scattering: Tables, Formulas and Applications, Academic Press, New York
(1980). Google Scholar
B.-H. Gao et al.,
“Efficient equation-solving integral equation method based on the radiation distribution factor for calculating radiative transfer in 3D anisotropic scattering medium,”
J. Quant. Spectrosc. Radiat. Transf., 275 107886 https://doi.org/10.1016/j.jqsrt.2021.107886
(2021).
Google Scholar
S. Y. Jeong et al.,
“Measurements of scattering and absorption properties of submillimeter bauxite and silica particles,”
J. Quant. Spectrosc. Radiat. Transf., 276 107923 https://doi.org/10.1016/j.jqsrt.2021.107923
(2021).
Google Scholar
S. A. Prahl, M. J. C. van Gemert and A. J. Welch,
“Determining the optical properties of turbid media by using the adding-doubling method,”
Appl. Opt., 32
(4), 559
–568 https://doi.org/10.1364/AO.32.000559 APOPAI 0003-6935
(1993).
Google Scholar
L. Wang, S. L. Jacques and L. Zheng,
“MCML-Monte Carlo modeling of light transport in multi-layered tissues,”
Comput. Methods Program Biol., 47
(2), 131
–146 https://doi.org/10.1016/0169-2607(95)01640-F
(1995).
Google Scholar
Y. Du et al.,
“Optical properties of porcine skin dermis between 900 nm and 1500 nm,”
Phys. Med. Biol., 46
(1), 167
–181 https://doi.org/10.1088/0031-9155/46/1/312 PHMBA7 0031-9155
(2001).
Google Scholar
X. Ma et al.,
“Bulk optical parameters of porcine skin dermis tissues at eight wavelengths from 325 to 1557 nm,”
Opt. Lett., 30
(4), 412
–414 https://doi.org/10.1364/OL.30.000412 OPLEDP 0146-9592
(2005).
Google Scholar
C. Chen et al.,
“A primary method for determination of optical parameters of turbid samples and application to intralipid between 550 and 1630 nm,”
Opt. Express, 14
(16), 7420
–7435 https://doi.org/10.1364/OE.14.007420 OPEXFF 1094-4087
(2006).
Google Scholar
P. Lemaillet et al.,
“Correction of an adding-doubling inversion algorithm for the measurement of the optical parameters of turbid media,”
Biomed. Opt. Express, 9
(1), 55
–71 https://doi.org/10.1364/BOE.9.000055 BOEICL 2156-7085
(2018).
Google Scholar
B. W. Xie et al.,
“Experimental study of the radiative properties of hedgehog-like ZnO–Au composite particles,”
J. Quant. Spectrosc. Radiat. Transf., 232 93
–103 https://doi.org/10.1016/j.jqsrt.2019.05.006
(2019).
Google Scholar
F. Bergmann et al.,
“Precise determination of the optical properties of turbid media using an optimized integrating sphere and advanced Monte Carlo simulations. Part 2: experiments,”
Appl. Opt., 59
(10), 3216
–3226 https://doi.org/10.1364/AO.385939 APOPAI 0003-6935
(2020).
Google Scholar
C. Y. Ma et al.,
“GPU-accelerated inverse identification of radiative properties of particle suspensions in liquid by the Monte Carlo method,”
J. Quant. Spectrosc. Radiat. Transf., 172 146
–159 https://doi.org/10.1016/j.jqsrt.2015.08.002
(2016).
Google Scholar
X. Liang et al.,
“Spectrophotometric determination of turbid optical parameters without using an integrating sphere,”
Appl. Opt., 55
(8), 2079
–2085 https://doi.org/10.1364/AO.55.002079 APOPAI 0003-6935
(2016).
Google Scholar
P. Tian et al.,
“Spectral determination of μa, μs and g from single and multiple scattering signals with one optically thick sample,”
J. Quant. Spectrosc. Radiat. Transf., 245 106868 https://doi.org/10.1016/j.jqsrt.2020.106868
(2020).
Google Scholar
P. Tian et al.,
“Multiparameter spectrophotometry platform for turbid sample measurement by robust solutions of radiative transfer problems,”
IEEE Trans. Instrum. Meas., 70 6003110 https://doi.org/10.1109/TIM.2020.3032184 IEIMAO 0018-9456
(2021).
Google Scholar
Y. Qin et al.,
“Robustness of inverse solutions for radiative transfer parameters from light signals measured with different detection configurations,”
J. Quant. Spectrosc. Radiat. Transf., 274 107883 https://doi.org/10.1016/j.jqsrt.2021.107883
(2021).
Google Scholar
M. Clerc and J. Kennedy,
“The particle swarm - explosion, stability, and convergence in a multidimensional complex space,”
IEEE Trans. Evol. Comput., 6
(1), 58
–73 https://doi.org/10.1109/4235.985692 ITEVF5 1089-778X
(2002).
Google Scholar
Z. Jones et al.,
“Study of inverse solution for multiparameter spectrophotometry by three photodiodes,”
Proc. SPIE, 12376 1237604 https://doi.org/10.1117/12.2651470 PSISDG 0277-786X
(2023).
Google Scholar
X. Chen et al.,
“Fast method for inverse determination of optical parameters from two measured signals,”
Opt. Lett., 38
(12), 2095
–2097 https://doi.org/10.1364/OL.38.002095 OPLEDP 0146-9592
(2013).
Google Scholar
P. Tian et al.,
“Quantitative characterization of turbidity by radiative transfer based reflectance imaging,”
Biomed. Opt. Express, 9
(5), 2081
–2094 https://doi.org/10.1364/BOE.9.002081 BOEICL 2156-7085
(2018).
Google Scholar
G. S. Fishman, Monte Carlo: Concepts, Algorithms, and Applications, Springer-Verlag(
(1996). Google Scholar
F. Marini and B. Walczak,
“Particle swarm optimization (PSO). A tutorial,”
Chemometr. Intell. Lab. Syst., 149 153
–165 https://doi.org/10.1016/j.chemolab.2015.08.020
(2015).
Google Scholar
H. Ding et al.,
“Determination of refractive indices of porcine skin tissues and intralipid at eight wavelengths between 325 and 1557 nm,”
J. Opt. Soc. Am. A, 22 1151
–1157 https://doi.org/10.1364/JOSAA.22.001151 JOAOD6 0740-3232
(2005).
Google Scholar
S. T. Flock et al.,
“Optical properties of intralipid: a phantom medium for light propagation studies,”
Lasers Surg. Med., 12
(5), 510
–519 https://doi.org/10.1002/lsm.1900120510 LSMEDI 0196-8092
(1992).
Google Scholar
BiographyJiahong Jin received his BS and MS degrees in mathematics and optics from Sun Yat-sen University, China, and his PhD in biomedical physics from East Carolina University, United States. He joined the Physics Faculty, Hunan Institute of Science and Technology (HNIST) in 2014 and the Institute for Advanced Optics of HNIST in 2016. He is conducting research in optical imaging of turbid materials and biological tissues and machine learning study on diffraction image data of cells. Zachary D. Jones received his BS degree in applied physics from Appalachian State University in the United States and is currently a PhD candidate for the degree in biomedical physics at East Carolina University. He has experience in medical imaging device research and development as a research assistant with Perfusio Corp. His research interests include tissue optics and Monte Carlo modeling of light transport in tissue phantoms. Jun Q. Lu received her BS and MS degrees in physics from Nankai University, China, in 1983 and 1986, respectively and her PhD in physics from the University of California at Irvine, United States in 1991. She joined the Faculty of Physics, East Carolina University, United States, in 1996. Her research areas included theoretical modeling, numerical simulations, and machine learning-based analysis of light scattering data acquired from biological cells and turbid samples. Xin-Hua Hu received his BS and MS degrees in physics from Nankai University in China, his MS degree in physics from Indiana University at Bloomington, and his PhD from the University of California at Irvine in 1991. He joined the Physics Faculty of East Carolina University in 1995 and established the Biomedical Laser Laboratory. He conducts and directs research programs in optical imaging and machine learning study of light scattering by cells and in turbid materials including biological tissues. |
Simulations
Particle swarm optimization
Photons
Spectrophotometry
Integrating spheres
Radiative transfer
Turbidity