Rapid calculation of diffuse reflectance from a multilayered model by combination of the white Monte Carlo and adding-doubling methods

: To rapidly derive a result for diffuse reflectance from a multilayered model that is equivalent to that of a Monte-Carlo simulation (MCS), we propose a combination of a layered white MCS and the adding-doubling method. For slabs with various scattering coefficients assuming a certain anisotropy factor and without absorption, we calculate the transition matrices for light flow with respect to the incident and exit angles. From this series of precalculated transition matrices, we can calculate the transition matrices for the multilayered model with the specific anisotropy factor. The relative errors of the results of this method compared to a conventional MCS were less than 1%. We successfully used this method to estimate the chromophore concentration from the reflectance spectrum of a numerical model of skin and in vivo human skin tissue.


Introduction
The reflectance of skin noninvasively provides information about the inner conditions, such as the scattering coefficients, absorption coefficients, and chromophore concentrations [1]. When deriving this information, a Monte Carlo simulation (MCS) has often been used as the standard [2][3][4][5][6][7][8][9], since a MCS precisely follows the behavior of each photon that is scattered from or absorbed in a medium, and the range of scattering and absorption coefficients to which it can be applied is not limited. In addition, a MCS can be applied easily to arbitrary multilayered structures. A multilayered structure is unavoidable when evaluating the reflectance of skin, which consists of an epidermis that contains melanin and an underlying dermis that contains hemoglobin. Because of its high reliability, a MCS is often used to evaluate approximation methods [2][3][4][5][6][7][8][9], although it has a high computation time. In a MCS, a large number of photons are generated, and each photon numerically propagates through the skin model following designated probabilities of scattering and absorption; from a statistical analysis of these results, the reflectance and other physical quantities are then calculated. Increasing the number of photons increases the precision of the MCS, but it also increases the calculation time, which may become too long to be used for the inverse problem, imaging, or interactive tools. A sufficiently short calculation time is crucial for those applications, since it is necessary to calculate the reflectance many times. Therefore, to shorten the calculation time, several methods have been considered for approximating the results of an MCS. For example, various studies have proposed diffusion approximation [3], hybrid diffusion and two-flux approximation [4], the path-integral method [2], and an empirical method aided by MCS [10]. However, in each of these methods, the difference from the MCS depends on the type of approximation, and often the applicable range of parameters is limited. The methods based on the diffusion approximation [3,4] limit the size of the ratio of the absorption coefficient μ a to the reduced scattering coefficient μ s '. In the path-integral method, although the trajectories of the light are represented by the classical path, it is difficult to find an appropriate path for the various absorption and scattering coefficients [2]. Finally, although the empirical method aided by a MCS can estimate the chromophore concentrations over a wide range, it is not optimized for the estimation of the reflectance [10].
Several studies have attempted to reduce the time required to calculate the reflectance while keeping the results equivalent to those of an MCS [5-7, 9, 11]. These studies were based on a method known as white MCS (WMCS) [5] or single MCS [6,7,9], which were primarily developed for solving time-domain problems [5][6][7]9]. The WMCS is based on the Beer-Lambert (B-L) law, which states that the absorbance of photons traveling through a medium along a certain trajectory can be calculated from the path length and the absorption coefficient. The WMCS also utilizes similarity: In a semi-infinite and homogeneous medium without absorption, if μ s ' grows to α times the original condition, then a homothetic trajectory that is 1/α times smaller can be associated with the original one [5,6]. With similarity, we only need to calculate the MCS one time without absorption, and then for any absorption and scattering coefficient of the medium, we can derive a histogram for the path length and determine the spatial distribution of reflectance in the time domain. Not apply only to semiinfinite and homogeneous media, a WMCS can be applied to an arbitrary composite structure composed of multiple media if the path length in each medium can be calculated for each trajectory [8,11]. However, in such cases, similarity cannot be fully utilized, and the application range for a certain precalculated data set is limited. For a multilayered structure like skin, the optical path-length matrix method (OPLM) has been developed [8]. In the OPLM, a MCS is used to determine the path lengths for each photon in the epidermis and in the dermis; this is recorded on a two-dimensional histogram. The absorption coefficient is determined, and the reflectance is then calculated from the two-dimensional histogram by using the B-L law [8]. Although this method works with a multilayer model, if we want to change any of the scattering coefficients or the thickness of any of the layers, we must recalculate everything.
To resolve this issue, we developed a faster method for estimating the reflectance of an arbitrary multilayered model; its results are in close agreement with those of the MCS. Here, we omit the spatial and time-resolved information. In this method, the WMCS was modified for application to a slab by the introduction of transition matrices for light flow without absorption from the incident angle to the exit angle, based on the distribution by path length. The method is called the layered WMCS (LWMCS). Once the value is given for μ a , the transition matrices for each layer can be calculated for that μ a . The transition matrices of all the layers are then combined by using the adding-doubling (AD) method [12]. From the previously calculated values, the transition matrices can be derived for all the layers by matrix arithmetic; the calculation time is thus substantially shorter than it is with the conventional MCS (cMCS). In addition, with this method, we can change the number of layers, the thickness, μ a , μ s ', and the refractive indices without having to repeat the initial calculations. To demonstrate the benefit of the quicker calculation, we used this method to estimate the chromophore content, which requires repeated calculations of the reflectance. Fig. 1. Range of direction of incident angle Θ i and exit angle Θ j . Here, Θ i and Θ j indicate the direction between the two respective cones of θ i-1 and θ i , and of θ j-1 and θ j , respectively.

Outline of the method
In this paper, we focus on the incoming angle of the light, measured in relation to the normal to the surface, and we exclude any analysis of lateral distribution or time. We begin by deriving the transition matrices of reflectance R and transmittance T, which are based on the transition probability relative to the incident and exit angles. The incident and exit angles are separated by the angles {θ 0 ,θ 1 ,…,θ n-1 ,θ n }, where θ 0 = 0° and θ n = 90°. For convenience, the ranges [θ i-1 ,θ i ] and [θ j-1 ,θ j ] are represented by Θ i and Θ j , respectively, where i = 1,2,…,n and j = 1,2,…,n (Fig. 1). R and T are expressed as: Here, R(Θ i →Θ j ) and T(Θ i →Θ j ) represents the discretized conditional probability of the exit angle being between [θ j-1 ,θ j ] when the incident angle is Θ i , for reflection and transmission, respectively. R(Θ i →Θ j ) and T(Θ i →Θ j ) can be expressed using the conditional probability densities: Here, r(Θ i →θ j ) and t(Θ i →θ j ) are, respectively, the conditional probability densities of the exit angle θ j when the incident angle is Θ i , for reflectance and transmittance. We should note that, for simplicity, these values are not normalized to a solid angle, but they can be converted to normalized values if needed. Although the expression does not consider the azimuth angle, the reflectance and transmittance can be derived from these matrices, under the condition that either or both the incident angle and measurement angle distributions are axially symmetric. This precondition will be satisfied in the case of normal or diffuse incidence, or when the reflected light is gathered with an integrated sphere. To derive R and T for a model of skin, we developed the LWMCS [ Fig. 2]. In the LWMCS, the transition matrices of reflection and transmission involving path lengths without absorption are calculated in advance by a MCS for several μ s 's for a specific thickness (phase 1, Sec. 2.2). In this case, each element of each of the transition matrices is a histogram of the path lengths. For an arbitrary μ s ', the thickness of a layer d can be associated with one from the precalculated base data by using a scale factor; this step is called resizing (phase 2, Sec. 2.3). To take account of μ a , the B-L law was used, and the absorbance for a certain path length was derived from μ a and the path length. The intensity for each path length was calculated, and these values were then accumulated. This phase is called coloring (phase 3, Sec. 2.4). Since the accumulation is done across all path length, the results are in the form of Eq. (1) and Eq. (2) (i.e., each element becomes a scalar). Then, from the derived transition matrices of each layer, the transition matrices of the multilayered model can be calculated by using the AD method [12]; this is called lamination (phase 4, Sec. 2.5). The refractive index is not considered in the first three phases, but it is considered in phase 4. The protocol described above is called the generalized LWMCS (gLWMCS).
Phase 1 requires repeated evaluations of a MCS; however, once the parameters of a multilayered model are defined, the reflectance can be calculated quickly. This allows us to shorten the trial-and-error process of adjusting the parameters. The method can also be applied to optimization problems that require recursive calculations of the reflectance; an example of this would be estimating chromophore concentrations.
To estimate the error of the interpolation as a function of ξ in gLWMCS, we also examined a specified LWMCS (sLWMCS). In an sLWMCS, the transition matrices from phase 1 are specified for a particular μ s ' and d in the model used; this means that the interpolation in phase 3 is not necessary. However, phase 1 must be repeated when μ s ' or d are changed in any layer.

Precalculation of the transition matrices of a single layer (phase 1)
The aim of phase 1 is to use the path length without absorption to calculate R and T for each layer and for several different scattering powers with the MCS. The results will then be used to quickly derive R and T for layers with an arbitrary absorption coefficient. Assuming a path length l and layer thickness d 0 , the normalized path length 0 l d ζ = can be separated as {ζ 0 ,ζ 1 ,…,ζ m-1 , ζ m } where ζ 0 = 0 and ζ m = ∞. For convenience, the range [ζ j-1 ,ζ j ] is represented by Ζ j , where j = 1,2,…,m. A single MCS produced one column each of R and T for path lengths without absorption and for a particular pair of Θ i and the normalized reduced scattering coefficient ξ = μ s '·d 0 , as follows: Here, R W,ξ (Θ i →Θ j ,Ζ k ) and T W,ξ (Θ i →Θ j ,Ζ k ) are the conditional probabilities that the exit angle and the path length are between [θ j-1 ,θ j ] and [ζ k-1 , ζ k ], respectively, when the incident angle is Θ i . The subscript W represents white (no absorption). Since those transition matrices are defined by the product of μ s ' and d 0 as mentioned in Sec. 2.3, d 0 can be set as an arbitrary constant in the MCS. Each row in R W,ξ (Θ i ) and T W,ξ (Θ i ) forms a histogram of the path length. Since there is no absorption in MCS of phase 1, the photons will never be attenuated in slabs. Therefore, the trace of a photon is terminated when the photon exits from a slab. In sLWMCS, the particular values of ξ in the skin model are used for the precalculation. However, in the gLWMCS, the value of ξ used in the precalculation was not the same as the value used in the model.

Resizing the thickness of a layer (phase 2)
Since, from similarity, the transition matrices can be associated with the product of μ s ' and the thickness d, it is not necessary that d has the same thickness as d 0 in the precalculation. Thus, the ratio α of d to d 0 is calculated and is used as a scaling factor. The concept of similarity is described as follows. The transition matrices for the case with reduced scattering coefficient s μ α ′ , absorption coefficient a μ α , and thickness d α ⋅ are the same as those for μ s ', μ a , and d, where α is an arbitrary factor. The probability that a photon follows an arbitrary trajectory magnified by α, with s μ α ′ , a μ α , and d α ⋅ is the same as that of a photon following the original trajectory in a layer with μ s ', μ a , and d ( Fig. 3) [5,6]. This means that the transition matrices in phase 1 can be expressed as a function of the normalized reduced scattering coefficient ξ = μ s '·d and the scale factor α = d/d 0 . The form of the transition matrices as a function of ξ do not change with α, while ξ changes not only with μ s ' but also with d. Therefore, if the transition matrices for the series of ξ = {Ξ 1 , Ξ 2 ,…, Ξ m } are prepared as in phase 1, appropriate transition matrices can be derived for arbitrary μ s ', μ a , and d. Although Ξ i has a discrete value, we can derive R and T for a given ξ by interpolating between the largest and smallest values of ξ in Ξ i . In phase 2, an appropriate set of transition matrices is chosen. For a particular layer with μ s ' and d, we choose the precalculated data set with ξ = μ s '·d. In the sLWMCS, the transition matrices for the particular ξ in the skin model are calculated in phase 1. In the gLWMCS, the largest and smallest values of ξ near to the value of Ξ i are obtained in phase 1, and the results are then interpolated in phase 3.

Coloring of transition matrices for a single layer (phase 3)
In the LWMCS, the absorbance for a particular absorption coefficient and path length can be calculated in accordance with the B-L law. Assuming that an absorption coefficient μ a , thickness d, and normalized path length ζ = l/d are given, the absorbance A can be expressed using the normalized absorption coefficient μ a ·d as Therefore, the total amount of reflected and transmitted light can be expressed as a function of the angle, as follows: This vector becomes a column of the transition matrix for the associated incident angle. If we calculate for all incident angles and then combine them, the transition matrices of a layer with an arbitrary absorption coefficient are derived as: In the sLWMCS, a particular ξ in the skin model is used in phase 1, therefore the transition matrices can be simply calculated from Eqs. (9) and (10). In the gLWMCS, the value of ξ in the model does not usually match the value of ξ used in phase 1; instead, R and T are derived by simple linear interpolation with Eqs. (11) and (12). For a particular ξ, assuming that the transition matrices for Ξ i were calculated in phase 1, the transition matrices for ξ are calculated by linear interpolation, as follows: where Ξ + and Ξare the smallest and largest Ξ i near ξ. We will number the layers from the top down. We will not allow the refractive index to change at the boundaries, and instead, any change in the refractive index will be considered to take place in a virtual layer. This means when there is a change in the refractive index, the boundary (in the usual sense) should be counted twice, and the effect of the refractive index ratio is treated as a characteristic of a virtual layer at that location. R, T and L are defined as in Fig. 4(a). The notation of the transition matrix is expanded to multiple layers, and the transition matrix of the combined layers, from boundaries k to l, are expressed as R k:l and T k:l [ Fig. 4

Lamination of layers with the adding-doubling method (phase 4)
Our aim in phase 4 is to derive the matrix R 1:N that contains the transition probability from an arbitrary incident angle to an arbitrary exit angle for the layers with boundaries 1 to N. We begin by combining two pairs of transition matrices of two adjacent layers into a single pair of transition matrices. Then, by using recursion, the transition matrices of an arbitrary number of layers can be combined into a pair of transition matrices.
The reflection and transmission transition matrices for a layer between the boundaries k and k + 1 (Fig. 4) can be written as: Here, R k:k+1 (Θ i →Θ j ) and T k:k+1 (Θ i →Θ j ) represent the discretized conditional transition probability of reflected and transmitted light, respectively, and they can be expressed as: in the same manner as in Eqs. (3) and (4). Here, r k:k+1 (Θ i →θ) and t k:k+1 (Θ i →θ) are the conditional probability density of the exit angle of reflected and transmitted light, respectively, when the incident angle is Θ i .
Assuming depolarized light, the transition matrix of a virtual layer k:k + 1 with a change in refractive indices (a single boundary in the usual sense) can be derived by using the Fresnel equation. If , where N 1 and N 2 are the refractive index of the reflected side and the transmitted side, respectively, the transition matrix for reflectance can be expressed as follows: Here, Β i is the transmission angle, which can be derived by using Snell's law ( , and δ i,j is the Kronecker delta. With regard to the transition matrix for transmission, the following relationship is useful: The transition matrix of transmission can be derived if we assume the transmitted light is distributed uniformly across the range of [β i-,β i+ ], where for the incident angles θ i-1 and θ i , we have the transmission angles β i-and β i+ , respectively. When The angular dependency of the light flow at a boundary k is expressed as , We can obtain similar relationships for the boundary 2:3. In addition, by considering two layers to be one composited layer, we can obtain the same kind of equation for the composited boundary 1:3. By eliminating L 2,↓ and L 2,↑ , noting that R and T are noncommutative, and noting that the equations should hold for arbitrary values of 1, L ↓  , 1, L ↑  , 3, L ↓  , and 3, L ↑  , R 1:3 and T 1:3 can be expressed in terms of R 1:2 , R 2:1 , R 2:3 , T 1:2 , T 2:1 , and T 2:3 as follows [12]: ( ) where E is the unit matrix. By interchanging 1 and 3, we can obtain similar equations for R 3:1 and T 3:1 . Once R 1:h and T 1:h have been derived, by substituting R 1:h , R h:h+1 , T 1:h , and T h:h+1 for R 1:2 , R 2:3 , T 1:2 , and T 2:3 , respectively, R 1:h+1 and T 1:h+1 can be obtained. We can obtain R h+1:1 and T h+1:1 in a similar way. Thus, the pairs of transition matrices can be derived recursively for an arbitrary number of layers. We adopted the adaptive bin with non-equal divisions of logarithm for binning of path length as shown in Fig. 5. With a particular constant σ R , the normalized path length ζ is transformed into η R = exp(σ R × ζ) for R. For T, since the path length is greater than or equal to the thickness of layer d, the path length was transformed into η T = exp(σ T × (ζ-1)), noting that ζ = l/d. Then, η R and η T (0 to 1) are divided by the following manner: a certain point η c is set between 0 and 1. The region from 0 to η c and the region from η c to 1 are then each divided separately into m l or m h sections, respectively. This binning can be characterized by four parameters (σ, η c , m l , and m h ). For R and T, the parameters are expressed as (σ R , η cR , m lR , m hR ) and (σ T , η cT , m lT , m hT ), respectively. By doing this, the bin widths of ζ for shorter path lengths became narrower than those for longer path lengths. Another benefit of the transformation is that the semi-infinite range of path lengths [0, ∞] is transformed to a finite range [0, 1], which enables fading of the binning.

Skin model
We used a two-layered model: an epidermis with uniform melanin, and an underlying dermis with uniform oxygenated and deoxygenated hemoglobin. Each parameters are shown in Table  1. Here, C m , C oh , and C dh are the concentrations in each layer of melanin, oxygenated hemoglobin, and deoxygenated hemoglobin, respectively. For ε m (λ), we used the average absorption coefficient of the monomer melanosome with a concentration of 1 mole/liter; this was approximated as 6.6 × 10 11 × λ -3.33 in cm −1 , where the unit of λ was nm [13]. For ε oh and ε dh , we used the extinction coefficients of oxygenated and deoxygenated hemoglobin, respectively, converted to a hematocrit concentration of 45 [14]. The scales of C m , C oh , and C dh were the ratio against the concentrations under which ε m (λ), ε oh (λ), and ε dh (λ), respectively, were derived. The scattering coefficient μ s (λ) was derived from the reduced scattering coefficient μ s '(λ) [15][16][17] and the anisotropy factor g [18] as μ s (λ) = μ s '(λ) /(1-g). For the error estimation presented in Sec. 3.4.1, μ a was set from 0 to290cm −1 in intervals of 10cm −1 for the epidermis and was set from 0 to 29cm −1 in 1cm −1 intervals for the dermis. The range of μ a from 0 to 290cm −1 of the epidermis covers 0 to 20% of C m, and that from 0 to 29cm −1 of the dermis covers 0 to 1.1% of C oh, and 0 to 1.0% of C dh for wavelengths in the range of 400 to 700 nm. These concentrations are the actual values for skin that is light to moderately pigmented [16]. Other than μ a , the parameters of the skin model in Sec. 3.4.1 were the same as those of the skin model described above.

MCS
For calculating the MCS, the program MCML [19] was modified to record the histograms of reflectance and transmittance against the path length, and to perform the calculations for an arbitrary incident angle. With the modified MCML, we executed the cMCS and phase 1 of the sLWMCS and gLWMCS. The number of photons for each condition was 10 5 [20]. For the standard condition for the skin model, the incident angle was set to 0°, and the reflectance was integrated over all exit angles.

Conventional MCS (cMCS)
As a standard, we used conventional MCS (cMCS). In the cMCS, the parameters μ s , g, μ a , d, and the refractive index of each layer, were set as described in Sec. 3.2, and photons were considered to travel in a multilayered model of skin.

sLWMCS
Following the procedure described in Sec. 2, we used a sLWMCS to calculate R and T for the skin model in Sec. 3.2. Phase 1 was calculated by modified MCML with the condition described in Table 2. The divisions of binning for the incident angle and exit angle θ i were both set to {0°,1°,…,89°,90°}, and the representative values Θ i were defined as {0.5°,1.5°,…,88.5°,89.5°}. Phases 2 to 4 were executed using MATLAB ® (Mathworks, Natick, MA, USA). The incident angle (for the whole model) was set to 0.5°, which approximated the 0° incidence angle of the standard condition. The thickness d of the slabs in the models of phase 1 were set to 0.06mm and 4.94mm for epidermis and dermis, respectively, and the μ s at each wavelength was used.
As mentioned above, the appropriate values of σ R and σ T depend on μ s ', g, and d. The value of σ was determined visually to be such that the elements of the histogram did not concentrate in the lowest and highest bins (around 0 and 1 for η).The values for {σ R , η cR , m lR , m hR } and {σ T , η cT , m lT , m hT } for binning of path length were empirically derived and summarized in Table 3. As the output of phase 1, the size of the transition matrices with path length was 90 × 90 × 100. As the output of phase 3, the transition matrices lost the information about path length, then the size of them became 90 × 90.
The MCS was executed 180 times (90 incident angles × 2 layers) for each of 31 wavelength from 400nm to 700nm in 10nm interval, 5580 times in total. The total process time was about 16 days with our computational environment. The size of the base data set became 63 MB for epidermis, and 74 MB for dermis as MATLAB data files.

gLWMCS
Following the procedure described in Sec. 2, R and T of the skin model in Sec. 3.2 were calculated by using the gLWMCS. Phase 1 was calculated by modified MCML with the condition described in Table 4. The same values that were used for the sLWMCS (Sec. 3.3.2) were used for the discretization of the incident and exit angles. The value of Ξ i was set to 10 (i/10-1) for i in the range of 0 to 30 in each interval, which means Ξ i varied from 0.1 to 100 with equal intervals on a log scale in gLWMCS. The range of Ξ i is comparable to 238-238095cm −1 in μ s for d = 0.06mm (where d is for the epidermis in our skin model), and to 3-2886cm −1 in μ s for d = 4.94mm (where d is for the dermis in our skin model). Since μ s is 1473.3cm −1 at 400nm and monotonically decreases to 273.3cm −1 at 700nm in our skin model, the variation covers μ s in the epidermis and dermis at wavelengths from 400 to 700 nm. Due to the similarity, d 0 does not affect the results basically, however, it has an effect on the appropriate values for binning parameters of path length. The MCS was repeated for each Ξ i ; this was achieved by setting d to 1cm and μ s to Ξ i /(1-g) in the modified MCML. Phases 2 to 4 were executed using MATLAB ® . The value of σ was defined visually to be such that the elements of the histogram did not concentrate in the lowest and highest bins (around 0 and 1 in the η); this was the same as for the sLWMCS. The parameters for binning of path length are described in Table 3. The size of the transition matrices was the same with sLWMCS.
The MCS was executed 90 times (90 incident angles) for each of 31 Ξ i from 10 (0/10-1) to 10 (30/10-1) , 2790 times in total. The total process time was about 12 days with our computational environment. The size of the base data set became 70 MB as MATLAB data files.

Error estimation
To estimate the error, the skin model in Sec. 3.2 was used. For μ s , we used the values at 400, 500, 600, and 700 nm, which were, respectively, 1473.3, 712.7, 414.9, and 273.3 cm −1 . The value of μ a for the epidermis was set to range from 0 to 290 cm −1 in intervals of width 10, and μ a for the dermis was set to range from 0 to 29 cm −1 in intervals of width 1. We calculated the relative error E for the cMCS, sLWMCS, and gLWMCS. The relative error of reflectance from a method evaluated was defined as follows: reflectance from a method evaluated reflectance from cMCS .
For the cMCS, the reflectance was calculated again under each condition besides the calculation as the reference, and the results from the second trial were used to evaluate the method. Since the difference between the first calculation and the second one is the sequence of random numbers, E for the cMCS demonstrates statistical error. The difference between sLWMCS and the gLWMCS demonstrates the effect of interpolation.
The reflectance spectra representing the actual human skin under the conditions of normal, occlusion, and post-occlusion reactive hyperemia were simulated with cMCS, sLWMCS, and gLWMCS. In those simulations, the concentrations of melanin, oxygenated hemoglobin, and deoxygenated hemoglobin were selected as follows: normal condition, C m = 6.15%, C oh = 0.10% and C dh = 0.10%; occluded condition, C m = 6.21%, C oh = 0.00% and C dh = 0.19%; post-occlusion reactive hyperemia, C m = 6.20%, C oh = 0.33% and C dh = 0.13%. The values of C m , C oh , and C dh were derived from the estimated values of human skin at 0, 300, and 510 s shown in Fig. 10 of Ref [10]. The MCML for the cMCS and phase 1 of each LWMCS was executed on Windows 8 ® (64bit). The PC had 4 GB of memory, and the CPU was Core i5-3230M ® , 2.60GHz (Intel, Santa Clara, CA, USA). For the calculation of phases 2 to 4 in the sLWMCS and gLWMCS, we used MATLAB ® 7.0.4 on a Windows XP ® emulator in a VMWare Player ® 5.0.2 on Windows 8 (64bit); note that MATLAB ® 7.0.4 does not work in Windows 8. Processing time was evaluated with the MATLAB profiler. Windows XP ® and Windows 8 ® are products of Microsoft (Redmond, WA, USA), and VMWare Player ® is the product of VMWare ® (Palo Alto, CA, USA). The memory allocated to the Windows XP emulator was 2GB, and 3 of the 4 CPU cores were allocated to the emulator.

Estimating chromophore concentrations from a spectrum
The inverse problem, such as the estimation of chromophore concentrations, is where rapid calculations become important, because the spectrum must be calculated iteratively under several different conditions. From a given spectrum, we searched for the chromophore concentrations as an optimization problem that minimized the evaluation function F. The optimization problem was solved with the built-in function fminsearch of MATLAB, in which the Nelder-Mead simplex method is implemented. The skin model of Sec. 3.2 was used, and the wavelength was set from 400 to 700 nm at intervals of 10 nm. However, 500 to 600 nm at intervals of 10 nm was used for the actual skin spectra in the experiment of Sec. 3.4.3. Because this range of wavelength is often used for chromophore estimation due to the distinguishable difference between oxygenated and deoxygenated hemoglobin concentrations. In this case, the ranges of 400 to 490nm and 610 to 700nm were extrapolated using respective LWMCS with the estimated chromophore concentrations. The evaluation function to be minimized in the optimization problem was defined as follows: S c c c λ is the estimated spectrum for given chromophore concentrations.
With the spectrum from the cMCS as the input, we used the sLWMCS and gLWMCS to search for the chromophore concentrations that minimized F. The results were compared with the input of the cMCS. For the input of cMCS, the concentrations of melanin, oxygenated hemoglobin, and deoxygenated hemoglobin were set to 1%, 0.1%, and 0.1%, respectively, which are in the range of Asian skin, and the initial values of the iteration were set to 2%, 0.5%, and 0.5%, respectively. The processing time and number of iterations were evaluated with the MATLAB profiler.

Experiments with in vivo human skin
To confirm the applicability of the proposed method to an actual human skin, we performed the time series measurements of melanin, oxygenated hemoglobin, and deoxygenated hemoglobin on in vivo human skin tissue during cuff occlusion. Experimental procedure was approved by the Ethical Committee for experiments involving human subjects at Tokyo University of Agriculture and Technology. Written informed consent was obtained from a Japanese male subject. A pressure cuff was fixed around the upper arms of the subject and pressured up to 250 mmHg was applied by using a rapid cuff inflator (E-20, D.E. Hokanson Inc., Bellevue, WA, USA).  Figure 6 schematically shows experimental setup for measuring diffuse reflectance spectra. A 150-W halogen-lamp light source (LA-150SAE, Hayashi Watch Works Co., Ltd, Tokyo, Japan), which covers the visible wavelength range from 400 to 700 nm, illuminates the skin surface of the forearm via a light guide and lens with a spot diameter of 4 mm. The diameter and focal length of the lens are 54 mm and 100 mm, respectively. In the time series measurement, forearm of the subject was fixed on the sample stage to prevent motion artifacts. The skin surface was fixed at the sample port of an integrating sphere (RT-060-SF, Labsphere Inc., North Sutton, NH, USA). The detected area of the skin surface was circular, with a diameter of 22 mm. Before making the measurements of reflectance spectra, we measured the intensity of the light that is incident on skin to collect the reflectance spectrum by using the optical power meter. In the 400-700 nm range, the maximum power was 495 μW at 510 nm, and the minimum power was 45 μW at 700 nm. Light diffusely reflected from this area was received at the input face of an optical fiber having a core-diameter of 400 μm placed at the detector port of the integrating sphere. The optical fiber transmits the received light into a multichannel spectrometer (SD-2000, Ocean Optics Inc., Dunedin, FL, USA), which measures reflectance spectra in the visible wavelength range under the control of a personal computer. A standard white diffuser with 99% reflectance (SRS-99-020, Labsphere Inc.) was used to correct the spectral intensity of light source and the response of spectrometer. In the in vivo human skin measurement, a single reflectance spectrum was obtained by averaging ten successive recordings of reflectance spectrum, in which one recording was made with the integration time of 200 ms. Therefore, the acquisition of the single reflectance spectrum needs 2 s in total. In the time series measurement of in vivo reflectance spectra before, during and after cuff occlusion, the acquisition of each single reflectance spectrum was repeatedly made at 30-s intervals for 14 min (840 s). According to the estimation procedure for chromophore concentrations described in Sec. 3.4.2, the concentrations of melanin C m , oxygenated hemoglobin C oh , and deoxygenated hemoglobin C dh for each time point were obtained. The concentration of total hemoglobin C th was simply calculated as sum of C oh and C dh . The tissue oxygen saturation StO 2 was calculated as StO 2% = 100 × {C oh /(C oh + C dh )}.

Error estimation
In Fig. 7, the values of E for the cMCS, sLWMCS, gLWMCS are plotted against the cMCS, for the various conditions of μ a in the epidermis and dermis. The average and standard deviation of E over the whole range of μ a are shown in Table 5. On average, the absolute values of E for both the sLWMCS and the gLWMCS were one order of magnitude larger than those of the cMCS. The standard deviation of E for both the sLWMCS and the gLWMCS were smaller than that for with the cMCS. It is seen that the values of E increase on average for the values of higher μ a in Fig. 7. This is probably due to the discretization error in binning of path length. On the other hand, the variations in E also increase for the values of higher μ a in Fig. 7. This is probably due to the statistical error originated from the Monte Carlo methods. The number of photons reemitted from the skin model will decrease as the value of μ a increases. In such a case, the error in reflectance will be statistically increased. This limits the accuracy of sLWMCS and gLWMCS as well as cMCS. Increase in the number of incident photons used in the phase 1 will overcome this problem. The simulated reflectance spectra representing the actual human skin under the conditions of normal, occlusion, and postocclusion reactive hyperemia obtained from cMCS, sLWMCS, and gLWMCS are shown in Fig. 8. The processing time of the cMCS was 25270 sec. for the derivation of a spectrum (31 wavelengths) on average. The processing time of phases 2 to 4 in the sLWMCS and gLWMCS were 3 sec. and 5 sec., respectively, for the derivation of each spectrum, which are, respectively, about 8000 times and 5000 times faster than that of the cMCS. As can be seen in Fig. 8, the skin reflectance spectrums from both the sLWMCS and the gLWMCS were practically identical to that from the cMCS. The difference between the cMCS and the sLWMCS or gLWMCS was more than two digits smaller than the reflectance, which reflects the precision shown in Table 5.

Estimating chromophore concentrations from the cMCS spectrum
The results are summarized in Table 6. With the sLWMCS and gLWMCS, the estimated concentrations of melanin, oxygenated hemoglobin, and deoxygenated hemoglobin are shown, for the spectrum from the cMCS with the concentrations of 1%, 0.1%, and 0.1%, respectively. The processing time and iteration count are also shown. The processing time of phase 2 was described as 0 sec., because only the homothetic ratio was calculated in the phase.   Figure 9(a) shows the comparisons between the typical measured reflectance spectrum obtained from the actual human skin under the normal, occluded and post-occluded condition with the reconstructed reflectance spectra from the chromophore concentrations estimated by the optimization method with the gLWMCS. The reconstructed reflectance spectrum with gLWMCS agrees reasonably with the measured reflectance spectrum in the wavelength range for the fitting (from 500 to 600 nm). Figure 9(b) shows the time courses of C m , C oh , C dh , C th , and StO 2 during cuff occlusion at 250 mm Hg. The average values of C m and C th were 6.24 ± 0.06% and 0.211 ± 0.008%, respectively, in pre-occlusion (normal), which are close to typical values for Japanese subjects reported in the literatures [10,21,22]. The average value of 44.3 ± 3.4% for StO 2 in pre-occlusion agrees with the mean blood oxygen saturation of normal human subjects (range 30.2-52.4%) reported in the literatures [23]. During cuff occlusion, C oh and C dh decreased and increased, respectively. The value of StO 2 exhibited the well-known deoxygenation curve, in which the oxygen saturation falls exponentially. The slight increase in C th probably has a physiological cause. This is because, during occlusion, the venous outflow is reduced more than the arterial inflow. After the occlusion, C th increased substantially due to the endothelial function. Despite the remarkable changes in C oh , C dh , C th , and StO 2 , the value of C m , which is independent of temporary hemodynamics, remained almost unchanged during the measurements. In this way physiological conditions of actual human skin tissue were successfully monitored by using the proposed method.

Discussion
As can be seen in Fig. 7 and Table 5, the values of the E of the sLWMCS and gLWMCS were less than 1% across the values of μ s ' and μ a that are included in the range of actual human skin, from light to moderately pigmented [16]. This shows that we achieved our aim of quickly obtaining results that are in good agreement with those of the cMCS on average. The standard deviations of these methods were smaller than that of the cMCS. This means that the statistical errors in the sLWMCS and the gLWMCS are smaller than that of the cMCS, when the same number of photons is used. In the MCS, the reflectance is calculated stochastically, and so statistical errors result from the finite number of photons. In the LWMCSs, the MCS was executed for each layer and for each incident angle, and so the total number of photons was much larger than it was in the cMCS. On the other hand, the average of the values of E from the LWMCSs was larger than that of the cMCS, which is probably due to the discretization error in the path length and the angles in both the sLWMCS and gLWMCS. The errors of the sLWMCS are close to those of the gLWMCS. This means ξ was successfully estimated by interpolation. The value of σ in the gLWMCS also seems to be appropriate, based on a comparison of the errors of the gLWMCS and the sLWMCS with the cMCS. In Fig. 9(a), there is a good agreement in the wavelength range from 500 to 600nm, but not that good for wavelengths near 400 to 450nm and 650 to700nm. The discrepancy between the measured reflectance spectrum and the fitted reflectance spectrum can be observed not only in the shorter wavelength region at which light is strongly absorbed by hemoglobin and melanin but also in the longer wavelength region. Therefore, the discrepancy between the measured reflectance spectrum and the fitted reflectance spectrum is probably due to the difference between the simulation model and actual human skin rather than the value of higher absorption coefficient. In real cases, the thickness of the layers and the other structural attributes of actual skin are not exactly the same as in the skin model that we used. We also did not consider polarization, although, according to Fresnel's law, transmitted and reflected light are partially polarized, and this affects reflectance. The calculations of the sLWMCS and gLWMCS were four orders of magnitude faster than that of the cMCS in our calculation model and computational environment. The calculation time will become even shorter if the latest MATLAB is used without an emulator. It is worth noting that the processing times of these methods will not be affected by the number of photons, the absorption coefficient, the scattering coefficient, or the anisotropy factor. In contrast, the calculation time of cMCS strongly depends on those parameters. Although the processing time of phase 1 can be varied by those parameters, phase 1 is only calculated once, and thus, it will not affect the processing time of phases 2, 3, and 4. The statistic errors can be minimized whatever we need without taking time to phase 2 to 4, by increasing the number of photons of MCS in phase 1.
The time to calculate the reflectance was sufficiently short that the chromophore concentrations could be successfully estimated by solving an optimization problem ( Table 6). The calculation time of the estimation was about ten to twenty minutes in the environment we used; this is comparable to that for a single cMCS. This may be acceptable for estimating chromophore concentrations, and it is better than the several weeks required for doing this with the cMCS, which has a processing time that is more than 5000 times longer than that of the gLWMCS. The optimization problem for the chromophore estimation was adapted to actual skin reflectance spectrum. The same time courses of changes in chromopohore concentrations with the result reported in literature [10,21].
We prefer gLWMCS in terms of flexibility of a base data set, although gLWMCS requires longer processing time than sLWMCS. gLWMCS can be extended to handle more complex skin model such as the multilayer model including a subcutaneous fat layer, without additional calculation of phase 1. Even when ξ of the new layer is larger than any Ξ i , transition matrices for larger ξ can be created by doubling a layer of Ξ i (or multiplying more). This is a significant characteristic of gLWMCS. Only in the case that all the parameters are considered as constant except μ a , sLWMCS is superior to gLWMCS because of the shorter processing time. However, if we want to treat any one of those parameters such as thickness or μ s as variable, the additional sets of base data must be prepared in sLWMCS, and thus, the calculation time for the preparation become crucial. The OPLM [8] was also developed in order to quickly derive the reflectance with values identical to those obtained from a cMCS of a multilayered model. However, the gLWMCS has apparent advantages over the OPLM in terms of the ease with which the structure of the model can be changed. With the OPLM, if there are any changes in the parameters or the number of layers, the baseline data must be recalculated. The gLWMCS, however, allows us to alter the μ s ', thickness, or refractive index of any layer, or to add more layers, without the need to recalculate the baseline data. In terms of calculation time, we note that for the OPLM, the number of dimensions of the path length matrix increases if we expand the number of layers. This causes the computational cost to increase exponentially. With the gLWMCS, the additional calculations are only required for creating the transition matrices and for combining the additional layer with the other layers, so the computational cost only increases linearly with the number of layers. According to a study of the OPLM [8], the errors were larger than those of the sLWMCS and gLWMCS, especially for longer wavelengths [8]. However, this is probably due to the equal divisions in the binning of the path length, and we assume it would be improved if the binning were optimized.
For some applications, such as imaging, the calculation time is still too long. Since most of this is due to arithmetic operations on matrices in phases 3 and 4, the time could be shortened by reducing the size of the matrices. However, if this is done by simply decreasing the number of bins in the path length histograms, the accuracy will deteriorate. Optimizing the binning strategy (i.e., finding a better binning strategy than equal logarithmic widths with dual partitioning) will be effective for improving both the calculation time and the accuracy. The incident and exit angles were discretized with equal angle intervals, but other methods, such as equal cosine intervals [12], could be considered as possible alternatives to maximize the performance. Increasing the number of photons in the MCS will improve the error; it will lengthen the precalculation (phase 1), but not the subsequent calculations (phases 2 to 4). When estimating the chromophore concentration, the initial values affect the number of iterations required. We intentionally set the starting values far from the actual values to demonstrate the robustness of the methods. However, if we determine appropriate starting values by using an approximation method, such as a diffusion model, the iteration count, and therefore the processing time, will be reduced. In addition, we could use fewer than 31 wavelengths, and the method used for interpolation could be improved. In the gLWMCS, the transition matrices are calculated twice for each layer (for Ξ + and Ξ -), but they are only calculated once in the sLWMCS. This is the primary reason for the difference in their processing times. If the transition matrices could be interpolated at step 2, then they would not need to be calculated twice. Upgrading the computer hardware and software would also reduce the calculation time, since the environment we used was not the best currently available.
The selection of the Ξ i could also affect the accuracy of the calculation. A narrower interval of the Ξ i will enlarge the amount of baseline data and lengthen the time to complete phase 1, but it will improve the results. In optimizing the interval of the Ξ i , the dominant type of scattering should be taken into account. For low Ξ i , single scattering is dominant, while multiple scattering is dominant for high Ξ i . In the intermediate region, the linear interpolation error is expected to become larger due to the phase transition. Therefore, for the most accurate calculations of reflectance and transmittance, the interval of the Ξ i in the intermediate region should be smaller than those of the other regions.
In the cMCS for the multi-layered medium, the photons propagate throughout the medium due to the scattering and these scattered photons can take a random path across multiple layers. On the other hand, in LWMCS, the tracking of scattered photons is done for a single layer only in phase 1. In the lamination process of multiple layers in phase 4, only the energy transfer of photons across the different layers is treated to calculate diffuse reflectance and total transmittance. This lamination process in phase 4 inherits the characteristic of AD method [12]. It cannot be used for time-domain analysis, and it cannot derive the reflectance as a function of the distance from the incident point. The model is limited to a homogeneous multilayer model in which the interfaces are parallel and infinite, and each layer must have uniform scattering and absorption properties, and a uniform refractive index. Either or both the incident light or the detected light must be axially symmetric. However, if we limit the model to the condition that g = 0, the method may be adapted to a nonsymmetric condition by separate treatment of scattered light and nonscattered light. This works because once the light has been scattered, the direction is independent of the azimuthal angle of the incident light.
This method also inherits benefits of the AD method [12], as follows. The characteristics of the scattering properties of each layer are derived separately, which is a remarkable characteristic and enables a physical interpretation of the results. In addition, The fluence rate at an arbitrary depth or absorption in a particular layer can be derived from R and T for the layers above the depth of interest as well as those for the whole medium based on the theory of AD method [12].
We only varied μ s to model the scattering characteristics. If we wish to consider varying the anisotropy, the base transition matrices can be pre-calculated for various values of the anisotropy coefficient, in addition to various values of μ s . In this case, the calculation time will increase since the interpolation will become two dimensional.

Conclusion
We developed a method to quickly derive results equivalent to those of the cMCS for the calculation of reflectance from a layered medium and with a broad range of applicable parameters. Once the baseline data set is prepared, it can be applied to various multilayer models, including human skin. The difference between the reflectance calculated by our method and that of the cMCS was less than about 1% over the range of μ s ' and μ a of light to moderately pigmented skin. In addition, we successfully used the method to estimate the chromophore concentration from a spectrum; this required iterative calculations of the reflectance. We expect that this method can substitute for the cMCS when the total amount of reflection is the concern. With this fast and accurate method, the parameters to be fit in the inverse problem can include the thickness, the scattering coefficient and anisotropy factor, which are usually assumed to be constant due to limitations in computation time and resources. We discussed ways in which the processing time can be further shortened while still maintaining accuracy, and this should be further studied in future work.