Efficient construction of robust artificial neural networks for accurate determination of superficial sample optical properties

: In general, diffuse reflectance spectroscopy (DRS) systems work with photon diffusion models to determine the absorption coefficient μ a and reduced scattering coefficient μ s ' of turbid samples. However, in some DRS measurement scenarios, such as using short source-detector separations to investigate superficial tissues with comparable μ a and μ s ' , photon diffusion models might be invalid or might not have analytical solutions. In this study, a systematic workflow of constructing a rapid, accurate photon transport model that is valid at short source-detector separations (SDSs) and at a wide range of sample albedo is revealed. To create such a model, we first employed a GPU (Graphic Processing Unit) based Monte Carlo model to calculate the reflectance at various sample optical property combinations and established a database at high speed. The database was then utilized to train an artificial neural network (ANN) for determining the sample absorption and reduced scattering coefficients from the reflectance measured at several SDSs without applying spectral constraints. The robustness of the produced ANN model was rigorously validated. We evaluated the performance of a successfully trained ANN using tissue simulating phantoms. We also determined the 500-1000 nm absorption and reduced scattering spectra of in-vivo skin using our ANN model and found that the values agree well with those reported in several independent studies.


Introduction
Diffuse reflectance spectroscopy (DRS) is an optical technique which offers a means to noninvasively and quantitatively characterize the macroscopic structural and physiological parameters of biological tissues, such as brain, breast, and skin [1-3]. In a typical DRS setup, visible to near-infrared light is injected into the tissue and one or several detectors are used to collect the diffuse reflectance. With a proper photon transport model such as the standard diffusion equation (SDE), derived from the radiative transport equation with P 1 diffusion approximation, or the Monte Carlo method, the optical absorption and scattering properties of samples could be recovered from the measured reflectance [3,4]. To interrogate tissue volumes lying between the tissue surface and 1 to 2 mm deep using DRS, measurements with source-detector separation shorter than 3mm are necessary [5,6]. In this regime, the SDE generally cannot be reliably used to describe photon propagation due to a series of approximations made in deriving the SDE. For example, SDE considers only the first moment of the phase function, but reflectance measured at short source-detector separations are sensitive to the phase function [7]. Besides, diffusion approximation requires that the medium has a high-albedo (predominant scattering, μ s ' >> μ a ) [8], but this requirement is often not satisfied for biological tissues at wavelengths shorter than approximately 600 nm and at wavelengths longer than 1000 nm [9]. Thus, SDE is not suitable for recovering sample optical properties at short source-detector separations [6,7]. Advanced diffusion models with extended order of approximation were proposed by several researchers, and they were verified to be more accurate than the SDE near boundaries and sources. For instance, the P3 approximation employed the first four Legendre polynomials (P0; P1; P2; P3) to precisely model the radiance [10]. Although the performance of the P3 approximation method was demonstrated to be better than that of SDE, analytical solutions for a pencil beam light source in the semi-infinite or layered sample geometries have not been found, which impede its use in practical applications.
On the other hand, the Monte Carlo method offers a flexible yet rigorous approach for simulating photon transport in turbid media [11]. Owing to its accuracy and its ability to handle heterogeneous optical properties, the Monte Carlo method is currently widely accepted as the gold standard approach for photon migration modelling. However, the need for intensive computation resources limits its applicability for real-time calculation and routine clinical applications. To enhance the computational speed, Monte Carlo simulations carried out in a parallel configuration by employing graphics processing units (GPUs) has been developed [12,13]. Although GPU-based Monte Carlo method has greatly reduced the computation time by almost three orders of magnitude comparing to the single CPU based Monte Carlo method [13], it is still not a practical forward model for recovering the optical properties of turbid media in real time.
Our goal of this study was to develop an efficient and robust model for quantitatively determining the absorption and scattering coefficients of superficial tissues in a broad wavelength range without applying spectral constraints. To this end, an artificial neural network (ANN) was generated to relate the sample absorption and scattering coefficients and the reflectance at a precision level comparable to that of a Monte Carlo model. The ANN technique is one kind of learning algorithms that can be trained to find the connections between input and output parameters. The ampler and more exact the information provided to the ANN training function, the more precise the connection between input and output parameters can be established. Researchers have demonstrated that well-trained ANNs could accurately map the input optical properties to the diffuse reflectance of media [14][15][16][17]. Farrell et al. first showed that using an trained ANN as the forward model in the iterative fitting process, the optical properties of samples could be derived from spatially resolved diffuse reflectance measured at 8 source-detector separations (SDSs) ranging from 3 to 15 mm [14]. However, in Farrell's study, SDE was used for the ANN training due to training time consideration which led to a compromised ANN accuracy; besides, the source-detector separations employed were not suitable for superficial tissue investigations. Later, Kienle et al. used an ANN that was trained with 120 Monte Carlo simulations to determine the optical properties from reflectance measurements at 9 or 11 SDSs ranging from 2 to 12 mm [15]. They reported that it was necessary to simultaneously employ two ANNs to determine the optical properties of sample in the range 0.2 < μ s ' < 2.5 mm −1 and 0.001 < μ a < 0.5 mm −1 due to the wide range of light attenuation introduced by turbid samples; besides, consistent results could not be obtained when using SDSs shorter than 2mm. Similarly, Jäger et al. utilized 6000 Monte Carlo simulated spatially resolved diffuse reflectance to train 23 ANNs for determining the optical properties of turbid media (11 SDSs ranging from 1 to 23 mm) [16]. In their setup, the measured reflectance were first classified according to their shapes and then sent to a certain ANN to recover the optical properties, and the recovered optical property errors could be within ± 20%. Wang et al. demonstrated that a two-stage ANN could be trained with data generated from relative computational efficient condensed Monte Carlo simulations (6 SDSs ranging from 0.25 to 2.25 mm) [17]. Their experimental measurements showed that the deviation of recovered optical properties were up to 44%.
From a thorough literature review, we have found that although the ANN method could potentially be utilized for rapid, accurate tissue optical properties recovery, obtaining a proper ANN model for a specific application is not a simple task. The difficulties encountered when constructing a useful ANN that were reported in the literature can be boiled down to three major aspects. 1. The number and the range of SDSs required varied in different studies. In some studies, a wide range of SDSs was employed which caused the magnitude difference between the largest and the smallest diffuse reflectance could be higher than 3 orders. This fact would lead to either an unsuccessfully ANN training or an ANN with poor performance. 2. The size of training set should be optimized to lead to the best ANN training result. Since the Monte Carlo method has been recognized as the best tool for producing the ANN training database, the size of the training database is usually constrained by the heavy computational loading required for carrying out Monte Carlo simulations. 3. To understand the theoretical performance of an ANN after the training process is finished, a careful validation process has to be performed. Some investigators did not indicate whether the validation process was performed, while others proposed various validation methods. There is no consensus regarding how an ANN validation should be carried out.
In this study, we overcame the three difficulties mentioned above and proposed a workflow for constructing an ANN for rapid, accurate determination of the optical properties of superficial samples from the spatially resolved reflectance without applying spectral constraints. To increase the efficiency of the ANN construction process, the GPU-based Monte Carlo model was utilized in this study. By using the GPU-based Monte Carlo model, a training database of desirable sizes, which should be optimized according to the range of sample optical properties of interest, could be generated in a reasonable time frame. Therefore, it was possible to iteratively adjust the size of training database and the number and the range of SDSs until an optimization was achieved. After an ANN was built, a validation process was carried out to ensure the validity of the ANN in the optical property range of interest. One of our successfully trained ANNs only needs to take reflectance measured at two SDSs as the input parameters for rigorously recovering the optical properties of tissue simulating phantoms and in-vivo skin in the wavelength range from 500 to 1000 nm. To our knowledge, this is the first study that reveals the method for efficiently constructing an ANN model that can rigorously recover the superficial sample optical properties by using reflectance collected at only two SDSs. It will be seen in our phantom measurement results that one of our successfully trained ANN model which took reflectance measured at SDSs of 1 and 2 mm as inputs can determine the phantom optical properties with less than 1% recovery error. The performance of this ANN model is much better than the ones presented by Wang et al. using a similar SDS range [17]. The experimental results and the analysis of the difference between the performance of our ANN and a SDE will also be demonstrated.

GPU-MCML and ANN
The basic structure of Monte Carlo modeling of light transport in multi-layered tissues (MCML) was initially developed by Wang et al. [11]. To enhance the computation speed, Alerstam et al. proposed a highly optimized implementation utilizing the Compute Unified Device Architecture (CUDA) [12]. It is the GPU version of the widely cited MCML code and was referred to as GPU-MCML [13]. In contrast to the MCML code where only one photon packet is simulated at a time using a CPU, the GPU-MCML code utilizes GPUs that have parallel throughput architectures to simultaneously simulate many photon packets in parallel using many scalar processors [13]. In general, GPU-MCML could be faster than MCML by near three orders of magnitude [13].To construct a database for the ANN training, we first used the GPU-MCML code to calculate the diffuse reflectance for various optical properties and various SDSs. Because our interest was in skin applications, the range of optical properties was carefully defined based on the values reported by literatures. We have found that absorption coefficients and the reduced scattering coefficients of skin are in the range between 0.07 and 0.5 mm −1 , and between 0.5 and 3.5 mm −1 , respectively, from 500 to 1000 nm [3, 18,19]. We slightly expanded the range described above and defined the ranges of μ a and μ s ' of interest in this study to be from 0.015 mm −1 to 0.51mm −1 and from 0.2 mm −1 to 5 mm −1 , respectively. The number of μ a and μ s ' combination was varied from 3000 to 10000 depending on the ANN training results. In each GPU-MCML simulation, one hundred millions photon packets were launched, and at least ten hundred photon packets could be received by the detector. Other simulation parameters are listed below: the refractive index of sample was 1.33 for liquid phantom studies and 1.43 for skin studies, the anisotropy factor was 0.9, and the phase function was the Henyey-Greenstein phase function. The refractive indices of fibers were not taken into account in the simulations.
The time needed for constructing the database depends on the size of database, and it took less than 20 hours to calculate the reflectance of 6100 combinations of μ a and μ s ' on a computer equipped with a NVIDIA GTX 780 Ti graphic card. Once the database was constructed and carefully examined, it could be used for the ANN training. The concept of the ANN training is illustrated in Fig. 1. In this example, an ANN model is to be built to relate 6100 sets of μ a and μ s ' (output) to the corresponding reflectance sets (input); each reflectance set consists of two values collected at two SDSs. A successful ANN training process could take a few hours depending on the complexity of the ANN structure. In this study, we employed the Matlab (Mathworks, MA) training function "trainbr" based on Bayesian regularization to update the weight and bias values according to Levenberg-Marquardt optimization. We chose the Log-sigmoid transfer function "logsig" for the hidden layer to train these network models. Depending on the size of training database, we used 30-80 neurons in the ANN training for different source-detector separation groups. The numbers of hidden layers and neurons had to be properly adjusted to optimize the ANN performance. Consequently, an ANN model with the 2-in-2-out configuration could be obtained, and it could take the reflectance measured at 2 SDSs as the input and determine μ a and μ s ' of the sample under investigation.

Liquid phantom
We fabricated homogeneous liquid phantoms by using Lipofundin 20% (B. Braun Melsungen AG, Germany) as the scatterer, the water-soluble dye Nigrosin (MP Biomedical Inc., Germany) as the absorber (0.18 g Nigrosin powder in 1000 ml de-ionized water), and de-ionized water as the solvent. Lipofundin is fat emulsion that produced from safflower oil, egg phospholipids and glycerin. It contains emulsified fat particles and frequently used in the research of light propagation in turbid media. Nigrosin is a powder mixing of synthetic black dyes (CI 50415, Solvent black 5), and it has a prominent absorption peak at 576 nm. The recipe of the two liquid phantoms used in this study is listed in Table 1. These phantoms were designed to mimic the optical property spectra of skin in the 500 to 1000 nm region and thus had relatively high absorption (~0.4 mm −1 ) around 550 nm and moderate reduced scattering.

Inverse adding doubling (IAD) method
The inverse adding doubling method has been shown to be able to accurately recover μ a and μ s ' from the reflectance and transmittance measurements of turbid slab samples [20]. We employed this method in this study to define the benchmark absorption and reduced scattering coefficients of the two liquid phantoms. The system configuration of our integrating sphere measurement system is described in the following. Our light source was a supercontinuum laser (SuperK COMPACT, NKT Photonics) which provided a wide output spectrum covering the 450-2400 nm region. An adjustable aperture was used to limit the diameter of light beam to be within 2.5mm. A quartz cuvette was used to contain the liquid phantom and a cuvette holder was attached to the integrating sphere (3P-GPS-053-SL, Labsphere, NH) ensuring that the cuvette was right against the integrating sphere entrance. A frame reducer (PFR-FM/M-250-100-SF, Labsphere, NH) was used to reduce the entrance port from 2.5" diameter to 1". An optical fiber was connected to the exit port to collect a part of light from the integrating sphere and deliver the light to a spectrometer (Acton PIXIS 256 and SP2300, Princeton Instruments, NJ) to record 500 to 1000 nm spectra. Measured sample reflectance and transmittance were processed with the inverse adding doubling algorithm to derive the sample optical properties.

DRS system
The schematic of our DRS system is depicted in Fig. 2. In the system, a spectrometer equipped with a back-thinned CCD (QE65000, Ocean Optics, FL) was used to collect the reflectance from the detector fiber in the wavelength range from 500 to 1000 nm. To achieve shallow interrogation depths for skin applications, an optical fiber probe consisted of one detector fiber and seven source fibers with a 1 mm core to core separation between any two adjacent fibers was fabricated. All optical fibers employed in the probe are multimode fibers with 400-μm core diameter and 0.2 numerical aperture. Among the fibers, we found that SDSs in the range from 1 to 3 mm were favorable. In this range of SDSs, we formed three groups of different SDS combinations: 1 and 2 mm, 2 and 3 mm, and 1, 2, and 3 mm. The optical fiber connecting to a supercontinuum source was bridged to one of the three source optical fibers of the probe using a 1x4 optical switch (Piezosystem Jena, Germany). The spectrometer and the optical switch were connected to a laptop computer and were coordinated and controlled by a graphic user interface developed based on Labview (National Instruments, TX). When we employed the DRS system for phantom and skin studies, five measurements for each sample were carried out so that the standard deviation of measurements could be calculated.

In-vivo skin measurement
To understand the differences between the performance of the ANN and SDE in recovering skin optical properties, we measured the diffuse reflectance spectra from the inner arm and the outer thigh skin of a subject with our DRS system. The protocol was approved by the Institutional Review Board (No. ER-100-332) and the written informed consent was obtained from the subject prior to the measurements.

ANN training
The goal of this study was to obtain an ANN that could accurately compute accurate sample optical properties based on the input reflectance collected at several SDSs. To understand how to efficiently obtain a useful ANN, we elaborate the key procedures we adopted in the following examples. In the first example, we were aiming for producing an ANN that took reflectance collected at SDSs of 2 and 3 mm as the input parameters. The first step was to generate a training database that contained the information of the reflectance at the two SDSs for various optical property combinations. The size of the database would affect the ANN training result, and in the optical property range of interest in this study, we typically divided μ a and μ s ' ranges into 100 and 61 equal steps, respectively, to construct a database containing 6100 reflectance sets. As stated earlier, the required size of database required depends on the range of the sample optical property covered. We tested several database sizes and found that using 6100 database size could minimize the network training time without sacrificing the performance of the trained network in the optical property range of interest in this study.
Next, it was important to thoroughly check the uniqueness of each reflectance set (2 spectra at SDSs of 2 and 3 mm) in the database. This was to ensure that there was no two or more different optical property combinations that had the same reflectance values at SDSs of 2 and 3 mm; had this condition be true, the ANN training would fail under this one-in-multi-out situation. Experimentally, the condition of uniqueness of a reflectance set is related to the measurement uncertainty of the DRS system. We found that the major source of the measurement uncertainty in our system came from the instability of the light source, and other system components had negligible contribution. We monitored the long term stability of our light source and found that the fluctuation of light intensity was within 0.5% in a 2-hour time frame. Thus we defined that any two sets of reflectance that had differences within 0.5% at both 2 and 3 mm SDSs were not distinguishable and therefore not unique. After exhaustively checking the uniqueness of all reflectance sets in the database, one map as shown in Fig. 3 can be plotted. To plot Fig. 3, the number of reflectance sets deviating from a certain reflectance set within 0.5% at both 2 and 3 mm SDSs for a certain optical property combination was determined and this action was repeated until all 6100 reflectance sets in the database had been processed. Fig. 3. The number of indistinguishable reflectance sets at various μ a and μ s ' combinations. Each point in the plot indicates the number of reflectance sets in the database that have difference within 0.5% from a certain reflectance set corresponding to a certain μ a and μ s ' combination. Each set of reflectance consists of reflectance collected at source to detector separations of 2 and 3 mm.
It can be seen in Fig. 3 that the largest number of indistinguishable reflectance sets are 4. If the experimental indistinguishable reflectance sets correspond to optical properties that are far away from each other, the ANN training would fail; because this one input (one reflectance set), multiple outputs (multiple sets of optical properties) condition does not meet our preset ANN training configuration. However, if the indistinguishable reflectance sets correspond to optical properties that are close enough, one still could obtain a useful ANN (one input, one output). To understand this point, we calculated the largest optical property variation between the indistinguishable reflectance sets, and the results are depicted in Fig. 4. For a certain sample optical property combination, its corresponding reflectance set could be indistinguishable from those of other sample optical property combinations; each data point in Fig. 4 indicates the maximal variation between these sample optical property combinations. We can observe in Fig.  4(a) that the maximal variation of μ a is up to nearly 200% in the low μ a , high μ s ' region. On the other hand, the maximal variation of μ s ' is up to 70% as shown in Fig. 4(b). The results shown in Fig. 4 imply that it is very unlikely to produce an useful ANN based on a database in which two sets of optical properties that have 200% difference in μ a or 70% difference in μ s ' could correspond to two sets of reflectance that are not distinguishable. This problem could be alleviated by reducing the optical property range of the database. From the skin optical properties reported by various groups, it can be found that μ a lower than 0.1 mm −1 and μ s ' higher than 3.5 mm −1 do not happen simultaneously for skin in the wavelength range from 500 to 1000 nm. By taking off the data corresponding to μ a lower than 0.1 mm −1 and μ s ' higher than 3.5 mm −1 from the database, the maximal optical property variation between the indistinguishable reflectance sets are plotted in Fig. 5. It can be seen in Fig. 5 that the maximal optical property variation are within 16% and 13% for μ a and μ s ', respectively. Although the maximal variations shown in Fig. 5 are much smaller than those shown in Fig. 4, they are still higher than 10%, an empirical threshold defined by our group. We hope that the recovered optical properties could have errors less than 10% given a 0.5% measurement uncertainty. Therefore, based on this criterion, it can be concluded that this database is not suitable for further use in the ANN training. To proceed, we followed the same procedures described above to investigate whether the diffuse reflectance collected at other SDSs, say 1 and 2 mm, could be employed for constructing a useful database. After a database containing 6100 reflectance sets was generated with SDSs of 1 and 2 mm, an exhaustively examination was carried out to check the uniqueness of all reflectance sets in the database. Similarly, any two sets of reflectance that have differences within 0.5% at both 1 and 2 mm SDSs can be defined as not distinguishable. The largest number of indistinguishable reflectance sets was 3 in the optical property range of interest. The maximal percent variation between the optical property combinations that corresponded to indistinguishable reflectance sets was also calculated, and the results are illustrated in Fig. 6. It can be observed in Fig. 6 that the maximal variations are less than 2% for both μ a and μ s ', which meet our previously stated criterion. This database was considered as a useful one and then sent to the Matlab Bayesian regulation backpropagation network training function 'trainbr' to generate an ANN. A typical training process would take a couple of hours. The performance of the resultant ANN was first evaluated using the training database. All reflectance sets in the database were taken as the input parameters of the ANN to compute the optical properties which were then compared to the benchmark optical properties to determine the percent errors. It was found that the errors were less than 0.72% and 0.03% for the recovered μ a and μ s ', respectively, as demonstrated in Fig. 7. Next, 1000 optical property combinations, excluding the optical property combinations employed in the training database, were randomly chosen in the optical property range of interest, and the reflectance sets were calculated using GPU-MCML simulations for these 1000 optical property combinations to form a validation database. By sending the 1000 reflectance sets to the ANN, 1000 optical property sets were recovered and their percent errors were determined. The percent errors of the recovered optical properties are depicted in Fig. 8. It can be seen in Fig. 8 that the recovery error for this ANN was less than 2.7% for μ a and 0.7% for μ s '. Based on the validation result, the training of this ANN was considered successful and could be employed for phantom and in-vivo skin optical property recovery. It is worth noting that we also successfully generated an ANN that took reflectance at SDSs of 1, 2, and 3 mm as input parameters to calculate optical properties at an accuracy level comparable to the ANN that took reflectance at SDSs of 1 and 2 mm as the input parameters. It can thus be inferred that the role of the reflectance at SDS of 1mm is very crucial for accurate optical property recovery in the SDS range discussed in this study.  In summary, from our experience, there are two important checkpoints in producing a useful ANN. The first one is to carefully examine the uniqueness of any reflectance set in the database. The definition of uniqueness depends on the stability or noise level of the measurement system. If the uniqueness condition does not hold well, either reducing the optical property range or choosing another SDS combination is necessary. The second checkpoint is to evaluate the performance of an ANN produced from the network training function. In this step, it is important to first check the performance of the ANN using the training database; if the performance of the ANN is not satisfactory, one could adjust the configuration of the network training function and run the network training process again until a satisfactory result is obtained. Next, it is also important to randomly select a number of optical property sets to form a validation database by utilizing GPU-MCML simulations to evaluate the performance of the ANN once more. Carrying out this step is to ensure that the ANN does not only perform well at data points included in the training database. If unsatisfactory result is obtained in this step, one could increase the size of training database and enter the network training process again until one finds the performance of the ANN acceptable. Depending on the optical property range of interest, we have used training database of size in the range from 3000 to 10000. The flowchart that includes the key procedures stated above for producing a useful ANN is illustrated in Fig.  9. It is expected that the procedures described in Fig. 9 could be applied to efficiently obtain useful ANNs for various SDS combinations and optical property ranges not discussed in this study.

Recovery of liquid phantom optical properties
To test the experimental performance of the trained ANN in recovering sample optical properties, we measured the diffuse reflectance at SDSs of 1 and 2 mm of two liquid phantoms whose recipe was listed in Table 1. The benchmark absorption and reduced scattering spectra in the wavelength range from 500 to 1000 nm of the two liquid phantoms were determined using the IAD method described in section 2.4. We employed phantom LP1 as the calibration phantom to calibrate the system response of the measured diffuse reflectance of the other liquid phantom. The calibrated diffuse reflectance spectra of LP2 were sent to the ANN to recover the sample optical property spectra. It took less than 0.1 second for the ANN to calculate the absorption and reduced scattering spectra of each phantom in the full wavelength range. The results are shown as black lines in Fig. 10. The maximal deviations of the recovered optical properties from the benchmark values are 0.3% and 0.6% for μ a and μ s ', respectively.
We also employed the SDE to recover the optical property spectra of LP2 from the diffuse reflectance measured at SDSs of 1 and 2 mm, and the results are shown as green lines in Fig. 10. The maximal recovery errors are 6.1% and 25.3% for μ a and μ s ', respectively. The advantage of our ANN model over the SDE can be clearly seen here. Moreover, it can be observed that the SDE overestimates the absorption and underestimates the reduced scattering of the liquid phantom; this phenomenon might due to the crosstalk between the absorption and the reduced scattering coefficients in the diffusion theory which had also been observed in our previous study [21]. Fig. 10. (a) μ a and (b) μ s ' spectra of the liquid phantom LP2 recovered using the ANN (black lines) and the SDE (green lines). The benchmark spectra (red lines) were recovered using the IAD method.

Determination of skin optical properties and chromophore concentrations
We employed our DRS system to measure two anatomical sites of one subject's skin. The measured diffuse reflectance at SDSs of 1 and 2 mm were then fit to the SDE and the ANN to determine the skin optical properties from 500 to 1000 nm in a 2 nm step. The absorption and reduced scattering spectra of the outer thigh skin and the inner upper arm skin are plotted in Fig.  11 and Fig. 12, respectively. The error bars indicating the standard deviation of five measurements at each site are included in the plots. We can see in Fig. 11(a) and Fig. 12(a) that the absorption spectra of the two sites recovered using the ANN are distinct in the 500-600 nm region; this difference was very likely introduced by different hemoglobin contents at the two sites. From the visual observation on the two measurement sites, the upper inner arm of the subject had a denser blood vessel distribution with higher redness than that of the outer thigh; this could explain the notable absorption feature of oxygenated hemoglobin in Fig. 12(a) in the 500-600 nm region. Moreover, the absorption spectra of the two sites are very similar at the 600-1000 nm region, and the absorption peak at 980nm introduced by water can be clearly seen in both Fig. 11(a) and Fig. 12(a). It can be noted that although the absorption spectra recovered using the SDE, which are shown as green lines in Fig. 11(a) and Fig. 12(a), have a trend similar to those recovered using the ANN from 500 to 600 nm, the absorption coefficients recovered using the SDE are consistently higher than those recovered using the ANN. Based on the phantom measurement results discussed in the previous subsection, we speculate that the skin absorption was overestimated by the SDE. In contrast, due to relatively low absorption of skin in the 600 to 1000 nm region, the validity of diffusion approximation was preserved, and the SDE and the ANN produced similar absorption recovery results.  Reduced scattering spectra of the two sites recovered using the ANN are depicted as black lines in Fig. 11(b) and Fig. 12(b), and these two spectra resemble each other closely both in trend and magnitude. Contrarily, reduced scattering coefficients derived from the SDE are lower than those derived from the ANN especially in the 500 to 600 nm range. This phenomenon agrees with our phantom measurement results, and we infer that the accuracy of the reduced scattering coefficients recovered using the SDE were impaired by the high sample absorption in this wavelength region. It can also be noted in Fig. 12(b) that the reduced scattering coefficients recovered using the SDE are slightly lower than those recovered using the ANN around 980 nm. This μ s ' underestimation of the SDE could be induced by the locally prominent water absorption at 980 nm which weakens the validity of diffusion approximation; on the contrary, similar phenomenon cannot be observed in Fig. 11(b) and Fig. 12(b) for the μ s ' spectra recovered using ANN, which suggests the robustness of our ANN model. Overall, the optical properties of in-vivo skin recovered in this study are comparable to those of in-vivo or ex-vivo skin reported in the literature [3, 18,19,21].
The skin absorption spectra were linearly fit with positive constraint to a set of chromophore spectra consisting of the extinction coefficients of oxygenated hemoglobin, deoxygenated hemoglobin, melanin, and water retrieved from the Oregon Medical Laser Center website [22]. The fitting results for the absorption spectra derived from the two different models are listed in Table 2. We can see that all chromophore concentrations recovered from the SDE derived spectra are higher than those of the ANN counterpart. Our result suggests that using the SDE to calculate skin absorption coefficients would lead to overestimated chromophore

Conclusion
In this study, we proposed the use of a GPU-based Monte Carlo model for efficiently generating the data sets required for constructing an ANN that could relate the sample absorption and reduced scattering coefficients with the reflectance in real time, and the model has a precision level comparable to that of a Monte Carlo model. In addition, we suggested a robust workflow for generating and validating a functional ANN in a certain sample optical property range of interest.
We have successfully trained an ANN that can take reflectance at SDSs of 1 and 2 mm as input and accurately determine the sample optical properties in the range of 0.015< μ a <0.51 mm −1 and 0.2< μ s ' <5 mm −1 without applying spectral constraints. In the liquid phantom study, we found that the ANN could recover sample optical properties with errors less than 0.6%. In contrast, as the SDE was used to determine the optical properties of the same phantom, the optical property recovery error as high as 25.3% was observed. Finally, we employed our ANN model to determine the optical properties and chromophore concentrations of in-vivo skin of a subject and the results agree with those reported in various independent studies. It was found that the SDE overestimated the absorption and underestimated the reduced scattering of samples.
Recent remarkable advances in GPU computing have made the efficient construction of an accurate ANN feasible. It is expected the GPU computing aided ANN training workflow proposed here could be applied to expeditiously construct a functional ANN for various DRS measurement geometries designed for specific applications in which diffusion approximation does not sustain or analytical solutions of the diffusion equations cannot be found.