1. Introduction

Diffuse reflectance spectroscopy (DRS) is a noninvasive technique with good temporal resolution and rich spectral information. DRS has been applied to study the absorption and scattering properties of biological tissues, which can be further translated to clinically useful tissue functional information, such hemoglobin, water, and lipid concentrations [1, 2]. The probing volume of DRS depends on the average photon travel lengths in tissue; thus, many studies manipulated the separation between source and detector to adjust the DRS probing volume. For instance, DRS with long and short SDSs (SDSs) was employed to study deep tissues, such as muscle and brain [3, 4], and superficial tissues, such as skin [2], respectively. In addition, DRS is a model based technique which involves solving an inverse problem by utilizing a nonlinear regression algorithm to fit the measurement data to a photon transport model to obtain the optical properties of samples. The adequacy of the photon transport model affects the accuracy of the recovered optical properties. For example, standard diffusion equation derived from radiative transport theory can be used to model the transport of photons that are in the diffusion regime [5]. To reach this regime, photons have to experience sufficient scattering events and thus have long travel lengths (typically greater than 10 transport mean free paths) [6]. In contrary, the standard diffusion model was shown to be ineffective in computing diffuse reflectance at SDSs shorter than 5mm [7]. Therefore, DRS works in conjunction with a standard diffusion equation at SDSs in the range from 10 to 30 mm has been successfully demonstrated for quantifying the optical properties of deep tissues [3, 8].

To investigate deep tissues using DRS, most studies employed a standard diffusion equation that assumed a semi-infinite sample geometry. Such semi-infinite assumption in modeling has been reported to introduce large errors for the determination of the optical properties of layered samples. Farrell et al. studied the influence of the top layer on the bottom layer optical properties recovery, and they reported that the semi-infinite assumption could lead to a recovery error higher than 100% [9]. Because most biological tissues have a layered structure, diffusion models that take layered sample structure into account have been proposed. Martelli et al. used the eigenfunction method to solve the time-dependent diffusion equation for two- and three-layered samples [10]. Liemert and Kienle found an analytical solution of the n-layered diffusion model in the Fourier domain [11]. Although these multi-layered diffusion models can precisely describe photon transport in the diffusion regime, several studies indicated that using these models to determine the properties of a layered sample from measurements was still a difficult task [3, 12, 13]. Pham et al. and Kienle et al. used layered diffusion models with frequency-domain and time resolved DRS system, respectively, to determine five parameters (absorption and reduced scattering coefficients of the first and second layers, and the thickness of the first layer) of two-layered samples. They both reported that simultaneous recovery of five parameters was possible only when the initial values provided for the inverse problem were very close to the true values; otherwise the errors of the recovered values would be unacceptably large [12, 13]. As stated by Pham, the major reason for this phenomenon is that the inverse problem for the five-parameter fitting does not have a unique solution [13]. The situation may become worse when the sample under investigation has a layer number larger than two.

The goal of this study was to find a method to enhance the ability of a typical multi-layered diffusion model in reliable recovery of layered sample optical properties. In this paper, we will reveal a robust iteration algorithm developed for solving skin related layered sample problems such as those with dermis-fat or dermis-muscle structures. We carried out Monte Carlo simulations to obtain the benchmark frequency domain diffuse reflectance for layered samples with various layer thicknesses and optical properties. The simulation results were utilized to evaluate the performance of the algorithm and the layered diffusion model in the frequency domain.

2. Theoretical models

2.1 Layered diffusion model

Liemert and Kienle derived the solution of diffusion equation for N-layered turbid medium in frequency domain [11]. The two-layered and three-layered photon diffusion models used in this study were adopted from their work and are summarized in the following. In a layered turbid medium system, the time dependent diffusion equation can be written as:

The fluence rate at the air-sample interface for a layered sample whose bottom layer is infinitely thick has the following form in the Fourier domain:

2.2 Layered Monte Carlo model

The Monte Carlo model used here for calculating the benchmark frequency domain diffuse reflectance is an extension of the general multi-layer Monte Carlo code developed by Wang et al. [14]. In this paper, three kinds of sample setup will be discussed, namely, dermis over fat, dermis over muscle, and dermis over fat and muscle [3, 9, 15]. These three setups represent the most utilized sample structures for skin-related diffuse reflectance studies. Epidermis was neglected here because the SDSs used in this study were in the range from 5 mm to 20 mm and typical human epidermis of 50-200 μm thickness would have minimal influence on the diffuse reflectance measured at this SDS range. In the simulations, the optical properties of the dermis, fat, and muscle were set at (μ_a = 0.05 mm⁻¹, μ_s' = 1.5 mm⁻¹), (μ_a = 0.002 mm⁻¹, μ_s' = 1.0 mm⁻¹), and (μ_a = 0.05 mm⁻¹, μ_s' = 0.5 mm⁻¹), respectively, by referring the values reported by various literatures [2, 3, 9, 15]. For two-layered samples, the thickness of dermis layer was set between 0.5 to 3 mm, and fat or muscle were assumed to be infinitely thick (10⁸ mm). In the three-layered sample setup, the thickness of the dermis and fat were 1 mm and 5 mm, respectively. In all simulations, we employed a Henyey-Greenstein phase function with an asymmetry parameter of 0.8 and index refraction of 1.43 for all layers. Frequency domain diffuse reflectance, including amplitude demodulation and phase shift, were determined in the frequency range from 0 to 500 MHz with 1 MHz step using the Monte Carlo model, and all frequency data were used in the inversion fitting unless otherwise noted. All Monte Carlo simulations ran continuously until a total number of 50,000 photon packets arrived at the detector to achieve amplitude and phase noise less than 1.0% and 0.5°, respectively. For comparison, typical amplitude and phase measurement uncertainties have been reported to be less than 5.0% and 0.6°, respectively [13].

3. Results and discussion

In this section, we will first explore the outcome of directly applying a layered diffusion model for solving the five parameters of two-layered samples in 3.1. To reliably solve the two-layer problem, an iterative algorithm was developed and will be described in detail in 3.2. Furthermore, the iterative algorithm was extended to solve the three-layer problem and will be discussed in 3.3.

3.1 Direct optical properties recovery of two-layered samples with the layered diffusion model

To test the layered diffusion model’s capability, we fit the Monte Carlo simulated frequency domain diffuse reflectance, including amplitude demodulation and phase delay, of two-layered samples to the diffusion model to simultaneously recover the five parameters of samples. In the Monte Carlo simulations, the dermis-fat structure was employed with the top dermis layer thickness of 1 mm. The five parameters of two-layered samples recovered using the layered diffusion model from 16 independent Monte Carlo simulations of 16 SDSs ranging from 5 to 20 mm are shown in Fig. 1. In this study, the initial value of the least-squares routine (MATLAB “lsqcurvefit” function with the Trust-Region-Reflective algorithm) at all SDSs was set to the benchmark values of the sample unless otherwise noted. It can be seen in Fig. 1(a) and 1(b) that, in general, the recovery percent errors of the bottom layer optical properties (maximal recovery percent errors are 20% and 32% for μ_a2 and μ_s2', respectively) are smaller than those of the top layer optical properties (maximal recovery percent errors are 71% and 49% for μ_a1 and μ_s1', respectively). The recovery error of the top layer thickness can be as high as 177%. It is found that the bottom layer optical properties and the top layer reduced scattering coefficient could be more accurately recovered at large SDSs than at short SDSs. Among the five parameters, the top layer absorption coefficient and the top layer thickness both have higher recovery errors than other parameters in general.

 

Fig. 1 Recovered (a) absorption coefficients and (b) reduced scattering coefficients of the top layer (squares) and bottom layer (triangles) at various source-detector separations. (c) Recovered top layer thickness L₁ (asterisks) at various source-detector separations. Solid lines represent benchmark values.

Download Full Size | PDF

The fact that the optical properties of a two-layered sample cannot be reliably recovered using the frequency domain two-layered diffusion model at a single SDS can be further understood by Fig. 2. To generate Fig. 2(a), we first ran a Monte Carlo simulation to calculate the frequency domain reflectance of a sample with dermis-fat structure (Y) using the following parameters: μ_a1 = 0.05 mm⁻¹, μ_s1' = 1.5 mm⁻¹, μ_a2 = 0.002 mm⁻¹, μ_s2' = 1.0 mm⁻¹, n = 1.43, top layer thickness = 1 mm, modulation frequency = 0 to 500 MHz with 1 MHz step, and the SDS = 6 mm. Next, by keeping the parameters other than μ_a1 and μ_s1' the same as those used in the Monte Carlo simulation, we utilized the two-layered diffusion model to calculate the frequency domain reflectance at various μ_a1 and μ_s1' combinations (Y_j). Here the subscript j is the number of the μ_a1 and μ_s1' combinations, and it was set in the range from 1 to 2500 for plotting Fig. 2. Finally, we compute χ² = (Y –Y_j)² at various μ_a1 and μ_s1' combinations and plot Fig. 2(a). In computing χ² displayed in Fig. 2 and 3, the amplitude magnitudes were magnified by the ratio of average phase to average amplitude, and we summed up the results at all modulation frequencies. It can be clearly seen in Fig. 2(a) that as there are only two variables to be fit (μ_a1 and μ_s1'), a least-squares fitting method can lead to a single solution for this two-layered sample problem since there is only one global minimum. It should be noted that the trough of the curve in Fig. 2(a) corresponds to (μ_a1 = 0.0449 mm⁻¹, μ_s1' = 1.724 mm⁻¹) which deviates from the benchmark top layer optical properties set in the Monte Carlo simulation. This suggests that there exists system differences between Monte Carlo method and the diffusion equation as indicated by Kienle et al. [12] and the diffusion equation has non-negligible errors even in the diffusion regime. Following similar procedure in creating Fig. 2(a), we generated Fig. 2(b) by changing top layer thickness to 1.5 mm and keeping other parameters in the diffusion model the same as those used for generating Fig. 2(a). Note that we employed the same benchmark Monte Carlo simulation data (Y) for generating Fig. 2(a) and 2(b). From Fig. 2(b), we can find that the value of (Y_i -Y_i)² at the global minimum is 0.4272 which is smaller than that of Fig. 2(a) (0.4558). This result implies that as the free parameters increase to 3 (μ_a1, μ_s1', and top layer thickness), a least-squares fitting routine may find a solution that greatly deviates from the benchmark values for this two-layer problem. In addition, we found that increasing the number of free parameters to 4 or 5 would lead to multiple local minima in the χ² space for this two-layered sample problem (data not shown) and thus a least-squares fitting routine employed for solving this problem would find different answers for different initial value settings.

 

Fig. 2 χ² computed at various top layer μ_a1 and μ_s1' combinations for (a) top layer thickness = 1 mm, and (b) top layer thickness = 1.5 mm, at source detector separation of 6 mm.

Download Full Size | PDF

 

Fig. 3 χ² computed at various top layer μ_a1 and μ_s1' combinations for (a) top layer thickness = 1 mm, and (b) top layer thickness = 1.5 mm, at source detector separation of 13 mm.

Download Full Size | PDF

We also investigated the effect of simultaneously employing multiple SDSs on the accuracy of the recovered optical properties of a dermis-fat structured sample. We tested various combinations of SDSs and part of the results are listed in Table 1. In Table 1, data recovery results listed were obtained with 2, 3, and 4 source-detector pairs. All three combinations include the relatively short (6 mm) and long (20 mm) SDSs, and two of them also include intermediate SDSs. It can be seen that μ_s1', μ_a2 and μ_s2' have recovery errors less than 10%, while the recovery errors for μ_a1 and top layer thickness are greater than 25.2% and 192.6%, respectively. In general, employing multiple source-detector pairs produced better recovery results than the single source-detector pair counterpart. For reference, the errors of (μ_a1, μ_s1', μ_a2, μ_s2', L₁) recovered at 6 mm SDS were −81.0%, 46.5%, 5.0%, −39.4%, and 33.3%, respectively. However, the benefit of increasing the source-detector pair number to more than 2 is found only for μ_s1' and μ_s2' but is not evident for other parameters especially the top layer thickness. Our results indicate that the two-layered sample problem cannot be reliably solved by simultaneously employing multiple source-detector pairs.

Table 1. Optical properties of top (μ_a1, μ_s1') and bottom (μ_a2, μ_s2') layers as well as top layer thickness (L₁) recovered at various source-detector separation combinations for a dermis-fat structured sample. Recovery errors are listed below the values.

View Table | View all tables in this article

3.2 Optical properties recovery of a two-layered sample by iteratively using a layered diffusion model at multiple SDSs

3.2.1 Iterative algorithm and selection of SDSs

By carefully reading the data shown in Fig. 1 and Table 1, we found that the bottom layer properties could be properly determined when the long SDS data was included in data fitting. This is reasonable since at larger SDSs, photons traveling longer path are more likely to be detected and these photons are more sensitive to the variation of bottom layer properties and less sensitive to the properties of a thin top layer. In the limiting case, the two-layered sample problem with a thin top layer approaches a homogeneous semi-infinite sample problem as the SDS becomes large and this problem can be unambiguously solved with the frequency domain diffuse reflectance measured at a single SDS [5]. Nevertheless, short SDS data cannot be employed to reliably determine the top layer properties when the bottom layer properties are also fitting parameters because simultaneously determining five variables of a two-layered sample is an ill-posed problem as explained in the previous section.

As investigating the two-layered sample problem, we discovered an interesting fact that the top layer thickness could be better recovered at some SDSs than others. This fact can be seen in Fig. 1(c) where the top layer thickness can be more accurately recovered at intermediate SDSs roughly ranging from 8 to 13 mm than at other SDSs. Besides, in this intermediate range, top layer absorption and scattering properties can sometimes be properly recovered. We generated plots in a similar way to those demonstrated in Fig. 2 to understand this phenomenon. The reference data was generated from our Monte Carlo model to simulate the dermis-fat structure with top layer thickness = 1 mm and SDS = 13 mm. We then calculated the χ² between the reference data and the reflectance determined from the two-layered diffusion model at various (μ_a1, μ_s1') combinations to generate Fig. 3. The top layer thickness in the diffusion model was set at 1mm for Fig. 3(a) and 1.5 mm for Fig. 3(b). It can be seen in Fig. 3 that χ² is smaller at 1 mm than at 1.5 mm which is in contrary to the trend demonstrated in Fig. 2. In fact, as the bottom layer optical properties were fixed at the benchmark values and there were three floating variables, the minimum of the χ² for SDS = 13 mm was found at L₁ = 1.028 mm, μ_a1 = 0.057 mm⁻¹, and μ_s1' = 1.89 mm⁻¹. This implied that top layer properties could be properly estimated at intermediate SDSs as the bottom layer optical properties were known.

Furthermore, it can be noticed in Fig. 2(a) that when the top layer thickness is close to 1 mm and the SDS = 6 mm, the best top layer optical properties that minimize χ² are near μ_a1 = 0.045 mm⁻¹, μ_s1' = 1.72 mm⁻¹, which are closer to the benchmark top layer optical properties than those determined at SDS = 13 mm (μ_a1 = 0.057 mm⁻¹, μ_s1' = 1.89 mm⁻¹). This suggests that if the top layer thickness is fixed at a value that is close to the true value, 6 mm SDS measurement data could lead to a better top layer optical properties recovery than the 13 mm counterpart. It can thus be inferred that measurements conducted at short and intermediate SDSs are more sensitive to the variations of top layer optical properties and top layer thickness, respectively, than those conducted at other SDSs.

Based on the facts discussed above, we developed an algorithm in an effort to robustly and accurately determine the optical properties of a two-layered sample from frequency domain diffuse reflectance. The essence of the algorithm is to employ data measured at long, intermediate, and short SDSs to provide the reference values of (μ_a2, μ_s2'), L₁, and (μ_a1, μ_s1'), respectively. New reference values are calculated by referring to previous reference values, and such iteration keeps running until the convergence condition of the five parameters is met. The configuration of the algorithm is described in detail in the following. The first step in our algorithm for determining the five parameters of a dermis-fat sample is to calculate the reference bottom layer optical properties at 20 mm SDS. (The other three parameters are also floating variables in this step.) Second, the top layer optical properties and the top layer thickness are calculated at the intermediate 13 mm SDS. It should be noted that in our program, the bottom layer optical properties are fixed at the reference values obtained in the first step here. Third, similar to the second step, the three parameters of the top layer are calculated at the 6 mm SDS with the bottom layer optical properties fixed at the values calculated in the first step and the top layer thickness and optical properties are allowed only to vary within 5% and 10%, respectively, of those calculated from the second step. Larger variation percentages for the top layer properties could be used in the third step to reduce the number of iterations but this may also lead to an instable system which produces values that diverge far from true values. Finally, in the fourth step, the bottom layer optical properties are calculated again at 20 mm SDS with the three parameters of the top layer fixed at the values calculated in the third step and then updated. This finishes the first round of iteration, and starting from the second round of iteration only step 2 to step 4 described above are carried out until stable results are obtained. A flowchart including all key steps in the iterative algorithm is illustrated in Fig. 4.

 

Fig. 4 Flowchart of the iterative algorithm for solving the two-layered sample problem.

Download Full Size | PDF

The algorithm described above was used to recover the properties of a dermis-fat sample with 1 mm top layer thickness, and the results are shown in Fig. 5. Three distances 6, 13, and 20 mm were chosen to cover the short, intermediate, and long SDS regimes, respectively. We can see that the bottom layer optical properties were quite stable throughout the iterative calculation process, while the other three parameters gradually approached the values that were near the benchmark values. The iteration algorithm was terminated under the following conditions: 1. One of the determined top layer properties of this iteration had less than 0.1% variation from those of the last iteration. (convergence termination) 2. One of the determined top layer properties showed a reverse trend in the iteration and the variation was larger than 1%. (oscillation/divergence termination) Termination condition 1 usually implies successful iteration process is achieved while termination condition 2 suggests that the calculated values are possibly oscillating or starting to diverge. When the variation percentage setting is judiciously adjusted in the step 3 of an iteration described above, termination condition 2 is satisfied most likely due to oscillating results rather than diverging results. To avoid the algorithm from premature termination, the algorithm was set to have a minimal iteration number of 5. The iteration process for obtaining the results demonstrated in Fig. 5 was terminated because the determined top layer thickness satisfied the termination condition 2. It is worth noting that the recovered values of the top layer properties in the first round of iteration were not close to the benchmark values which suggested that our algorithm was robust. In addition, we also found that the algorithm was not sensitive to the initial value settings.

 

Fig. 5 Recovered (a) absorption coefficients and (b) reduced scattering coefficients of the top layer (squares) and bottom layer (triangles), as well as (c) recovered top layer thickness L₁ (asterisks) versus iteration number for a dermis-fat structure. Solid lines represent benchmark values.

Download Full Size | PDF

The final recovery results of Fig. 5 are listed in the first row of Table 2. To understand the effect of measurement noise on the recovered results, we superimposed additional 5% amplitude and 0.6° phase noises to the Monte Carlo data to simulate the measurement noise of real frequency domain reflectance systems [13]. The percent errors of the recovered μ_a1, μ_s1', μ_a2, μ_s2', and L₁ were −11.5%, 15.7%, 5.0%, −4.0%, and 11.2%, respectively, which are comparable to the values listed in the second row of Table 2. This result indicates that the iterative algorithm is still robust under typical measurement noise.

Table 2. Optical properties of top (μ_a1, μ_s1') and bottom (μ_a2, μ_s2') layers as well as top layer thickness (L₁) recovered at various source-detector separation combinations for a dermis-fat structured sample. Recovery errors are listed below the values.

View Table | View all tables in this article

Moreover, we investigated the influence of the choice of SDSs on the recovered results. The sample setup here was again a dermis-fat structure with 1 mm top layer thickness. In making the combination of the SDSs, we kept the longest SDS fixed at 20 mm and adjusted the shortest SDS. The intermediate SDS was approximately set at a position which was close to the average of the longest and the shortest SDSs. The results are listed in Table 2 and it can be seen that all recovered values are similar as the shortest SDS varies from 6 to 8 mm, and the (8, 14, 20) mm combination produced the best top layer results among them. As the shortest SDS increases to 12 and 14 mm, the errors of the recovered top layer properties grow accordingly. This means that the shortest SDS has to be short enough to enhance the accuracy of the top layer properties recovery. It is worth noting that the bottom layer optical properties were properly recovered for all combinations with 20 mm SDS. In the last set of data shown in Table 2, we reduced the longest and intermediate SDSs. The recovery errors for the bottom layer properties of this data set were larger than those of the other data sets. Larger bottom layer optical properties recovery errors as well as a not optimized intermediate SDS setting led to the great deviation of the recovered top layer properties from the benchmark values in the last data set.

To further understand that whether the sample’s optical properties would affect the performance of the iterative algorithm, we applied the algorithm to calculate the properties of another commonly seen two-layered skin setup, the dermis-muscle structure. The Monte Carlo reference data was generated using following parameters: μ_a1 = 0.05 mm⁻¹, μ_s1' = 1.5 mm⁻¹, μ_a2 = 0.05 mm⁻¹, μ_s2' = 0.5 mm⁻¹, and top layer thickness = 1 mm. We employed the (6, 13, 20) mm SDS combination in the algorithm. In Fig. 6(c), at the first iteration, the percent error of the recovered top layer thickness is 95.4%, and we can observe that the top layer thickness approaches the benchmark value as the iteration goes. This iteration process was terminated because the variation of the recovered top layer thickness met the convergence termination condition described earlier. The percent errors of the final recovered μ_a1, μ_s1', μ_a2, μ_s2', and L₁ were reduced to −13.2%, 13.7%, −5.4%, −2.0%, and 23.9%, respectively. Under a similar setup, we tested our algorithm on two-layered skin structures with various optical property settings (μ_a1 = 0.01-0.1 mm⁻¹, μ_s1' = 1.0-3.0 mm⁻¹, μ_a2 = 0.002-0.1 mm⁻¹, μ_s2' = 0.5-1.5 mm⁻¹) and all demonstrated similar converging trends as those shown in Fig. 5 and 6.

 

Fig. 6 Recovered (a) absorption coefficients, (b) reduced scattering coefficients, and (c) layer thickness of the top layer versus iteration number for a dermis-muscle structure. Solid lines represent benchmark values.

Download Full Size | PDF

3.2.2 Influence of the layer thickness

The thickness of the human skin dermis is known to be site and subject dependent and is typically in the range from 0.5 to 3 mm. To test the capability of our algorithm in recovering the dermis thickness of the dermis-fat and dermis-muscle structures, in the Monte Carlo simulations, the top layer thickness was varied from 0.5 to 3 mm in a 0.5 mm step. It should be noted that as the transport mean free path is larger than the top layer thickness, alternative solutions for the two-layered diffusion model have to be employed [16]. The SDS combination used in the algorithm was (6, 13, 20) mm. The recovered top layer optical properties and thickness for the dermis-fat sample structure are shown in Fig. 7(a), 7(b), and 7(c). It can be seen in Fig. 7(a) and 6(b) that the determined top layer optical properties have weak dependence on the top layer thickness, and in Fig. 7(c) the recovered top layer thicknesses and the benchmark values are very close and the slope of the linear regression line of the data is 0.961. Similar trends can be observed for the dermis-muscle sample structure as shown in Fig. 8(a), 8(b), and 8(c), and the slope of the regression line of the recovered top layer thicknesses is 0.959. Results shown in Fig. 7(c) and 8(c) indicate that the iterative algorithm can be used to determine the top layer thickness of a two layered sample with errors within 19% if the L₁ = 0.5 mm data is excluded and within 39% if included. Similarly, the recovered top layer optical properties shown in Fig. 7(a), 7(b), 8(a), and 8(b) have percent errors within 20% if the L₁ = 0.5 mm data is excluded and within 36% if included. We found that the deviation of the reflectance determined using the layered diffusion model from the benchmark values at the SDS of 6 mm increased as top layer thickness decreased from 1 to 0.5 mm. An alternative photon transport model that is more accurate in determining the reflectance at short SDSs could be employed to work with the iterative algorithm to improve the results when top layer thickness is smaller than 1 mm.

 

Fig. 7 Recovered (a) absorption coefficients, (b) reduced scattering coefficients, and (c) layer thickness of the top layer at various top layer thickness for a dermis-fat structure. Solid lines represent benchmark values.

Download Full Size | PDF

 

Fig. 8 Recovered (a) absorption coefficients, (b) reduced scattering coefficients, and (c) layer thickness of the top layer versus top layer thickness for a dermis-muscle structure. Solid lines represent benchmark values.

Download Full Size | PDF

3.2.3 Hemodynamics

Several researchers studied the muscle optical properties as well as muscle hemodynamics using a semi-infinite diffusion model and reported compromised accuracy of the recovered muscle optical properties due to the presence of dermis or fat [3, 15]. This problem can be approximated to a dermis-muscle structure problem and solved with a two-layered diffusion model. We investigated that whether the bottom layer absorption coefficient recovered using the iterative algorithm can faithfully reflect the variation of the muscle layer absorption property. To this end, in the Monte Carlo simulations, the bottom layer absorption coefficient was varied from 0.05 to 0.07 mm⁻¹ to simulate muscle hemodynamics and the other parameters were kept fixed. The layer thickness of dermis was 1 mm. The recovery results are shown in Fig. 9. In Fig. 9(a), the recovered muscle absorption coefficients are systematically lower than the benchmark values by 0.0039 mm⁻¹ in average (6% error in average) and the slope of the regression line of the two sets of data is 1.069. Our results suggest that muscle hemodynamics could be tracked using our method with a modest error. The other four parameters plotted in Fig. 9(a) to 9(c) do not show clear dependence with the bottom layer absorption coefficient.

 

Fig. 9 Recovered (a) absorption coefficients and (b) reduced scattering coefficients of the top layer (squares) and bottom layer (triangles), as well as (c) recovered top layer thickness L₁ (asterisks) versus top layer thickness for a dermis-muscle structure. Sample optical properties recovered using a semi-infinite diffusion model are depicted as empty circles. Solid lines represent benchmark values.

Download Full Size | PDF

In addition, using the same set of two-layered Monte Carlo data, we calculated the optical properties of sample using a frequency domain semi-infinite diffusion model at SDS of 20 mm, and the results are shown as empty circles in Fig. 9(a) and 9(b). It can be seen that although the recovered absorption coefficients recovered with a semi-infinite diffusion model vary in accordance with the benchmark values, they remarkably deviate from the benchmark values by 0.023 mm⁻¹.On the other hand, the reduced scattering coefficients recovered using the semi-infinite model do not show strong dependence on the variation of bottom layer absorption, but they are less accurate than those recovered using the iterative algorithm.

3.3 Optical properties recovery of a three-layered sample by iteratively using a layered diffusion model at several SDSs

3.3.1 Iterative algorithm and selection of SDSs

The skin of trunk and limbs usually sits on a fat-muscle structure which forms a three-layered sample model. The thickness of fat layer is highly body-site, gender, and subject dependent, and can vary between 1 to 70 mm [17]. If the fat layer is very thin, the three layer problem can be approximated as a dermis-muscle two-layer problem. On the other hand, if the fat layer is extremely thick, the three-layer problem can be treated as a dermis-fat two-layer problem. When the fat layer has a moderate thickness, the three-layer problem is practical that one needs to consider when studying muscle hemodynamics using DRS techniques.

As shown previously, the two-layer problem can be solved by choosing an appropriate combination of SDSs. It is reasonable to expect that by extending the iterative algorithm developed for the two-layer problem, the three-layer problem could be solved. The concept of choosing different SDSs to determine the reference values of sample parameters of different layers mentioned earlier is still used for solving the three-layer problem but slightly modified. In a three-layered sample problem, there are eight parameters to be determined including (μ_a1, μ_s1', L₁) for the top layer, (μ_a2, μ_s2', L₂) for the middle layer, and (μ_a3, μ_s3') for the bottom layer. We utilize a long SDS to calculate the reference values of (μ_a3, μ_s3'), two intermediate SDSs to calculate the reference values of (μ_a2, μ_s2', L₂), and two short SDSs to calculate the reference values of (μ_a1, μ_s1', L₁). For instance, our algorithm consists of following four major steps when the SDS combination is (6, 10, 13, 17, 20) mm: 1. The reference bottom layer optical properties as well as the other six parameters are first calculated at 20 mm SDS. 2. The top layer optical properties and the top layer thickness are calculated at 6 and 10 mm SDS with the middle and bottom layer optical properties fixed at the values calculated in the previous step. In this step, 10 mm data is used to calculate the top layer thickness and optical properties are allowed to vary within 10% and 5% of those calculated from the first step. The updated results are then passed to the 6 mm data processing where the top layer thickness and optical properties are allowed to vary within 5% and 10% of the updated reference values. 3. Similar to the second step, the three parameters of the middle layer are calculated at the 13 and 17 mm with the top and bottom layer properties fixed. 4. The bottom layer optical properties are calculated again at 20 mm SDS with the other six parameters fixed at the values calculated in the second and the third steps. After the first round of iteration is completed, the second and further iterations only need to carry out steps from 2 to 4. The iteration of algorithm stops when the two termination conditions described in the previous section are satisfied. We tested the algorithm on a Monte Carlo simulated dermis-fat-muscle data. The thicknesses of dermis and fat layers were 1 and 5 mm, respectively. The recovery results and their percent errors are listed in Table 3. It can be seen that the bottom layer optical properties were less accurately recovered than the other parameters. By looking up the recorded photon trajectories, we found that this phenomenon could be caused by the fact that some photons arriving at the 20 mm SDS detector did not ever enter the bottom layer and thus did not carry the bottom layer information. This situation became worse when the middle layer thickness increased and could be alleviated by increasing the distance of the longest SDS pair. However, measurement noise increases as the SDS increases. The attainable longest SDS is system dependent and is typically less than 30 mm [1, 15].

Table 3. Recovered properties of top (μ_a1, μ_s1', L₁), middle (μ_a2, μ_s2', L₂), and bottom (μ_a3, μ_s3') layers for a dermis-fat-muscle structured sample using the iterative algorithm. Recovery errors are listed below the values.

View Table | View all tables in this article

The SDS combination of (6, 10, 13, 17, 20) mm does not work for all kinds of three-layered sample structures such as the structure with thick middle layer. For such cases, larger intermediate SDSs may be required to obtain reasonable results for the middle layer. For example, we ran another set of Monte Carlo simulation where the middle fat layer thickness was increased to 10 mm and kept other parameter values untouched. By using the SDS combination of (6, 10, 13, 17, 20) mm, the recovered (μ_a3, μ_s3') had percent errors of −62.7% and 93.9%, respectively. When the SDS combination was adjusted to (6, 10, 20, 25, 35) mm, the errors of recovered (μ_a3, μ_s3') were reduced to −36.0% and −34.1%, respectively. Therefore, picking a proper combination of SDS based on the structure of the site under investigation is the prerequisite for obtaining accurate properties of three-layer samples.

3.3.2 Hemodynamics

As mentioned earlier that the three-layered diffusion model has its significance in studying muscle hemodynamics. Thus, it is crucial to understand that whether the iterative algorithm with a three-layered diffusion model can be used to track the hemodynamics of muscle tissue in a three-layered tissue structure. Monte Carlo simulations were carried out to generate the experimental data. In the Monte Carlo simulations, the thicknesses of dermis and fat layers were 1 and 5 mm, respectively, and the absorption coefficient of muscle layer was varied from 0.05 to 0.07 mm⁻¹. The recovered muscle layer absorption and reduced scattering coefficients using the SDS combination of (6, 10, 13, 17, 20) mm are shown in Fig. 10 as filled triangles. In Fig. 10(a), the recovered muscle absorption coefficients increase with the benchmark values and the linear regression slope is 0.538. Also depicted in Fig. 10 as empty circles are the recovered muscle optical properties using the semi-infinite diffusion model at 20 mm SDS. It can be seen that the muscle optical properties determined using the iterative algorithm are closer to the benchmark values than those determined using the semi-infinite model. The absorption coefficients calculated using the semi-infinite model barely follow the variation trend of the benchmark values and have a linear regression slope of 0.145. Our simulation results suggest that the semi-infinite model is not suitable for tracking muscle hemodynamics in a three-layered tissue model.

 

Fig. 10 Recovered (a) absorption coefficients and (b) reduced scattering coefficients of the muscle layer (triangles) at various muscle layer absorption for a dermis-fat-muscle structure. Crosses are the results obtained with fixed top and middle layer thicknesses. Sample optical properties recovered using a semi-infinite diffusion model are depicted as empty circles. Solid lines represent benchmark values.

Download Full Size | PDF

Although the absorption coefficients calculated using the three-layered model are closer to the benchmark values than those calculated using the semi-infinite model, their linear regression slope is only 0.538. This could be the result of the masking effect of the middle layer mentioned in the previous sub-section and could be alleviated by increasing the SDS. We carefully analyzed the data and found that the recovered muscle absorption coefficients were affected by the small fluctuations of the recovered top layer and middle layer thicknesses. The top layer and middle layer thicknesses recovered from the five independent Monte Carlo simulated experiments varied in the range from 1.0278 to 1.0559 mm and from 4.7356 to 5.0324 mm, respectively. By keeping the top layer and the middle layer thicknesses at the values recovered from the first experiment (L₁ = 1.0448 mm, L₂ = 4.9716 mm) and then performing the data recovery for the rest of experiments (with 6 floating parameters), we obtained bottom layer optical properties as those illustrated in crosses in Fig. 10. The linear regression slope of the crosses in Fig. 10(a) is 1.262, which means the sensitivity to the bottom layer absorption variation is ameliorated by locking the top and the middle layer thicknesses. For most clinical applications, it is reasonable to keep the top layer and middle layer thicknesses fixed at the values determined at a certain time point and then track the time evolution of optical properties. Our algorithm with a three-layered diffusion model provides a robust and sensitive way for tracking the hemodynamics of muscle tissue.

4. Conclusion

We observed that the measurements conducted at short and long SDSs were more suitable for recovering the optical properties of top and bottom layers of a layered sample, respectively. Based on our observation, we developed an algorithm in an attempt to accurately determine the properties of layered samples. The iteratively algorithm revealed in this study is capable of solving the optical properties and layer thicknesses of layered samples. We tested this algorithm using Monte Carlo simulations and found that it was robust and could reliably recover the optical properties and layer thicknesses of s two-layered or three-layered sample structures. For a two-layered dermis-fat or dermis-muscle structure with maximum top layer (dermis) thickness of 3 mm, we found that (6, 13, 20) mm was a good SDS combination for the iterative algorithm. For all the two-layered skin cases we tested, the recovered optical properties and layer thickness had percent errors less than 25%. On the other hand, for the three-layered skin structure, it is difficult to determine the best SDS combination since the thickness of the middle fat layer is highly site and subject dependent and has a great variation. For example, as the top and middle layer thicknesses were 1 and 5 mm, respectively, we found that the reflectance data measured at a SDS combination of (6, 10, 13, 17, 20) mm was sufficient for accurately deriving sample optical properties and thicknesses. However, when the middle layer thickness was increased to 10mm, the longest SDS had to be increased to at least 35 mm to obtain reasonable bottom layer optical properties.

In this study, it typical required 1000-3000s to obtain the final results on an Intel i7-4930K based computer. We found that the computation time could be reduced to less than 60s by only employing 11 frequencies in the 0-500 MHz range; however, the recovery accuracy was slightly compromised. Due to the speed limitation, this algorithm may not be suitable for clinical settings that require real-time information of tissues.

We employed a frequency domain layered diffusion model to work with the iterative algorithm in this study; nevertheless, it is reasonable to expect that other photon transport models in frequency or time domain can also work with the algorithm. For example, photon diffusion models are generally known to be incapable of describing the photon transport in highly absorbing materials or at short SDSs, while Monte Carlo methods do not have such limitations. Employing Monte Carlo models together with the iterative algorithm could help determine the optical properties of highly absorbing layered samples or the thickness of samples with very thin layers. Furthermore, the algorithm could be extended to study samples that have layer number larger than three such as human or animal head models.