Analysis of diffuse reflectance spectroscopy by means of Bayesian inference and separation of the parameters for scattering strength and spectral dependence of the scattering

For many applications in tissue optics, knowledge of the scattering coefficient or at least the reduced scattering coefficient is essential. In addition, its spectral dependence is also an important feature, as this provides information about scatterer sizes. While the characterization of the spectral dependence should be a simple fit, experimental results show strong fluctuations even for the same tissue type. For in‐vivo measurements, this problem is even greater. Therefore, it is the aim of this study to analyze the instabilities of the scattering characterizations and find a solution for it by means of Bayesian inference. In this study, this behavior is investigated using the example of diffuse reflectance spectroscopy. It can be shown that the currently used fitting functions are unstable for the fitting as both parameters for characterizing the reduced scattering coefficient describe the spectral dependence as well as the scattering strength. However, by a simple coordinate transform, a stable mathematical description of the scattering is derived. By the fact that a posteriori probability of the reduced scattering coefficient narrows down significantly with the Bayesian inference, the new fitting function is verified.


| INTRODUCTION
For many applications in tissue optics, knowledge of the scattering coefficient or at least the reduced scattering coefficient is essential. For some diseases, scattering might be used as a disease-related parameter. In addition, knowledge of the scattering is essential for manufacturing phantoms to test different optical imaging modalities for calibration. Also, the design of completely new approaches such as shifted position-diffuse reflectance imaging [1] requires reliable knowledge of the optical properties of tissues.
Nevertheless, the measured optical properties for tissues vary greatly for the same tissue type [2]. The strong variations are to a large extent caused by the inherent variations of tissues from different individuals. However, the differences are also caused by different approaches or measurement strategies from different groups. This problem is discussed in detail for the characterization of intralipid (IL) as one of the most important optical phantoms [3,4]. Based on previous work, the optical properties of IL could be characterized with an integration sphere setup with very high precision without dependent scattering [5] and with dependent scattering [6]. Despite the high accuracy of their work, there are still some minor errors in the final optical properties, especially for scattering [7], which is discussed in the next paragraph.
As they use a fitting formula to describe the whole range of the optical properties. For low concentrations, it often happens that the description is not so well due to the minimization of the root mean square error by the fitting. This effect is expected to happen also in the studies from Aernouts et al. [5,6]. It should also be noted that the scattering coefficient of IL changes with temperature [8]. However, this temperature dependence is lower than the difference of the scattering found in our former study [7]. Moreover, the properties should not change much with age ("Measurements repeated on samples from the same batch showed that there are no appreciable variations over a period of about 2 years, but the small batchto-batch variations we observed suggest that the optical properties remain probably stable over a much longer time (the oldest batch on which we have carried out measurements has been manufactured about eleven years ago).") [3]. This also shows that batch-to-batch variations should be low [3].
In summary, the problem of different results from different groups is still existent. One possible way of reducing this problem is averaging the results by fitting. Typically, the spectral dependence of the reduced scattering coefficient is fitted [2] by the following equation: where μ 0 s is the reduced scattering coefficient, λ is the wavelength and the fitting constants are a À and b À . Alternatively, it is fitted by [2]: In Equation (2), the scattering is split into a fraction of Rayleigh scattering and a fraction of Mie scattering, described by f Ray . The fitting constant b À is therefore replaced by the exponent from the Mie scattering b Mie and a À is replaced by a 0 . Equation (2) has the advantage that Mie scattering can be described more in detail as well as that the fraction of the Mie scattering can be derived. However, the fitting approaches also have disadvantages. The fitting function does not seem to be very stable. This can be seen from the fact that the constants in Equations (1) and (2) vary strongly even for the same tissue type [2] (Table 1). Therefore, it is possible or even likely that the fitting parameters do not represent the physically correct parameters. Furthermore, the arbitrary normalization to 500 nm causes artifacts [9] (fig. 5).
In general, the instability of the fits results in larger issues when these functions are inserted in a diffuse reflection spectroscopy (DRS) model. This can be partly compensated by changing Equation (1) to the following: where a is the new fitting constant. It should be noted that the unit is not strictly correct. In this equation, a cannot have a unit. For this reason, the unit is multiplied. Furthermore, there are strong hints in our previous study [9] that the introduction of Bayesian inference allows the representation of the instability of the fits. In our former study [9], the model from Perelmann et al. [10] was used. However, the presented function is rather complicated and shows systematic errors. Later extensions even complicate the description of the diffuse reflectance more [11]. Hence, they are not suited for a stability analysis.
There are further models relying on diffusion theory [12], which only works for the diffusion approximation. Moreover, these models require certain source-detector distances and a priori known parameters [12]. Thus, the results are seen as not generalizable enough. A spatially resolved approach can be realized by absolute measurements of the diffuse reflectance as input to trained neural network (NN) to estimate the optical properties [13,14]. However, a black-box approach such as NN is not suited for stability analysis. Empirical modeling is another potential approach realized by Bargo et al. [15]. This approach works. However, formal expression is even more complicated than the Perelmann et al. [10] approach.
The following solution was found by Zonios et al. [16] allowing an easy description of the diffuse reflectance with a precise reconstruction of the optical properties: where R λ ð Þ is the reflectance, μ a λ ð Þ is the absorption coefficient and k 1 and k 2 are parameters describing the geometry of the experimental probe.
In this study, we use Bayesian inference as a tool for a stability analysis of diffuse reflectance spectra. Thereby, a new formula for the reduced scattering coefficient is presented which is stable and has only one dominant variable, which describes universally the scattering strength independent of the presence of Rayleigh or Mie scattering. Figure 1 shows the overview of the methodology of this study. First, the diffuse reflectance samples are measured. The samples are used to calibrate the model from Zonios et al. [16]. Afterward, the spectrum is fitted by means of Bayesian inference for both models (the standard model and model with different coordinate system). The scatter plots from the Bayesian inference are generated and afterward the a posteriori probabilities of the scattering and absorption coefficient are shown exemplarily.

| METHODS
Moreover as in this study, the focus is on the stability description, the equations do not have the same units on the left and right side. However in the whole study for the (reduced) scattering coefficient as well as the absorption coefficient, the unit is 1 cm . To improve the readability, the additional naming of the unit as in Equation 3 is omitted.

| Experimental parameters and set-up
Set-up Figure 2 shows the DRS set-up. The main set-up consists of a 50μm bifurcated fiber (2 m, QBIF50-VIS-NIR, Ocean Optics, Dunedin, USA). One end is used for illumination. Two planoconvex lenses with a focal length of 50 mm are used to adjust the light beam from the light source. The adjustment is performed to reach maximal light intensity. The other side of the bifurcated fiber is connected to a spectrometer (QE65000 Spectrometer, Ocean Optics, Dunedin, USA). For all measurements, the fiber is placed inside the phantom to minimize effects from imperfections of the alignment such as distance surface to probe, angle probe to surface. A BaSO 4 plate is used as a calibration sample. The final reflectance is calculated as follows: where R λ ð Þ meas is the measured spectrum for the further analysis, M λ ð Þ is the measured raw signal, D λ ð Þ is the dark signal and R λ ð Þ BaSO 4 is the measured signal from the BaSO 4 plate.

Samples
As sample, different concentrations of MontBlanc burgundy red ink (MontBlanc, Hamburg, Germany) mixed with intra lipid (Intra lipid 20%, Fresenius Kabi, Germany) are used. The used ink concentrations are 0.25%, 0.5%, 1.2%, 1.7% and 2.4%. The extinction coefficient shown in Figure 3 is measured by a spectrophotometer (UV3600, Shimadzu, Japan). The extinction coefficient is described as follows: where c ink is the ink concentration and ε ink is the measured extinction coefficient of the ink. The intra lipid concentrations are 1%, 2%, 4% and 5%. Thereby, the concentration describes the fat concentration and the IL has a stock concentration of 22.67%. The reduced scattering coefficient is taken from Aernouts et al. [5]: where c IL is the concentration of IL.

| Calibration of the model
The model from Zonios et al. [16] has to be calibrated to known samples to derive the parameters k 1 and k 2 .
For this, in Equations (4), (6) and (7) are inserted with the known ink and IL concentrations: where R λ, k 1 , k 2 ð Þ c IL ,c ink is the reflection for a known sample to calibrate the geometry. When Equation (8) is fitted, the calibration model from Zonios et al. [16] is performed. This is done for all 20 combinations of c IL and c ink .
The constants for k 1 and k 2 are calculated by standard nonlinear least square fitting. The final values for k 1 and k 2 are chosen by the geometric mean. This is necessary because k 1 and k 2 depend on the optical properties. In this case, k 1 and k 2 are considered as fixed priori information. This contradiction that k 1 and k 2 are not completely independent of the optical properties is the weakness of model from Zonios et al. [16]. Despite this The main part of the set-up is the bifurcated fiber. The light from the light source is guided in one fiber while the other end is connected to the spectrometer. The tip of the probe is inside the sample F I G U R E 1 Flowchart showing the overall methods F I G U R E 3 Extinction coefficient of MontBlanc burgundy red ink weakness, the model from Zonios et al. [16] still provides reasonable good results.

| Diffuse reflectance model
In general, the absorption and reduced scattering can be described as follows: where a and b are the fitting constants for the reduced scattering coefficient. The diffuse reflectance is described by the model from Zonios et al. [16] shown in Equation (4) resulting in the following final equation for describing the diffuse reflectance with Equations (6) and (9): Since the fit tends to be unstable, a different representation of the reduced scattering coefficient is used as well. It is known that within the diffusion approximation the total scattering strength is described exactly. However, the phase function does not play a significant role any more. Therefore, the description of the scattering strength should work well, while effects from the phase function cannot be described without major errors. Therefore, the variation in the Bayesian results introduced by the scattering strength should be small, while the variation introduced by the phase function should show a strong variance in the Bayesian results. Thus, a plot of the parameters a vs b could lead to a certain direction of variance. Hence, the aim of the proposed coordinate transformation is similar to that of a principal component analysis. Rotate the direction of greatest variance in one direction: where the coordinates of a ϕ and b ϕ are rotated by the angle ϕ. This results in the following equation for the diffuse reflectance: For the optimal selected angle ϕ, it is expected that one of the parameters a ϕ and b ϕ should describe the scattering strength while the other one describes the effects of the phase function. The data analysis is done for Equation (10) as well as Equation (12). Thereby, it should be noted that the angle ϕ can be derived easily by looking at the results from the Bayesian inference for the parameters a and b.

| Data analysis Bayesian approach
Bayesian statistics is based on the Bayesian interpretation of probabilities. In this case, probability expresses the degree of belief in an event.

Bayesian fitting with MCMC
To briefly explain: The MCMC algorithm selects a random neighboring point. Then, its likelihood is calculated and compared with the likelihood of the previous parameter set. If the likelihood of the new parameter set is greater, the new parameter set is accepted. If the likelihood of the new parameter set is smaller, the new parameter set is accepted with a certain probability. All accepted parameter sets are stored. This process is repeated until the intended number of points has been generated.
For the analysis, the Gaussian MCMC is chosen. Hence, the likelihood is calculated as follows based on the study from Kawahara et al. [17]: where l a, b, c ink ð Þ σ is the likelihood, R λ,a, b ð Þis the calculated spectra by Equations (10) or (12) and σ is standard deviation of the Gaussian function. Normally, σ can be associated with the noise from the measurements [17]. In this study, σ is used as a tool for setting the variations allowed by MCMC and with this the percentage of accepted parameter sets.
The first parameter set (initial values) of the fit are calculated by a prior least square fit. A new parameter set will be selected by randomly choosing a neighboring point for each parameter with the help of Gaussian distributed random number. As the area under the curve is normalized to one, the Gaussian function only depends on σ. The exact parameter values are explained in the next section.
If the newly selected parameter has a lower likelihood, it is still accepted when the following condition is met: where R 0::1 ½ is a uniform random number from zero to one, and l a, b, c ink ð Þ new σ and l a, b, c ink ð Þ old σ are the likelihood of the new and the old parameter set, respectively.
After the desired amount of intended number of parameter sets has been generated, the calculation finishes. To minimize distortions due to the fact that the initial parameter set is not the global minimum, the first 1000 accepted values are removed as burn-in.

Bayesian analysis parameter
As mentioned earlier, in this part of the analysis k 1 and k 2 are taken as single values resulting from the geometric mean for all calculated values for k 1 and k 2 . This is done because it speeds up the calculation process considerably. Nevertheless, the most important results can be derived. First, the fit is performed with Eq. 10. From this fit, the angle ϕ is calculated and used to repeat the analysis with Eq. 12. The parameter σ is set to 1 for the fit of Equations (10) and (12). These values are chosen as they have been shown to give stable results. The new random data point is derived by a Gaussian random number with a standard deviation of 1% of each old parameter. In total, 2 500 000 parameter sets are generated for each Bayesian fit.
Effect of a ϕ ð Þ and b ϕ ð Þ on the spectrum The analysis of the effect of a ϕ ð Þ and b ϕ ð Þ on the spectrum is done by simply plotting the functions for different values of one parameter while the other one is constant.
As a ϕ ð Þ and b ϕ ð Þ do not change their behavior for different values, arbitrary parameters are chosen.

| Variation
This section describes the variation of the results from the fits from Equations (10) and (12). Thus, the results are presented as scatter plots. To show all variations in a single plot, the IL concentration is encoded as color and the ink concentration as symbol. Table 1 summarizes the encodings used in this section. Figure 4 shows the scatter plot of the results from the MCMC algorithm for Equation (10). It can be seen that the variation of the results is larger for low IL concentrations. This is due to the fact that in this range the sourcedetector distance is too small to fully satisfy the conditions for the diffusion approximation. Furthermore, it can be seen that between the constants a and b the variation is along one direction. Therefore, the idea of coordinate transformation seems to be reasonable. Furthermore, it can be seen that the orientation and, thus, the angle (ϕ) is constant regardless of the concentration of IL or ink. For the absorption coefficient, the results are quite heterogeneous.
It should be noted that the angle ϕ does not seem to be dependent on g, k 1 , k 2 , the distance sample to probe or the angle between sample and probe. The last two parameters were tested experimentally. Variations of the input parameters k 1 , k 2 did not change the angle ϕ. The influence of g was tested indirectly by multiplying a spectrum by a function proportional to λ À4 . None of these changes had any influence on the angle ϕ. No other probe could be tested because no other probes were available. At this point, ϕ seems to depend only on the base of the exponent.
The scatter plots from the fits from Equation (12) are shown in Figure 5. The angle ϕ is set to ϕ ¼ À0:334 ¼ À19:14 . In general, the results are much better sorted and all three parameters are almost independent. It should be noted that the results for 1% IL (blue) might be contradictory as the conditions for diffusion theory are not fulfilled. Therefore, these results are not considered in this discussion. Furthermore, in Figure 5 it can be seen that b ϕ clearly describes the scattering strength and, thus, the IL concentration. In Figure 5 (right), it can be seen that the parameters b ϕ and c ink are well distinguishable for b and c ink in contrast to Figure 4 (right). Thus, only these two parameters are relevant for measuring the absorption and scattering for a known substance. In contrast to b ϕ , a ϕ does not seem to be affected by scattering. It is slightly shifted by absorption. Finally, a ϕ and c ink are not independent from each other: The variance of the parameter c ink depends partly on a ϕ . It should also be noted that multiplying R λ ð Þ by a function, which tilts the spectrum does not significantly change the fitting results for b ϕ and c ink while a ϕ increases from around a ϕ ¼ 2 À 5 to around a ϕ ¼ 9 À 12. This is tested with the function 1 þ 0:6 Á λ=400 ð Þ À4 . This function is only chosen as it tilts the spectra and by this simulates more Rayleigh scattering. The scatter plot for a ϕ , b ϕ and c ink is shown in the attachment ( Figure S1). Nevertheless, this happens due to the fact that the source detector separation and the absorber influence the perceived spectral dependence. To illustrate the effect of the new parameters, their variation is visualized in Figure 6. It can be clearly F I G U R E 4 Scatter plots of the fit parameters from Equation (10). The large, bold symbols represent the mean position of all accepted parameter combinations F I G U R E 5 Scatter plots of the fit parameters from Equation (12). The large, bold symbols represent the mean position of all accepted parameter combinations seen that a ϕ influences the spectral dependence. Thus, it might also describe the effect of the phase function. However, it is also possible that a ϕ is influenced by the geometry of the probe and the spectral range of the measurement. At the same time, the parameter b ϕ influences the scattering strength. This is consistent with the prediction that the parameter with the small variance is responsible for the effects of the scattering strength. The parameter a ϕ shows the stronger variance compared to the parameter b ϕ . It may even be possible to use the parameter a ϕ to derive information about the scatterer size and thus the phase function by means of DRS. However, as already shown, the noise (variance) for this parameter is comparatively high. This could limit this possibility. Furthermore, it should be noted that the use of Bayesian statistics automatically provides these constraints correctly. Therefore, it is a powerful tool for reliable data analysis for DRS.

| Reconstruction of the absorption and reduced scattering coefficient
The final reconstruction of the absorption coefficient and the reduced scattering coefficient are shown in Figures 7  and 8, respectively. Both figures show the a posteriori probability density function. The exemplary results for an IL concentration of 2% and an ink concentration of 2.4% are shown. It should be noted that the results match well for the reduced scattering, while for the absorption there is a discrepancy of up to a factor of 2. This error occurs due to the model from Zonios et al. [16]. It is not caused by the Bayesian approach as the standard nonlinear least square fitting shows the same effect.
The absorption coefficient shown in Figure 7 is well reconstructed for both cases. It can be seen that the a posteriori probability density function of the absorption is a little bit narrower for the turned coordinate system.
The reduced scattering coefficient shown in Figure 8 shows strong differences using Equations (10) and (12) for the fitting. It can be seen that the standard function for analyzing the reduced scattering coefficient is inherently unstable. This can be compensated by the proposed coordinate transformation. By this, the results become more realistic as the maximal probability is closer to the real reduced scattering coefficient shown as red line.
In summary, the usage of Equation (12) with the coordinate transformation has no effect on the absorption coefficient compared to Equation (10). For the reduced scattering coefficient, its effect is tremendous. The a posteriori probability distribution from the fitting function gets extremely narrow. Moreover, the maximal a posteriori probability matches much better with the real values.
F I G U R E 6 Effect of variation of the parameters in Equation 11 on the reflectance derived with Equation (12). On the left side, parameter b ϕ is set to À0.334 and parameter a ϕ is varied from 0 to 6.3. On the right side, parameter a ϕ is set to 3.5 and parameter b ϕ is varied from 0 to À0.60 Overall, the analysis using Bayesian inference is suitable for analyzing the stability of the fitting functions. Furthermore, this analysis was able to identify a systematic error in the description of the reduced scattering coefficient, which makes the fitting function unstable. Therefore, strong fluctuations in the final fitting parameters are to be expected which are also present in the current literature [2]. In addition, the instability can be corrected by a coordinate transformation of the function describing the reduced scattering coefficient.
Furthermore, the methodology of Bayesian inference in this study might be applied to many similar problems F I G U R E 7 Exemplary probability density of μ a depending on the wavelength for an IL concentration of 2% and an ink concentration of 2.4%. Left with Equation (10) and right with Equation (12). The red line represents the real value F I G U R E 8 Exemplary probability density of μ s 0 depending on the wavelength for an IL concentration of 2% and an ink concentration of 2.4%. Left with Equation (10) and right with Equation (12). The red line represents the real value as it is a replacement of a nonlinear least square fit. Practical limitations might arise when the amount of fitting parameters gets too high. In this case, the required calculation time will be too long. The second result of the different description for the reduced scattering coefficient from Equation (11) should be able to be used for different descriptions of the diffuse reflectance. However, the parameter ϕ might not be the same as in this study.
With the newly introduced coordinate system, one parameter describes the scattering strength and the other the spectral dependence and thus the effects of the phase function. This allows both effects on the final spectrum to be clearly separated despite the limitations of the diffusion approximation. Furthermore, the effects of the phase function correctly show a high proportion of noise. Nevertheless, a prediction of the phase function could be possible with the presented approach.