A linear gradient line source facilitates the use of diffusion models with high order approximation for efficient, accurate turbid sample optical properties recovery

: In this paper, we demonstrate that a scanning MEMS mirror can be employed to create a linear gradient line source that is equivalent to a planar source. This light source setup facilitates the use of diffusion models of increased orders of approximation having closed form solution, and thus enhance the efficiency and accuracy in sample optical properties recovery. In addition, compared with a regular planar light source, the linear gradient line source occupies much less source area and has an elevated measurement efficiency. We employed a δ -P 1 diffusion equation with a closed form solution and carried out a phantom study to understand the performance of this new method in determining the absorption and scattering properties of turbid samples. Moreover, our Monte Carlo simulation results indicated that this geometry had probing depths comparable to those of the conventional diffuse reflectance measurement geometry with a source-detector separation of 3 mm. We expect that this new source setup would facilitate the investigating of superficial volumes of turbid samples in the wavelength regions where tissue absorption coefficients are comparable to scattering coefficients.


Introduction
In recent years, Diffuse Reflectance Spectroscopy (DRS) techniques have been developed rapidly and widely used for many biomedical applications [1][2][3][4]. In a typical DRS configuration, light is injected into the sample under investigation and the diffuse reflectance is noninvasively detected at a distance from the source. The optical properties of the sample, including absorption (μ a ) and reduced scattering (μ s ') coefficients, can be deduced from the measured diffuse reflectance, and these parameters can in turn provide ample information on a variety of physiological processes [1,2,[5][6][7]. For most cases, DRS techniques rely on a photon transport model to relate the optical properties of samples to the diffuse reflectance. In general, the standard diffusion equation (SDE) derived from the radiative transport equation (RTE) is used to describe photon transport for its simplicity and computational efficiency [8,9]. However, the validity of SDE is held only when the scattering is much larger than the absorption and the photon travel length surpasses five to ten transport mean free paths [l t = 1/(μ a + μ s ')] [10][11][12]. Therefore, DRS techniques that work with the SDE are usually applied to study highly scattering thick tissue, such as breast or brain, in the 600 to 1000 nm wavelength region and at large source-detector separations [4,13,14].
On the other hand, to investigate superficial volumes of tissues with DRS, it is essential to limit the source-detector separation to be less than 5 mm, but SDE has been proven ineffective in this regime [10][11][12]. Several researchers have proposed diffusion models with increased orders of approximation to the RTE for the purpose of studying superficial tissues. For example, Hull and Foster employed the P 3 approximation to derive the Green's function of RTE and demonstrated successful reconstruction of sample optical properties at 0.43 mm source-detector separation [15]. By introducing the δ-Eddington approximation to the source function, Star et al. and Prahl demonstrated that this approach could improve the depiction of the photon energy transport near boundaries or sources more precisely than did a typical diffusion equation [16,17]. Hayakawa et al. used the δ-P 1 approximation (essentially the same with δ-Eddington approximation) based diffusion equation to successfully recover the optical properties of infinite geometry samples having μ s ' to μ a ratios ranging from 0.33 to 290 [18]. It should be noted that these advanced diffusion models do not necessarily lead to closed form solutions in the classical single source-detector pair(s), semi-infinite sample DRS measurement setup and thus have higher programming complexity and lower computational efficiency than the SDE. This has impeded their widespread use in superficial tissue studies.
Nevertheless, several groups have reported that the aforementioned advanced diffusion models had closed form solutions in the planar source illumination (PSI) geometry as illustrated in Fig. 1(a). For instance, Spott and Svaasand as well as Carp et al. provided detailed descriptions of the closed form solutions of P 3 and δ-P 1 diffusion equations in the PSI geometry [19,20]. However, because an infinitely wide planar light source has been approximated by expanding a point source to a large dimension, it therefore has relatively low source intensity resulting in low measurement efficiency and accuracy. Besides, the PSI geometry requires a large illumination area to satisfy the infinitely wide source geometry assumption in the theoretical models, which greatly reduces its feasibility for clinical measurements.
In this paper, we demonstrate that a scanning MEMS mirror can be utilized to create a linear gradient line source pattern that is equivalent to a regular planar light source. Compared with the regular PSI geometry, the linear gradient line source illumination (LGLSI) geometry occupies much less area as depicted in Fig. 1. In addition, as illustrated in Fig. 1(b) and 1(c), given a laser source of power P 0 is used for producing PSI and LGLSI, the irradiance of LGLSI is higher than that of PSI. Depending on the illumination radius and laser beam size, experimentally we could get about 5-20 times of detected signal strength in the LGLSI geometry than in the PSI geometry. This new source geometry can work in conjunction with an advanced diffusion model that has a closed form solution to recover the optical properties of samples. The theoretical background of an advanced diffusion model employed in this study will be provided in section 2. We carried out Monte Carlo simulations and phantom studies to understand the characteristics and evaluate the performance of this new method. The experimental setup will be described in section 3, and simulation and experiment results will be shown and discussed in section 4.

Theoretical background
A regular planar illumination geometry is shown in Fig. 1(a). We can further assume that this geometry has a circular shape with a radius r (r→∞) and is symmetric with respect to the origin O, as illustrated in Fig. 1(b). In such a planar illumination geometry with radial symmetry, the ratio of the number of photons entering the sample at r i to that at r j is r i /r j , where r i and r j represent the distances of points i and j to the origin. By utilizing this fact, we can design a line source whose intensity increases proportionally with the distance from the origin as shown in Fig. 1(c), and it is clear that such illumination geometry would be equivalent to a regular planar illumination. Several groups have proposed advanced diffusion models in the planar illumination geometry that have closed form solutions [19,20]. In this study, we adopted the δ-P 1 diffusion model developed by Carp et al. [20] and we will briefly describe its theoretical background in 2.1. The Monte Carlo model used as the benchmark method in this study for evaluating the performance of the diffusion models will be described in 2.2.

δ-P1 Diffusion model with a planar illumination source
Carp et al. considered the δ-Eddington phase function in the derivation of diffusion model from the radiative transport equation so that the fluence rate close to source and surface could be described more accurately [20]. The δ-Eddington phase function provides an additional degree of freedom well suited to accommodate collimated sources and highly forwardscattering media [20]. The modified model based on this phase function is usually called the δ-P 1 diffusion model. Carp et al. have provided the detail derivation of the δ-P 1 diffusion model used in this study. We will recapitulate the derivation in the following.
The δ-Eddington phase function was first introduced in 1976 by Joseph et al. [21]: where ŝ and ŝ′ represent the unit vectors of the light propagation before and after scattering, respectively. This phase function includes a Dirac delta function to model the collimated, unscattered component of light. The remainder of light, 1-f, is diffusely scattered according to a P 1 phase function and is composed of an isotopic term and an anisotropic term. The parameters g* and f can be obtained by approximating Eq. (1) to the Henyey-Greenstein phase function, which yields: where g is the anisotropy factor typically used in the P 1 approximation. Substituting the δ-Eddington phase function into the RTE, the RTE can be simplified to a δ-P 1 diffusion equation for a planar source in frequency domain: , ω is the modulation frequency, and c is the speed of light in the medium. By applying the boundary condition at the sample-air interface for Eq. (3), one can obtain the diffuse fluence rate Φ d . * where E 0 is the source strength, R s is the specular reflection coefficient, and μ eff = (3μ a μ t ') 1/2 . The coefficients α and β are defined as: and where h = 2μ t '/3, and A can be approximated as A(n) = −0.13755n 3 + 4.339n 2 -4.90366n + 1.6896 with n representing the refractive index of the sample. Finally, the diffuse reflectance at the surface can be determined using the following equation: The experimentally measured or Monte Carlo simulated reflectance can be fit to this δ-P 1 diffusion model using the MATLAB least square fitting method 'lsqcurvefit' (MathWorks, MA) to calculate absorption and reduced scattering coefficients of samples. The anisotropy factor g was set to 0.8 in the δ-P 1 diffusion model. The results will be demonstrated and discussed in section 4.

Monte Carlo simulation
Our Monte Carlo code was developed from the one proposed by Wang et al. [22]. In the Monte Carlo model, the source and the sample are defined in the cylindrical coordinate. In the code, the diameter of detector and the refractive index of sample were set to 400 μm and 1.33, respectively. The planar source radius was 30 mm and the radius of the sample was 10 6 cm to simulate infinite lateral sample dimension. Photon packets propagating in the media followed the Henyey-Greenstein phase function with the anisotropy factor g set to 0.8. All photon packets were launched with unity weight. Each simulation ran continuously until a total number of 10,000 photon packets arrived at the detector to ensure the standard deviations of amplitude and phase were within 1.0% and 0.5°, respectively, in the frequency range from 0 to 1500 MHz. In addition, the probing depth of each simulation was determined by calculating the weighted average of the maximum penetration depth of each detected photon packet and was defined as following: Let W i and (z max ) i be the final weight and the maximum penetration depth of the i-th detected photon packet, respectively, and n be the number of photon packets detected. The probing depth of a simulation was defined as

Frequency domain photon migration system and linear gradient line source illumination geometry
A frequency domain photon migration (FDPM) instrument was employed to measure the amplitude and phase of photon density waves from 10 to 200 MHz. The setup of the FDPM instrument is depicted in Fig. 2. In the FDPM system, a network analyzer (N5230C, Agilent Technologies, CA) was used to generate 14 dBm RF signal at frequencies spanning from 10 to 200 MHz. A 780 nm laser diode with output power of 18 mW was mounted on a TEcooled laser diode mount (TCLDM9, Thorlabs, NJ) biased and temperature controlled by a laser diode controller (LDC-3908, Newport, CA). The RF signal from the network analyzer provided current modulation to the laser diode to produce intensity-modulated light. The light from the laser diode was collimated using an aspheric lens (C330TME-B, Thorlabs, NJ) and was then deflected by an aluminum coated MEMS (Micro Electro Mechanical Systems) mirror (Mirrorcle Technologies, CA) which had 1.2 mm diameter and 2.6 kHz resonant frequency. The light beam could be steered to various locations on a phantom by the MEMS mirror, and it could be manipulated by the MEMS mirror to generate the linear gradient pattern as shown in Fig. 1(c). We designed the stay time of the light beam at 50 evenly distributed points within the 10 mm scanning radius to be proportional to the distance to the origin. The MEMS mirror had a sample rate of 50000 s −1 and it took around 0.02 s to complete one round of pattern scanning. In one round of pattern scanning, the projected energy on the line source was in the range from 0.36 to 18 μJ. The detected light was converted to a RF signal by the avalanche photodiode detector (APD) (Module C5658, Hamamatsu photonics K.K., Japan) and the RF signal was sent to the network analyzer to determine the amplitude demodulation and phase delay introduced by the sample under investigation. The amplitude data represent a ratio between the sample and reference signal amplitudes, while the phase data represent the phase difference between the sample and reference signals. This system was controlled by LabVIEW 8.5 (National Instruments, TX), which allowed users to save the measurement data in the laptop for further data processing.

Liquid phantom
We made three homogeneous liquid phantoms with known optical properties in this study. These phantoms were fabricated by mixing proper amount of scattering agent Lipofundin MCT/LCT 20% (B. Braun Melsungen AG, Germany), water-soluble dye Nigrosin (MP Biomedicals, Inc., Germany), and de-ionized water. Table 1 shows the benchmark absorption and reduced scattering coefficients of these three phantoms at the wavelength of 780 nm.

Theoretical performance of the δ-P 1 diffusion model
Carp et al. demonstrated that by taking the collimated photon propagation into account the δ-P 1 diffusion model could accurately predict light distribution even when the μ s ' to μ a ratio is low [20]. To evaluate the performance of the δ-P 1 diffusion model, Monte Carlo simulations were conducted to obtain the benchmark reflectance of two samples with very different μ s ' to μ a ratios. In addition, the standard diffusion model in the PSI geometry derived by Pham et al. [23] was also included for comparison. The optical properties of high and low μ s '/μ a samples were set to μ a = 0.01 mm −1 , μ s ' = 1.67 mm −1 and μ a = 1.00 mm −1 , μ s ' = 1.67 mm −1 , respectively. Amplitude and phase generated using the Monte Carlo simulations, δ-P 1 diffusion model, and the standard diffusion model are shown in Fig. 3(a) and 3(b). It can be seen in Fig. 3 that for the sample with μ s '/μ a = 167.00, the amplitude deviations of the two diffusion models from the MC data are similar (both less 3.0%), while the phase deviations of the δ-P 1 diffusion model and the standard diffusion model are −1.3° (−5.7%) and −2.1° (−9.3%), respectively, at the modulation frequency of 1500 MHz. The performance of the δ-P 1 diffusion models is slightly better than the SDE for the high μ s '/μ a sample. The advantage of the δ-P 1 diffusion model becomes obvious when the sample μ s ' to μ a ratio reduces to 1.67. The deviations of amplitude and phase of the δ-P 1 diffusion model and the standard diffusion model from the Monte Carlo data are 21.1%, 0.2° (8.5%), and 59.0%, −0.64° (−27.1%), respectively, at 1500MHz. It can be concluded from these simulation results that the δ-P 1 diffusion model is advantageous over the SDE for both high and low μ s ' to μ a ratio samples.

Probing depth of the planar source illumination geometry
Since in the PSI geometry the launch position of photon packets is not definite, it is not easy to intuitively infer the probing depth of such setup. To resolve this problem, we calculated the maximum probing depths of the PSI geometry at various sample optical properties using Monte Carlo simulations. In the simulations, we selected a wide spectrum of sample optical properties based on tissue optical properties previously reported by several groups in the wavelength range from 500 to 1000 nm [14,[24][25][26][27]. The probing depths of 1071 samples with different optical properties were determined using Monte Carlo simulation to generate the plot shown in Fig. 4(a). It can be seen in Fig. 4(a) that the probing depth of the PSI geometry has a large variation across the spectrum of optical property investigated; it changes from 6.0 mm for low μ s ' and μ a samples to less than 0.5 mm for high μ s ' and μ a samples. This phenomenon is reasonable since for high μ a samples only photons launched at distances close to the detector, which have shallow interrogation depths, would have a good chance to get to the detector. On the other hand, for low μ a samples, photons launched at distances far from the origin with deep interrogation depths still have some chances to reach the detector. Thus, the probing depth of the PSI geometry varies greatly with sample optical properties. It is found that the human skin optical properties are generally greater than 1 mm −1 and 0.04 mm −1 for μ s ' and μ a , respectively, in the wavelength range from 500 to 1000 nm [24,27]; thus, for most skin applications, it can be seen from Fig. 4(a) that the maximum probing depth of the PSI geometry is typically less than 2.0 mm. We further carried out another series of Monte Carlo simulations to calculate the maximum probing depths of a classical DRS geometry with 3 mm source-detector separation. The results determined from 1071 Monte Carlo simulations are illustrated in Fig. 4(b). We can observe in Fig. 4(b) that the probing depth varies from 3.5 mm for low μ s ', low μ a samples to 0.5 mm for high μ s ', high μ a samples. Moreover, it can be estimated that the probing depths of the classical 3 mm source-detector separation geometry are in the range from 0.5 to 2 mm as this geometry is employed to study human skin. We have also run Monte Carlo simulations for a classical DRS geometry with 2 mm source-detector separation and found the probing depths were in the range from 0.5 to 1.5 mm for skin measurements. Therefore, it can be reasonably inferred that from the aspect of skin interrogation depth, the PSI geometry is comparable to the classical DRS geometry with 3 mm source-detector separation.

Sample optical properties recovery
Monte Carlo method was employed to calculate the benchmark frequency domain diffuse reflectance of two semi-infinite samples illuminated with the planar source with modulation frequency ranging from 10 to 200 MHz. The sample optical properties were μ a = 0.01 mm −1 , μ s ' = 1.67 mm −1 (μ s '/μ a = 167.00) and μ a = 1.00 mm −1 , μ s ' = 1.67 mm −1 (μ s '/μ a = 1.67), respectively. Normal distribution random noise of 0.3% in amplitude and 0.08° in phase delay were added to the Monte Carlo data to simulate measurement noise. The data were then fit to the δ-P 1 diffusion model to extract the sample optical properties. The recovery errors are listed in Table 2. The recovery errors are all within 3% for the two samples having μ s ' to μ a ratios different by two orders of magnitude. This simulation result suggests that the δ-P 1 diffusion model is capable of precisely describing photon transport for a wide μ s '/μ a range, which agrees with those reported earlier by Spott and Svaasand, and Hayakawa et al. [18,19]  We also conducted phantom measurements using our FDPM system illustrated in Fig. 2. Due to the limited performance in producing strong intensity modulated light and the limited modulation frequency range of our laser diode, we could only measure samples with high μ s ' to μ a ratio at modulation frequencies from 10 to 200 MHz. We measured two liquid phantoms LP_1 and LP_2 whose μ s ' to μ a ratios were 167 and 136, respectively. We employed LP_3 to calculate the system response so that the system response calibration for other phantom measurements can be performed. A typical fitting results for LP_1 are illustrated in Fig. 5. The average optical property recovery errors of three measurements for each phantom are shown in Table 2. It can be seen that the recovery errors are all within 7%. Pham et al. and our group have previously reported that employing a SDE for sample optical property recovery in the PSI geometry could lead to around ± 9% and 50% recovery errors, respectively [23,28]. Compared to the results reported earlier, our results show that the LGLSI geometry combined with an advanced diffusion model is advantageous in accurate sample optical property recovery. Besides, it is reasonable to expect that this method could be applied to measure samples with comparable μ s ' and μ a , given a laser diode with better modulation capability is employed in the FDPM system.

Illumination radius of linear gradient line source
Our Monte Carlo model was capable of recording the incident positions of all photon packets reached at the detector locating at the origin in the PSI geometry. By analyzing the incident positions of all detected photon packets, a plot of detected photon packet number versus the incident position as shown in Fig. 6(a) can be generated. In addition, Fig. 6(a) illustrates the percentile cumulative distribution function (CDF) of the detected photon packet number. The sample optical properties used in the Monte Carlo simulation for generating Fig. 6(a) were set the same as those of the liquid phantom LP_1. It can be observed in Fig. 6(a) that the incident positions of about 90% of total detected photon packets are within 10 mm from the origin. We thoroughly looked up our simulation data and found that for measuring samples with low μ s ' and μ a in the PSI geometry, a source radius as large as 30 mm could be required. For example, the detected photon packet number versus the incident position for the sample with lowest μ s ' and μ a in our simulation data (μ a = 0.001/mm and μ' s = 0.5/mm) is illustrated in Fig.  6(b). It can be seen in Fig. 6(b), there are still some photon packets launched at around 30 mm could arrive at the detector. The CDF shown in Fig. 6(b) indicates that 90% of detected photon packets were launched within 15 mm from the detector. To ensure that the incident positions of most of detectable photon packets could be recorded at a proper simulation efficiency, we had used 30 mm illumination radius for all PSI Monte Carlo simulations demonstrated so far. Since in the derivation of the diffusion model the source radius was assumed to be infinite, the finite illumination radius of the PSI geometry in the measurements or Monte Carlo simulations would lead to inevitable sample optical property recovery errors. To understand the influence of the illumination radius on the recovered sample optical properties, we gradually reduced the illumination radius from 30 to 3 mm in the Monte Carlo simulations for the sample with μ a = 0.01 mm −1 , μ s ' = 1.67 mm −1 . The simulated data were then fit to the δ-P 1 diffusion model to determine the sample optical properties and the results are shown in Table 3. The corresponding cumulative distribution values at the illumination radii are also listed in Table 3. It can be seen in Table 3 that the μ a recovery error increases as the source illumination radius decreases; however, the μ s ' recovery error does not follow a similar trend. We speculate this phenomenon was due to the fact that the absorption mean free path (~100 mm) and the reduced scattering mean free path (~0.6 mm) are very different for this sample. As photons' incident positions are far from the detector in the LGSLI geometry, they could have relatively good chances of encountering absorption events; thus as shown in Table 3, when the illumination radius is larger than 6.5 mm, the μ a recovery errors are smaller than 10%. As the illumination radius reduces to less than 6.5 mm, photons may have relatively poor chances of encountering absorption events, and thus the μ a recovery errors are elevated.
In contrast, because the reduced scattering mean free path is very short, photons have good chances of encountering scattering events at all illumination radii listed in Table 3; therefore, the μ s ' recovery errors are not greatly influenced by the illumination radius. Furthermore, for this sample, the value of cumulative distribution has to be greater than 85% to keep the optical property recovery error less than 10%. We have checked many cases and found that as long as the illumination radius was set at a value ensuring the CDF reached 85%, the requirement of recovery errors of optical properties within 10% could be satisfied in general. Thus we defined the minimum illumination radius of the PSI geometry at which the 85% CDF was reached for a certain sample. By using the definition of minimum illumination radius stated above, we determined the minimum illumination radii for all sample optical property combinations investigated in this study using Monte Carlo simulations. The results are illustrated in Fig. 7. The data shown in Fig. 7 suggests that 10-12 mm illumination radius of the PSI geometry would be sufficient for most skin applications. The LGLSI radius of the measurements discussed in section 4.3 was 10 mm. We adjusted the LGLSI radius to 5 mm in the experiment and found that the recovery errors of μ a and μ s ' were −8.74%, −5.19% and −6.56%, −9.92%, for LP_2 and LP_1, respectively. Compared to the values listed in Table 2, the recovery error increases as the LGLSI radius decreases. This trend obtained from experiments agrees well with that of Monte Carlo simulation results shown in Table 3. 18 Fig. 7. The minimum scanning radius in mm of PSI geometry at various sample optical properties.

Conclusion
We have proposed a novel line source setup in this study. The line source has a linear gradient intensity variation which is conceptually equivalent to a planar source. Compared with a regular planar source, this LGLSI geometry increases the measurement efficiency by at least 5 to 20 times. In addition, because the LGLSI geometry has a more localized probing volume compared to the PSI geometry, it would be less affected by the tissue surface flatness and be more sensitive to the presence of an abnormal inhomogeneity. Besides, the detection fiber does not block a part of illumination in the LGLSI geometry as it does in the PSI geometry; therefore the repeatability in the LGLSI geometry is better than that in the PSI counterpart. We experimentally demonstrated that the LGLSI geometry can work with advanced diffusion models such as the δ-P 1 diffusion model to accurately recover the sample optical properties. The maximum probing depth of this new LGLSI geometry was found to be less than 2 mm for measuring human skin, similar to that of a classical DRS setup with 3 mm source-detector separation. In addition, from Monte Carlo simulations, the minimum LGLSI radius for most skin applications was estimated to be around 10 mm. With shallow interrogation depth, small source illumination area, and the capability of employing advanced diffusion models with closed form solutions, it is expected that the LGLSI geometry can be applied to study various superficial tissues in a wide wavelength range covering Visible to Near IR.