2.1.

Signal Measurement by MPS

An experimental system has been constructed to measure $R_d(\lambda )$, $T_d(\lambda )$, and $T_f(\lambda )$ 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.¹⁸ 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 $f_0=370\mathrm{Hz}$ 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 $I_0$ is monitored by a photodiode of $D_1$ (FDS1010, Thorlabs, Inc.). Three photodiodes (FDS100, Thorlabs, Inc.) of $D_2$ to $D_4$ are used to measure respectively $I_{\mathrm{Rd}}$ for diffusely reflected, $I_{\mathrm{Td}}$ for diffusely transmitted and $I_{\mathrm{Tf}}$ 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 $R_d=I_{\mathrm{Rd}}/I_0,T_d=I_{\mathrm{Td}}/I_0$ and $T_f=I_{Tf}/I_0$. Figure 1(a) shows the detection configuration for acquisition of measured signals from the sample assembly.

Fig. 1

(a) Configuration of signal detection with incident beam indicated by the red line: $D$: sample thickness; $D_1$: for monitoring $I_0$; $M$: mirror; $D_2$ to $D_4$: for measurement of $I_{\mathrm{Rd}},I_{\mathrm{Td}}$, and $I_{\mathrm{Tf}}$; ${\theta }_0$: angle of incident beam from $z$-axis; $d_R,d_T$, and $d_f$: distance between the origin and front center of $D_2$, $D_3$, and $D_4$; and ${\theta }_R,{\theta }_T$, and ${\theta }_f$: angle of sensor surface normal (blue dash lines) of D₂, D₃, and D₄ from the z-axis. (b) Work-flow chart of signal measurement in brown boxes and inverse calculation in green boxes.

2.2.

Signal Calculation by iMC

An in-house developed individual photon tracking MC (iMC) code was employed to calculate signals as $R_{dc},T_{dc}$ and $T_{fc}$ from given $\mathbf{P}$ for a phantom of the same shape as the sample inside a spacer between two glass slides by tracking N₀ photons, which imports the parameters of sample size and detection configuration as input data.^16,19,20 The code injects each of the N₀ 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 $N_{\mathrm{Rd}},N_{\mathrm{Td}}$, and $N_{\mathrm{Tf}}$ by the detector $D_2,D_3$, and $D_4$, respectively. The above process repeats until the total number of injected photons reaches $N_0$ and the calculated signals are given by $R_{\mathrm{dc}}=N_{\mathrm{Rd}}/N_0$, $T_{\mathrm{dc}}=N_{\mathrm{Rd}}/N_0$ and $T_{\mathrm{fc}}=N_{\mathrm{Tf}}/N_0$. Because of independence in trajectory among the ${\mathrm{N}}_0$ photons, the numbers of detected photons follow Poisson distributions.²¹ To estimate the variance in these photon numbers, one can draw a random number $q$ of Poisson distribution with $f(q,q_m)$ as probability mass function and $q_m$ as the mean. In the case of $N_{\mathrm{Rd}}$ calculated by an iMC simulation, $q_m$ equals to $N_0{R}_{\mathrm{dc}\infty }$ if distribution of $N_{\mathrm{Rd}}$ follows $f(N_{\mathrm{Rd}},q_m)$ and $R_{\mathrm{dc}\infty }$ yields the variance-free value of $R_{\mathrm{dc}}$ by tracking “infinite” number of photons. One thus can use $f(q,q_m)$ to quantify the effect of variance on calculated signals and optimize the objective function for variance reduction.

2.3.

PSO Algorithm

A stochastic algorithm based on PSO has been developed to guide iMC simulations and solve for ${\mathbf{P}}_s(\lambda )$ 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 $\mathbf{P}(\lambda )$, and optimization of this function reduces the difference to solve for ${\mathbf{P}}_s(\lambda )$ from an initial choice of $\mathbf{P}(\lambda )$. 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 $\mathrm{i}\in [1,\mathrm{I}]$ and $I$ as the number of threads or particles. Here, a particle $i$ symbolizes a thread of positions in the RT parameter space along which it moves from $\mathbf{P}(i,j)$ to $\mathbf{P}(i,j+1)$ as²²

3.1.

Effect of Variance in Calculated Signals on Solving ISPs

MPS requires solving one ISP for each value of $\lambda$ from the measured signals. The first step is to perform forward calculations of signals from given $\mathbf{P}=({\mu }_a,{\mu }_s,g)$ 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 ${\mathbf{S}}_{\mathbf{c}}(\mathbf{P})=(R_{\mathrm{dc}}(\mathbf{P})/R_d,$ $T_{\mathrm{dc}}(\mathbf{P})/T_d,T_{\mathrm{fc}}(\mathbf{P})/T_f)$. An inverse algorithm is to iterate iMC simulations toward the objective of making ${\mathbf{S}}_{\mathbf{c}}(\mathbf{P})$ as close to $\mathbf{1}=(\mathrm{1,1},1)$ as possible. A function of squared error sum $\delta (\mathbf{P})=|{\mathbf{S}}_{\mathbf{c}}(\mathbf{P})-\mathbf{1}{|}^2$ is often selected as the objective function to guide iteration and decide when stops. Ideally, $\delta ({\mathbf{P}}_s)$ vanishes for a perfect ISP solution given by ${\mathbf{P}}_s$. In practice, one deems an ISP as solved if $\delta (\mathbf{P})\le {\delta }_{th}$ with ${\delta }_{th}$ 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

3.2.

Comparison of Two Objective Functions for Variance Suppression

Based on above analysis, local averaging on $\delta (\mathbf{P})$ 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

Fig. 2

(a) The dependence of ${\rho }_{\delta }(\mathbf{P})$ and ${\rho }_{1/\delta }(\mathbf{P})$ on ${\delta }_f(\mathbf{P})$ with red line for ${\delta }_f(\mathbf{P})$ and blue line for ${\delta }_f(\mathbf{P}{)}^{-1}$. (b) The mean values of relative differences $\mathrm{\Delta }({\rho }_{\delta })$ between ${\rho }_{\delta }(\mathbf{P})$ and ${\delta }_f(\mathbf{P})$ and $\mathrm{\Delta }({\rho }_{1/\delta })$ between ${\rho }_{1/\delta }(\mathbf{P})$ and ${\delta }_f(\mathbf{P}{)}^{-1}$ against $q_m$ as the mean of Poisson distributions for calculated signals by tracking “infinite” photons.

To gain further insight on fluctuation in calculated signals, the relative difference Δ between ${\rho }_{\delta }(\mathbf{P})$ and ${\delta }_f(\mathbf{P})$ or ${\rho }_{1/\delta }(\mathbf{P})$ and ${\delta }_f(\mathbf{P}{)}^{-1}$ are calculated and averaged as

3.3.

Measurement of 20% Intralipid Sample

To validate the new approach, RT parameters of 20% intralipid (I141-100ML, Sigma-Aldrich) have been determined from three measured signals of $R_d,T_d$, and $T_f$ for $\lambda$ from 520 to 1000 nm in steps of 20 nm. Figure 3 presents the RT parameters of a sample with thickness $D=102\mu \mathrm{m}$ obtained by the PSO based algorithm with ${\rho }_{1/\delta }(\mathbf{P})$ 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 $n_r(\lambda )$ of the 20% intralipid were used to obtain interpolated values in iMC simulations at the wavelengths of measurement for the calculated signals and $n_r(\lambda )$ is presented in Fig. S2 in the Supplementary Material.²³ The same measurements were repeated on another sample of $D=81\mu \mathrm{m}$ and the inversely determined RT parameters agree well with those in Fig. 3 within the experimental errors.

Fig. 3

The wavelength dependence of RT parameters of one 20% intralipid sample of $D=102\mu \mathrm{m}$ in thickness. The error bars were inversely determined from different combinations of the mean and standard deviation values of the measured signals shown in Fig. S1 in the Supplementary Material on five selected wavelengths. The red line is a power law fitting of ${\mu }_s=C{\lambda }^{-2.340}$ with $C=1.626\times {10}^8$.

The search region ${\mathrm{\Gamma }}_s$ for PSO based inverse algorithm at each wavelength was set between $1.00\times {10}^{-4}$ and $2.00{\mathrm{mm}}^{-1}$ for ${\mu }_a$, 10.0 and $200{\mathrm{mm}}^{-1}$ for ${\mu }_s$, 0.10 and 1.00 for $g$. In additional to choosing $I=27$ for the particle number and $J=100$, we set $w_0=0.900$, $w_p=w_g=2.05$, and $\chi =0.730$ for the PSO parameters defined in Eq. (1) for optimized performance of inverse calculations. The initial positions of $\mathbf{P}(i,1)$ for all particles were randomly distributed in each of nine equal partitions of ${\mathrm{\Gamma }}_s$ to ensure globally optimized solutions of ISPs. We found that the convergence toward the final solution was not affected by the choice of $\mathbf{P}(i,1)$ under the above condition. After calculation of ${\rho }_{1/\delta }$ with $\mathbf{P}(i,1)$, iterations as defined in Eq. (1) were performed in ${\mathrm{\Gamma }}_s$. Specifically, ${\mathbf{P}}_p(i)$ was set as $\mathbf{P}(i,1)$ for particle $i$ and ${\mathbf{P}}_g$ as ${\mathbf{P}}_p(i^{\prime })$ that has the largest ${\rho }_{1/\delta }$ value among all particles for $j=1$ to obtain $\mathbf{v}(i,2)$ and $\mathbf{P}(i,2)$ by Eq. (1). The search then continued as j increases to 2, 3,… and stopped when either ${\rho }_{1/\delta }({\mathbf{P}}_g)$ exceeds ${\rho }_{th}$ for $j\le J$ or $j>J$ with ${\mathbf{P}}_g$ saved as ${\mathbf{P}}_s(\lambda )$ and output together with ${\rho }_{1/\delta }({\mathbf{P}}_g)$ and $\delta ({\mathbf{P}}_g)$. The threshold ${\rho }_{\mathrm{th}}$ corresponds to the value of ${\rho }_{1/\delta }$ in Eq. (5) with ${\delta }_{\mathrm{th}}$ set to 0.8% and $m_P$ to 7.

It took about eight iterations on average to solve an ISP of ${\mathbf{P}}_s(\lambda )$ 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 $N_0$ was gradually increased from $5\times {10}^5$ to $5\times {10}^6$ as $\delta$ became 1% or less. On average, over different values of RT parameters and $N_0$, each iMC simulation took about 0.2 s and solving ${\mathbf{P}}_s(\lambda )$ per wavelength took about 160 s by executing the code on one GPU board (Nvidia, GeForce RTX 2080 Ti). The results of ${\mathbf{P}}_s(\lambda )$ 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 ${\mu }_s$ and g agree well with those in Ref. 8 for the 20% intralipid sample while the range of ${\mu }_a$ is between those of Ref. 24 for 10% intralipid and Ref. 8. Since the 20% intralipid is high turbid with scattering albedo $a={\mu }_s/{\mu }_t>98\%$ over the measured wavelengths, the large differences in the values of ${\mu }_a$ can be attributed to the very large errors in determining ${\mu }_a$ with values smaller than ${\mu }_s$ by three orders of magnitude as previously discussed.⁸ We note further that the values of optical thickness $\tau$ and albedo $a$ of the intralipid sample reach the maximum values of 7.1 and 99.98%, respectively for $\lambda =520\mathrm{nm}$. Results of our study on samples of 20% intralipid with sufficiently large $D$ values showed that current MPS method with the inverse solver reported here can yield unstable solutions of ${\mathbf{P}}_s$ 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 $N_0$ much larger than those used for ISPs of $\tau \le 7.1$ for iMC simulations. A study is underway to further improve the PSO based inverse solver by taking into account the $\delta (\mathbf{P})$ distribution obtained from completed iterations, which is expected to enable accurate measurement of RT parameters for highly turbid samples with $\tau$ 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 $T_c$.