Open Access
Add to BibSonomy
Abstract
Analysing diffuse reflectance spectra to extract properties of biological tissue requires modelling of light transport within the tissue, considering its absorption, scattering, and geometrical properties. Due to the layered skin structure, skin tissue models are often divided into multiple layers with their associated optical properties. Typically, in the analysis, some model parameters defining these properties are fixed to values reported in the literature to speed up the fitting process and improve its performance. In the absence of consensus, various studies use different approaches in fixing the model parameters. This study aims to assess the effect of fixing various model parameters in the skin spectra fitting process on the accuracy and robustness of a GPU-accelerated two-layer inverse adding-doubling (IAD) algorithm. Specifically, the performance of the IAD method is determined for noiseless simulated skin spectra, simulated spectra with different levels of noise applied, and in-vivo measured reflectance spectra from hyperspectral images of human hands recorded before, during, and after the arterial occlusion. Our results suggest that fixing multiple parameters to a priori known values generally improves the robustness and accuracy of the IAD algorithm for simulated spectra. However, for in-vivo measured spectra, these values are unknown in advance and fixing optical parameters to incorrect values significantly deteriorates the overall performance. Therefore, we propose a method to improve the fitting performance by pre-estimating model parameters. Our findings could be considered in all future research involving the analysis of diffuse reflectance spectra to extract optical properties of skin tissue.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Diagnostic optical imaging techniques detect light propagating from the imaged sample towards the light detector. Reflected, transmitted or fluorescent light exiting the sample such as biological tissue contains important information on the optical properties of tissue within the interrogated volume. Among others, information about tissue absorption (e.g., blood and melanin content) and scattering (e.g., distribution of scatterers such as collagen fibres) properties can be determined. However, a relation between the sample properties and the light propagation in biological tissue must be established to extract quantitative optical properties from the detected light.
Light propagation in biological tissue is generally modelled using different iterative approaches such as diffusion approximation (DA), adding-doubling (AD) and Monte Carlo (MC) [1]. The latter is widely used for its flexibility and high accuracy but is computationally intensive and therefore time-consuming. DA is simple, fast and can be solved analytically but is highly inaccurate and only applies to strongly scattering media. On the other hand, AD can achieve an arbitrary level of accuracy much faster than MC and is adaptable enough to model various optical properties such as scattering anisotropy and phase function. In contrast to MC and DA, it can only model layered turbid media in a slab geometry [2]. Therefore, it is applicable only to detection schemes where spatially resolved information is not obtained (e.g., diffuse reflectance spectroscopy by integrating spheres, hyperspectral imaging using full-field illumination). In all these approaches, tissue optical properties are implicitly encoded in measured tissue properties such as reflectance [1,3].
In skin tissue, the main absorbers are blood (haemoglobin) and melanin in the visible spectrum, proteins and amino acids in the UV spectrum, and water and lipids in the IR spectrum [4]. The most prominent scatterers include collagen fibres and cellular organelles such as mitochondria and lysosomes [1]. Human skin is a layered tissue consisting of three layers: the epidermis, dermis, and subcutis, each with unique characteristics. To accurately extract the optical properties of skin, its layered structure must be considered in the tissue model, and the main absorbers must be included. However, an open question remains on how accurately the skin should be modelled, i.e. the number of layers that should be simulated and which (and how many) parameters representing optical properties should be accounted for in each layer.
Traditionally, many authors employed a human skin model based on a single homogeneous layer [5,6], two layers [3,7–13] (e.g., a thin upper layer and a semi-infinite bottom layer), or three layers [14,15], better representing the skin structure. Verdel et al. [16] applied a four-layer skin model, where the dermis was subdivided into two structurally different layers, the papillary and reticular dermis. Furthermore, research groups by Meglinski divided human skin into six [17], seven [18–21] or eight [22] layers, resulting in an advanced but exceptionally complex model. Intrigued by the problem, Karlsson et al. [23] showed that a homogenous one-layer optical skin model is incapable of explaining spatially resolved diffuse reflectance spectra at multiple source-detector separations in the visible wavelengths range. Similarly, Hennessy et al. [3] proved that using a one-layer model leads to significant errors in extracted optical properties. Fredriksson et al. [14] demonstrated that a three-layer multi-parameter tissue model is needed to explain spectra at two different source-detector separations for realistic simulated tissue models. However, Bjorgan and Randeberg [24] proved that by exploiting the scale-invariance of the reflectance modelling, it is possible to consider just the upper layer (epidermis) of the skin without completely modelling the lower layers. Importantly, their Monte Carlo simulations indicate that the fitted layer transmittance and reflectance spectra are unique, but an infinite number of physiological skin optical parameters yield almost identical spectra.
Besides partitioning skin tissue models into multiple anatomically justified layers, modelling skin involves considering various absorption and scattering properties. The latter are usually assumed to be the same in all layers [3,9,14,24,25]. On the other hand, absorption properties considered in different studies vary significantly. Generally, the model must include all essential tissue absorbers (also called chromophores) or else the included absorbers will compensate for the missing ones [14]. For example, Lindbergh et al. [6] showed that diffuse reflectance spectra of a human myocardium surface are modelled better with the inclusion of absorbers such as cytochrome aa3 and methemoglobin. Nevertheless, there is no general approach, and the authors seem to adapt the models to their needs in specific cases.
The same is true for another critical question that remains open in the field: which (and how many) tissue parameters integrated within the skin model could be assumed from the literature. Commonly, the authors refer to thorough reviews of reported tissue optical properties by Bashkatov et al. [1] or Jacques [4] to fix the values of selected parameters that the model does not estimate. While doing so, which is necessary to speed up and increase the robustness of the spectra fitting process, it may lead to inaccurate estimation of the fitted optical parameters.
To our knowledge, no previous study investigated the effect of fixing various tissue parameters in the skin model on the accuracy of estimating the remaining free parameters and the robustness of the fitting process. Therefore, one of the main goals of this study was to find the optimal sets of fixed parameters, i.e. the ones that yield the most accurate and precise results with the smallest number of fixed parameters, to speed up the analysis of diffuse reflectance spectra with a GPU-accelerated two-layer inverse AD algorithm (IAD). In order to do this, 15 sets of 11 tissue parameters were selected, where each set of parameters represented a noiseless skin spectrum simulated using the two-layer AD method. Then, these 15 spectra were fitted using the Levenberg-Marquardt (LM) algorithm for 27 different fitting cases, where various combinations of optical parameters were fixed to find the optimal set of fixed parameters. The two-layer skin geometry was used because multiple previous studies demonstrated that two layers are sufficient to simulate skin spectra [3,7–13]. Moreover, Bjorgan and Randeberg [24] proved that reflectance from one-layer dermis models identical to the reflectance of the multi-layer dermis models could be constructed. Thus, skin reflectance spectra can be investigated without fully modelling the deeper layers.
Another chief aim was to show that by assuming the optimal set of fixed parameters, the performance of our IAD algorithm improves significantly. Specifically, the optical parameters estimated by the model are less scattered around the mean values, resulting in much more stable and robust fits and agree better with the actual values, yielding higher accuracy. This hypothesis was tested by adding different noise levels to the simulated skin spectra and studying the effect of noise on the algorithm's performance.
Ultimately, we demonstrated that the two-layer IAD method is an accurate and robust technique for analysing in-vivo measured diffuse reflectance spectra by determining optical parameters from hyperspectral images of human hands.
2. Materials and methods
2.1 Skin spectra simulation
Skin spectra were simulated using the AD algorithm introduced to biomedical optics by Prahl et al. [2]. Inverse AD (IAD) is an iterative technique for numerical solving the radiative transport equation (RTE) in homogeneous turbid media of slab geometry. The method provides vast flexibility in modelling turbid media with any optical properties, such as albedo, refractive index, and optical thickness, considering anisotropic scattering within each medium and light reflection and transmission at the interface of different media. Compared to MCML [26], which is widely used for light transport modelling in multi-layered tissues in the biomedical optics community, the IAD method can model the same geometries if laterally homogenous illumination is considered and light is normally incident upon the sample – which is the case of hyperspectral imaging presented in this study. Generally, IAD allows for any desired accuracy to be achieved, given sufficient time and computing resources [1,2]. Specifically, Prahl et al. [2] verified IAD against inverse MC (IMC) for multiple combinations of absorption and scattering coefficients and anisotropy factors. They found that estimated parameters mainly were within 2-3% of those obtained by IMC, as confirmed in a recent study by Vincely et al. [27].
The IAD algorithm was adopted in this study to the graphics processing unit (GPU) to accelerate the calculation significantly. A two-layer skin model shown in Fig. 1 was considered comprising epidermis and dermis and defined by 11 tissue parameters. The forward AD algorithm simulated light propagation in the skin model and ultimately calculated the reflectance and transmittance spectra. The developed GPU IAD algorithm was extensively tested against MCML for a broad range of optical properties characteristic for biological tissues (absorption coefficient 0.1–1 cm-1, scattering coefficient 25–150 cm-1, anisotropy factor 0.55–0.96, refraction index 1.3–1.5). For 20 fluxes, the agreement between IMC and IAD was absolute down to the fourth decimal place. More information about the MCML and IAD can be found in [28].
Fig. 1. A two-layer skin model with a total of 11 tissue parameters was used in this study. The parameters are as follows: ${f_\textrm{m}}$ – melanin volume fraction, ${f_{\textrm{Hb}}}$ – deoxyhaemoglobin volume fraction, ${f_{\textrm{Hb}{\textrm{O}_\textrm{2}}}}$ – oxyhaemoglobin volume fraction, ${f_{\textrm{brub}}}$ – bilirubin millimolar concentration, ${f_{\textrm{CO}}}$ – reduced cytochrome C oxidase millimolar concentration, ${f_{\textrm{CO}{\textrm{O}_\textrm{2}}}}$ – oxidised cytochrome C oxidase millimolar concentration, $as$ – scattering coefficient, $bs$ – scattering power, ${f_{\textrm{Ray}}}$ – fraction of Rayleigh scattering, ${d_\textrm{e}}$ – epidermis thickness, ${d_\textrm{d}}$ – dermis thickness.
LM algorithm was adopted to GPU to perform non-linear least-squares fitting of the AD simulated reflectance spectra to the measured reflectance spectra. LM was used for non-linear least-squares fitting as it converges faster than other optimisation methods such as the Gauss-Newton method and can find an optimal solution even if the initially guessed parameters are far from the solution [29].
The absorption coefficient of the epidermis and dermis was calculated as [4]:
Parameter ${f_\textrm{m}}$ represents the volume fraction of melanin, ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_\textrm{2}}}}$ are volume fractions of deoxy- and oxyhaemoglobin, ${\mu _{\textrm{a,}\,\textrm{Hb}}}$ and ${\mu _{\textrm{a,}\,\textrm{Hb}{\textrm{O}_\textrm{2}}}}$ are associated absorption coefficients, ${f_{\textrm{brub}}}$ and ${\mu _{\textrm{a,}\,\textrm{brub}}}$ are millimolar concentration and absorption coefficient of bilirubin, whereas ${f_{\textrm{CO}}}$ and ${f_{\textrm{CO}{\textrm{O}_\textrm{2}}}}$ are respective millimolar concentrations of reduced and oxidised cytochrome C oxidase and ${\mu _{\textrm{a,}\,\textrm{CO}}}$ and ${\mu _{\textrm{a,}\,\textrm{CO}{\textrm{O}_\textrm{2}}}}$ the corresponding absorption coefficients. ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_\textrm{2}}}}$ are calculated from the blood volume fraction, B, and blood oxygen saturation, S [4]:
Furthermore, the reduced scattering coefficient was adapted from [4]:
Finally, the refractive index of the epidermis and dermis was calculated as [1]:
2.2 Fitting cases
A total of 15 skin spectra were simulated in this study from the predetermined sets of tissue parameters presented in Table 1. The baseline spectrum (Spectrum 2) represents skin with a 2 % melanin volume fraction, 0.65 % blood volume fraction (B) and 40 % oxygen saturation (S), low bilirubin concentration, absence of oxidised and reduced cytochrome C oxidase, the values of scattering coefficient and scattering power similar to the mean values reported for skin by Jacques [4], no Rayleigh scattering, and standard thicknesses of epidermis and dermis reported by Bashkatov et al. [1]. The tissue parameters for other skin spectra were varied within the specified physiological ranges: melanin volume fraction was varied between 1 % and 5 % [4], blood volume fraction between 0.2 % and 3.2 % [4], oxygen saturation on the interval of 15-90 % [4], bilirubin volume fraction between 1e-7 and 1e-2, as between 30 and 70 [1], bs between 0.7 and 2.1 [1], and epidermis thickness between 50 and 200 µm [1]. The red and green bolded text in Table 1 highlights the optical parameters with values lower and higher than the baseline spectrum, respectively. The values of deoxy- and oxyhaemoglobin, ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_2}}}$, are calculated from the blood oxygenation, B, and oxygen saturation, S, using Eq. (5).
Table 1. Tissue parameters values for each simulated skin spectrum. Legend: * Baseline spectrum, red - value lower than baseline, green - value higher than baseline
Table 2. Free (${\times} $) and fixed (▪) parameters for each fitting case.
Collectively, 27 different cases of fixed parameters selections were studied in this work. Table 2 presents the free (fitted) and fixed parameters for each fitting case, where the free parameters are marked with the cross sign (${\times} $) and the fixed parameters with the filled square sign (▪). In the first case, denoted with the letter a, all parameters were free, whereas in other cases, up to 7 parameters were set to a fixed value. Tissue parameters such as melanin volume fraction, ${f_\textrm{m}}$, deoxy- and oxyhaemoglobin volume fractions, ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_2}}}$, and bilirubin concentration, ${f_{\textrm{brub}}}$, were not fixed in any case because they significantly affect the reflectance spectra features, and inadequate fits are obtained.
2.3 Finding the optimal set of fixed parameters
We investigated the impact of different sets of fixed parameters by repeatedly fitting the 15 selected skin spectra for all 27 fitting cases ten thousand times. This allowed us to find the right balance between the time needed to complete the task and the desired statistical accuracy of the results.
The AD algorithm used to simulate the skin spectra was implemented for GPUs in MATLAB R2020b (Mathworks, MA). Diamond initialisation (DI) was utilised to simulate an initial thin skin layer, and incoming and outcoming light was divided into 20 fluxes to ensure reasonable accuracy. The maximum number of iterations of the LM algorithm was limited to 200, but the mean number of iterations was 39 ± 17. Fitting was performed on a computer consisting of an AMD Ryzen 7 1700X CPU, 16 GB RAM, and an Nvidia Titan Xp graphics card with 12 GB RAM, where the batch size (simultaneous number of GPU threads) was set to 1000. Total fitting time per spectrum was 0.14 to 1.8 s, and the mean value was 0.36 s ± 0.62 s.
Firstly, the coefficient of variation, ${R^2}$, was calculated to evaluate the goodness-of-fit for each fitted spectrum. The deviations in ${R^2}$ originate from different initial parameters, which leads to skin spectra converging to local minima instead of global. Therefore, the spectra for which ${R^2} < 0.99$ were discarded for diverging too much from the simulated ones, while the remainder was used in the analysis. Secondly, the determined parameters were normalised by dividing them with the true parameters corresponding to the predetermined skin spectra. Then, root mean square error (RMSE) for each normalised parameter of the remaining spectra was calculated, resulting in one numerical value per parameter per simulated spectrum per case. Additionally, the RMSE of fitted spectra was computed, obtaining one numerical value per simulated spectrum per case.
Furthermore, squares of RMSE values of all predefined spectra within each case were summed for each parameter. These cumulative values were characterised as measures of the accuracy of the fitting process for the particular optical parameter within individual fitting cases. The summation was repeated for all RMSE of the fitted spectra. As a result, the accuracy of the fitting process for each fitting case was described by 12 numerical values: 11 cumulative $\textrm{RMSE}_{\textrm{param},\textrm{i}}^2$ values for the individual fitted parameters ($\textrm{i}\, = \;\,1,\;2\ldots \,11$), and one cumulative value for the fitted spectra, $\textrm{RMSE}_{\textrm{spec}}^2$. Importantly, for the fixed parameters, the calculated RMSE value was 0.
Finally, a function, ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$, was defined to estimate the performance of the fitting process for different fitting cases:
The protocol for calculating ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ is visually presented in Fig. S1 (see Supplement 1). A similar approach was adopted to derive another estimation function based on the calculation of mean absolute error (MAE):
The first estimation function, ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$, is based on the calculation of RMSE because it is a standard measure of the differences between the observed values and values predicted by a model. Unlike RMSE, a quadratic scoring rule, MAE is a linear score that measures the average magnitude of the errors and is known to be less sensitive to outliers than RMSE. Thus, another estimator function based on MAE was utilized to improve the analysis.
The accuracy and robustness of the IAD algorithm for simulated spectra were evaluated as the agreement of mean (or median) values of fitted parameters and actual values and as a dispersion of fitted parameters, respectively. The performance was also evaluated using the two estimation functions, ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$.
2.4 Fitting skin spectra with noise
In addition to fitting noiseless skin spectra, predetermined simulated skin spectra were augmented by adding different noise levels and fitted using the IAD method to show that the proposed approach works with realistic, noisy spectra. Specifically, three different noise levels with a constant signal-to-noise ratio (SNR) were added, namely 30, 50 and 100 dB. The middle noise level was estimated to correspond average SNR of the experimental hyperspectral imaging system, whereas the first and the last SNR values represent the cases with a high and low degree of experimental noise, respectively. We avoided adding more noise to the simulated spectra (decreasing SNR) to comply with the experimental values obtained on our hyperspectral system presented in 2.5.
Noisy spectra, ${S_n}$, were generated using the following equation:
where ${S_0}$ is the noiseless spectrum, n is normally distributed random noise, and $\alpha $ is a scalar coefficient that yielded a predefined SNR.2.5 Fitting in-vivo skin spectra
The last step involved fitting measured in-vivo skin spectra recorded with our hyperspectral imaging (HSI) system thoroughly described in [30]. Briefly, the system consists of two modules, namely the push-broom hyperspectral imaging module and the optical profilometry (OP) module. The former is used to record hyperspectral images in the 400-1000 nm wavelength range with 2 nm spectral resolution, whereas the latter provides the 3D surface shape of the imaged sample needed for the curvature correction of the recorded images [30].
Raw hyperspectral image, ${I_\textrm{raw}}$, was first converted to reflectance, ${I_\textrm{ref}}$, using equation [31]:
Six healthy volunteers aged 23 to 24 were included in the study. Three hyperspectral images of fingers were recorded for each participant: before, during, and after arterial occlusion [16]. During the arterial occlusion, the blood pressure cuff placed around the upper arm was inflated to 200 mmHg and left for approximately 3 minutes. However, one image obtained during arterial occlusion was discarded due to insufficient OP data, resulting in 17 hyperspectral images analysed with the IAD algorithm. The experimental protocol conforms to the principles expressed in the Declaration of Helsinki and was approved by the Slovenian National Medical Ethics Committee (0120-629/2016-3; KME 66/01/17). Informed consent was obtained from all subjects included in the study.
Since the ground truth values of the experimental skin parameters are unknown, two approaches were tested. In the first approach, the values of all model parameters were adopted from the literature (Table 3). In another approach, optimal initial parameter values were pre-estimated before attempting the actual fitting process of the whole image. Specifically, for each processed hyperspectral image of a hand, a 20 × 20 px region of interest (ROI) was selected from which a mean spectrum was calculated. The mean spectrum was then repeatedly fitted 1000 times using different randomly selected initial parameter values and the parameters from the single best fit (the one with the highest ${R^2}$ value) were taken as the optimal parameters, and the parameters were either fixed or further estimated (free parameters) in the actual image fitting process accordingly.
For experimental spectra, the robustness was assessed similarly to the simulated data, whereas accuracy was estimated by comparing fitted values to other studies since the actual values of parameters are unknown.
Table 3. The values of initial model parameters used in the first approach of fitting in-vivo measured hyperspectral images of human hands.
3. Results
3.1 Spectra simulation
A total of 15 skin spectra simulated from tissue properties presented in Table 1 with the two-layer AD algorithm used in this study are shown in Fig. 2 Reference source not found.. The baseline spectrum (Spectrum 2) is coloured in a thick light blue line. Notably, different values of tissue parameters affect spectra shape significantly. For example, a 5 % melanin content in Spectrum 3 (purple line) markedly reduces the reflectance on the entire wavelength interval and obscures other spectral characteristics.
Fig. 2. Skin spectra simulated with a two-layer AD algorithm. Numbers 1 to 15 represent the serial number of the spectra, and 2 is the baseline spectrum (thick light blue solid line).
3.2 Fitting noiseless simulated spectra
3.2.1 No fixed parameters (fitting case a)
Four examples of noiseless simulated skin spectra fitted with all 11 parameters free (fitting case a in Table 1) are shown in Fig. 3. Figure 3(a-d) correspond to Spectrum 1, 5, 9, and 14, respectively. The simulated skin spectra (blue line) are presented alongside the mean fitted spectra (orange line) and the standard deviation of the fitted spectra (orange shade). Notably, the mean fitted spectra agree well with the simulated spectra and are generally within the ±1 STD range of the simulated spectra on the entire wavelength interval. However, slight deviations are present for wavelengths up to 450 nm. Figure 3(e-h) show the associated values of the fitted parameters plotted conveniently in box plots. On each blue box, the central mark indicates the median value, and the bottom and top edges indicate the 25th and 75th percentiles, respectively. One can observe that despite fitted spectra corresponding to the measured spectra (Fig. 3(a-d)), the values of estimated parameters in Fig. 3(e-h) are somewhat scattered, as also indicated by many outliers (red plus marks). Generally, the median values of estimated parameters such as ${f_\textrm{m}}$ are higher than the actual values (green diamond marks). However, the actual values are generally within the 25th to 75th percentile of the fitted values except for some parameters such as $bs$ in Fig. 3(g), meaning that the fitted values of optical parameters agree well with the actual values. To support this statement, numerical data for these four cases are presented in Table 4. Evidently, the mean values of the fitted parameters are within ±1 STD from the actual ones, but STD for most fitted parameters is of the same order of magnitude as the mean fitted parameters values.
Fig. 3. Noiseless simulated skin spectra fitted without any fixed parameters (fitting case a in Table 2). Solid blue lines represent noiseless simulated skin spectra, solid orange lines represent the mean fitted spectra, and orange shaded areas represent one standard deviation region. The examples presented are a) Spectrum 1, b) 5, c) 9, and d) 14. Figure 3(e-h) represent the corresponding values of the fitted parameters in boxplot diagrams. On each blue box, the central mark indicates the median, and the bottom and top edges indicate the 25th and 75th percentiles, respectively. The whiskers extend to the most extreme data points not considered outliers, and the outliers are plotted individually using the red plus (+) marker symbol. The actual value for each parameter is plotted with a green diamond.
Table 4. Mean values of the fitted parameters and the corresponding standard deviations for four different spectra (Spectrum 1, 5, 9, and 14) in the case of no fixed parameters (fitting case a). All fitted parameters are within the range of the actual values.
The correlation coefficients between the parameters were calculated and are graphically presented in Fig. S2 (see Supplement 1) for the selected four spectra. Multiple pairs of the parameters express a high degree of either positive (e.g., ${f_{\textrm{Hb}}}$ and ${f_{\textrm{brub}}}$ in the case of Spectrum 9 in Fig. S2(c) or negative correlation (e.g., $as$ and ${d_\textrm{e}}$ in the case of Spectrum 1 in Fig. S2(a), whereas no correlation is found for some parameters (e.g., ${f_{\textrm{CO}}}$ and ${f_{\textrm{Ray}}}$ in the cases of Spectrum 9 and 14). Overall, a high degree of negative correlation exists between ${f_\textrm{m}}$ and ${d_\textrm{e}}$ due to the same effect on the reflectance spectra (both a high ${f_\textrm{m}}$ and a thick epidermis reduces the reflectance). Therefore, a high ${f_\textrm{m}}$ results in a thin epidermis, and vice-versa. However, the degree of correlation between different pairs of parameters varies significantly for simulated spectra, indicating that individual values of parameters vastly affect the fitting process. As a result, no typical pattern for the degree of correlation between parameters was found for different simulated spectra.
3.2.2 Fixed $bs$, ${f_\textrm{Ray}}$, ${d_\textrm{e}}$, ${d_\textrm{d}}$ (fitting case q)
In the following case, two scattering parameters and the layer thicknesses were fixed to demonstrate the effect of parameters fixation on the fitting performance. Figure 4(a-d) show the simulated spectra (Spectrum 1, 5, 9 and 14), the mean fitted spectra and the fitted spectra standard deviation. Similarly to the previous subsection, the mean spectrum generally agrees well with the actual simulated spectrum on the entire 430-750 nm interval, except for minor deviations below 450 nm. Figure 4(e-h) compare the actual parameters values of the simulated spectrum (green diamond) and the values obtained by fitting (box plots). Evidently, some of the fitted values differ from actual values for more than ±1 STD, as presented in Table 5 (red text). Specifically, $as$ is inaccurate for three of the four presented spectra, possibly as a results of other scattering parameters being fixed. Nevertheless, this case still yields more precise results than case a, as demonstrated by generally lower STD values of estimated parameters in Table 5 than Table 4. In Fig. 4(e-h), the fixed parameter’s values are not presented since they match the actual values.
Table 5. Values of the mean fitted optical parameters and the corresponding standard deviations for four different spectra (Spectrum 1, 5, 9, and 14) for the fitting case q. The mean values of fitted parameters that differ for more than ±1 STD from the actual values are highlighted in red.
Fig. 4. Noiseless simulated skin spectra fitted with four parameters fixed: (fitting case q in Table 2). Solid blue lines represent noiseless simulated skin spectra, solid orange lines represent the mean fitted spectra, and orange shaded areas represent one standard deviation region. The examples presented are a) Spectrum 1, b) 5, c) 9, and d) 14. Figure 4(e-h) represent the corresponding values of fitted parameters in boxplot diagrams. The actual value for each parameter is plotted with a green diamond.
The correlation matrices for the selected spectra and all pairs of fitted parameters are shown in Fig. S3. Visual comparison to the matrix in Fig. S2 shows that different parameters correlate depending on the fixation. For example, parameters ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ express a high degree of negative correlation in the fitting case q, but are much less correlated in the fitting case a where all parameters are free – likely due to these two parameters compensating the fixation of others in case q, effectively acting contrary to one another.
3.2.3 Finding the optimal set of fixed parameters
Following the procedure outlined in 0, we evaluated the effect of fixing different sets of optical parameters on the accuracy and robustness of the fitting process, employing two different estimators, Eq. (9) and Eq. (11), and applying different weights ${w_1}$ and ${w_2}$. The value of weight ${w_1}$ was incrementaly changed from 0 to 1, with a 0.1 step, whereas the value of weight ${w_2}$ was calculated from Eq. (10). Figure 5(a-b) show the values of the two estimation functions for different fitting cases, while Fig. 5(c) shows the absolute difference between the two. It can be seen from Fig. 5(a) that an increase in ${w_1}$ (decrease in ${w_2}$) leads to an increase in the value of ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ – except in case h, where the contribution of both RMSE is similar. It implies that the RMSE of the parameters is generally higher than the RMSE of spectra. The RMSE of the parameters decreases substantially (Fig. 5(a), w1 = 1) as ${f_\textrm{CO}}$ and ${f_{\textrm{CO}{\textrm{O}_2}}}$ are fixed starting from case s. As the estimated values of these two parameters have high STD values (see Table 4 for case a and Table 5 for case q), fixing them reduced the RMSE of parameters considerably. A similar observation holds for ${f_\textrm{Ray}}$. On the other hand, the RMSE of spectra (Fig. 5(a), w1 = 0) does not show such a trend and is dependent on the specific selection of the fixed parameters. The graphs in Fig. 5(b) do not follow a similar trend but instead fluctuate, indicating that the MAE of parameters and spectra have comparable values and are thus both sensitive to reduction of the free parameters. Moreover, MAE values are less sensitive to specific parameters than RMSE values.
Fig. 5. The values of two estimation functions: a) ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and b) ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$. c) The absolute difference of ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$. Weights ${w_1}$ and ${w_2}$ occupy discrete values between 0 and 1, with a step of 0.1.
The minima of both fRMSE and fMAE are present for the same fitting cases (w and y – zz), indicating that both measures favour a lower number of free parameters. Finally, the absolute difference of the two functions seen in Fig. 5(c) is generally the lowest for ${w_1} = 1$ (${w_2} = 0$) and the highest for ${w_1} = 0$ (${w_2} = 1$), suggesting that the RMSE and MAE values of the parameters are similar, but the RMSE and MAE values of spectra are not. Based on the above observations, MAE is better for assessing fitting performance.
The absolute minimum values of ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ are obtained for cases z and zz (Fig. 5(a)). As the value of ${w_1}$ increases (${w_2}$ decreases), the values of ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ for case z become lower than those of case zz – the shift occurring at around ${w_1} = 0.2$ (${w_2} = 0.8$). A similar pattern can be seen in Fig. 5(b) for the values of ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$ for cases z and zz, only this time the shift occurs at approximately ${w_1} = {w_2} = 0.5$.
Figure 6 shows the leaderboards for the values of two estimation functions, namely a) ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and b) ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$. The purpose of these images is to aid the process of finding the optimal set of fixed parameters (the most optimal fitting case). Images visually present the cumulative number of times each fitting case ranks 1 (the highest) to 27 (the lowest) in terms of the value of the estimation function as the weights are gradually changed. For example, as seen in the top right corner of Fig. 6(a), fitting cases z and zz are the only ones ranking the highest on the leaderboard, meaning they had the lowest value of the ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and thus performed the best in terms of accuracy and robustness. Specifically, case z ranked the highest in 9 instances (colour-coded in yellow) and case zz in 2 instances (bluish-purple colour). A similar pattern is apparent in Fig. 6(b), where case z ranked the highest five times (in teal), and case zz ranked the highest a total of 6 times (in light green). The worst performing fitting cases are h in terms of ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$, which ranked the lowest in 10 instances (in yellow), and k in terms of ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$, which ranked the lowest 8 times (in orange). Notice that for both estimators, leaving all parameters free as in case a produces better results than fixing inappropriate sets of parameters, such as cases b, f, h and l.
Fig. 6. The leaderboards for estimation functions a) ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and b) ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$ demonstrate the success of the fitting process for different fitting cases and all combinations of weights ${w_1}$ and ${w_2}$. Colour-coded is the number of occurrences, i. e. the number of times each fitting case takes a particular place in the leaderboard when the values of estimation functions are ranked from 1 – highest to 27 – lowest.
The normalised cumulative RMSE and MAE values of estimated optical parameters (first terms in Eq. (9) and Eq. (11), respectively) were compared to normalised cumulative RMSE and MAE values of fitted spectra (second terms in Eq. (9) and Eq. (11), respectively) with respect to the number of fixed optical parameters as seen in Fig. 7. Specifically, normalised cumulative RMSE and MAE values of model parameters were plotted with orange circles, and normalised cumulative RMSE and MAE values of fitted spectra were plotted with blue circles, where each circle represents one of the fitting cases (a-zz). Filled circles represent fitting cases where both ${f_\textrm{CO}}$ and ${f_{\textrm{COO}_{2}}} $ or ${f_\textrm{Ray}}$ are fixed. Solid blue and orange lines with error bars show the corresponding mean spectra and model parameters RMSE and MAE values.
The purpose of calculating the normalised cumulative RMSE and MAE values is to assess the contribution of parameters error and fitted spectra error in Eq. (9) and Eq. (11) separately and find out how the contributions change as the number of fixed parameters increases. Importantly, graphs in Fig. 7(a) show that mean RMSE values for parameters remain constant as the number of fixed parameters increases from 0 to 4, followed by a significant decrease as more parameters are fixed. Meanwhile, the RMSE of spectra peaks at two fixed parameters, then gradually decreases. A similar trend is observed for mean MAE values in Fig. 7(b), but the mean MAE values of spectra culminate at four fixed parameters before gradually dropping. It demonstrates that as the number of fixed parameters increases to a particular value, the fitted spectra agree less with the actual spectra, whereas the variance of parameters remains relatively unchanged. Then, as the number of fixed parameters is further incremented, the fitted spectra are in a much greater agreement with actual ones, and the variance of optical parameters decreases, resulting in a more accurate and robust fitting process.
Fig. 7. a) Normalised cumulative a) RMSE and b) MAE values of fitted spectra (blue circles) and estimated parameters (orange circles) with respect to the number of fixed model optical parameters. Solid orange and blue lines with error bars represent the corresponding mean values. Filled circles represent cases where both ${f_\textrm{CO}}$ and ${f_{\textrm{COO}_{2}}} $ or ${f_\textrm{Ray}}$ are fixed.
Filled circles representing fitting cases (where both ${f_\textrm{CO}}$ and ${f_{\textrm{COO}_{2}}} $ or ${f_\textrm{Ray}}$ are fixed) confirm our previous findings that these parameters systematically improve the fitting performance. Specifically, fitting these parameters resulted in the lowest normalised cumulative RMSE of parameters and spectra for two to four fixed parameters.
3.3. Fitting noisy simulated spectra
Figure 8(a) shows Spectrum 1 with three different noise levels added, i.e. noise with constant SNR values of 30, 50 and 100 dB, respectively. It can be seen that spectra with the average experimental SNR of 50 dB and a low degree of noise (SNR of 100 dB) agree very well, coinciding with the noiseless skin spectrum (black dashed line with crosses) almost perfectly. On the other hand, a low SNR value results in a much higher noise, altering the simulated skin spectrum but preserving its original shape.
Two representative fitting cases, a (all free parameters) and z (five free parameters), were further studied by applying the three different noise levels. The latter is the overall top-performing case from the previous section with six fixed parameters, whereas the former is a relatively well-performing case with all free parameters and the most degrees of freedom. Figure 8(b-c) show the simulated spectrum (Spectrum 1) with added 30 dB noise, the mean fitted spectrum, and the standard deviation of the fitted spectra for fitting cases a and z, respectively. Notably, the standard deviation in case z is much lower (on average for around 50%) than in case a, suggesting better matching of fitted and noisy simulated spectra. Figure 8(d-e) compares actual and estimated values of all model parameters. The relative deviation of mean parameters from the actual values for the first three parameters in case a is 98%, 50%, and 165%; it is equal to 44%, 7% and 3% in case z, respectively. Parameter estimations are therefore more accurate in case z than in case a. Moreover, the parameters are less scattered in case z, as demonstrated by a lower standard deviation of parameters than in case a. Specifically, STD of ${f_\textrm{m}}$, ${f_{\textrm{Hb}}}$ and ${f_{\textrm{Hb}{\textrm{O}_\textrm{2}}}}$ was 5, 8 and 3 times higher in case a than in case z.
Fig. 8. a) An example of a simulated skin spectrum (Spectrum 1) with different noise levels added. Specifically, noise levels with a constant SNR value were added, namely 30, 50 and 100 dB. b, c) Simulated skin spectrum (Spectrum 1) with added 30 dB noise, the mean fitted spectrum and the standard deviation of fitted spectra in cases a and z, respectively. d, e) The comparison of actual and estimated optical parameters in cases a and z, respectively.
As in the previous section, we evaluated the fitting performance of cases a and z with different noise levels applied. Figure 9 shows the values of estimation functions ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$. Note that the values of the two estimation functions follow a common trend. More precisely, when only the error of spectra is considered (${w_1} = 0$), both estimators have a high value for cases a and z with the highest level of noise added (30 dB), implying the most considerable mismatch between simulated and fitted spectra. As ${w_1}$ is gradually increased, and parameters error is taken into account, the values of the estimation functions are becoming significantly lower for case z than case a and generally decrease as less noise is present in the simulated spectra. These results show that reducing the number of free parameters improves the performance of the fitting process significantly, also in the case of noisy spectra.
Fig. 9. The values of the two estimation functions for fitting cases a and z with different noise levels added: a) ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and b) ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$. Weights ${w_1}$ and ${w_2}$ occupy discrete values between 0 and 1, with a step of 0.1.
This claim is supported by Fig. 10, which shows the leaderboards for both estimation functions. As can be seen, the majority of occurrences lie on the diagonal between the bottom left (all free parameters with the highest level of noise added rank the lowest) and the top right (six fitted parameters with the lowest noise level rank the highest).
Fig. 10. The leaderboards for estimation functions a) ${f_{\textrm{RMSE}}}({{w_1},{w_2}} )$ and b) ${f_{\textrm{MAE}}}({{w_1},{w_2}} )$ for fitting cases a and z with different noise levels added and all combinations of ${w_1}$ and ${w_2}$.
3.4. Fitting in-vivo measured skin spectra
3.4.1 First approach
Shown in Fig. 11 are in-vivo measured (solid blue line) and fitted (solid orange line) skin spectra from selected pixels in the middle finger of a human hand (participant 2) for different phases of the arterial occlusion test. Measured spectra shown in the top row were fitted using case a (all free parameters), whereas the spectra in the bottom row were fitted using case z. In both cases, the first fitting approach was taken where initial parameters from Table 3 were input to the IAD algorithm. Notably, the measured and fitted spectra match very well for all stages of the arterial occlusion test: before (a, d), during (b, e), after (c, f).
Fig. 11. In-vivo measured (solid blue line) and fitted (solid orange line) skin spectra fitted using case a (a-c) and case z (d-f) of participant 2 in different stages of the arterial occlusion test: before (a, d), during (b, e), and after (c, f).
Figure 12 shows colour maps of three optical parameters extracted from a hyperspectral image of a human hand (participant 2) recorded before the arterial occlusion test: ${f_\textrm{m}}$ (a, d), ${f_\textrm{Hb}}$ (b, e) and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ (c, f). Color maps in the top row were extracted using fitting case a and those in the bottom row using case z. The mean values of the estimated parameters are: the mean values of ${f_\textrm{m}}$ are 0.71% ± 0.27% (case a) and 0.87% ± 0.50% (case z), the mean values of ${f_\textrm{Hb}}$ are 1.79% ± 0.69% (case a) and 2.04% ± 0.63% (case z) and the mean values of ${f_{\textrm{Hb}{\textrm{O}_2}}}$ equalled to 1.94% ± 1.04% (case a) and 1.70% ± 1.00% (case z). All mean values were calculated for the inner part of the hand and fingers of participant 2 recorded before the vascular occlusion test, where the surface inclination angle was lower than 45°. In this way, we removed the hand edges, where the spectra are altered significantly due to signal loss, resulting in incorrect parameter estimation.
Fig. 12. Colour maps of ${f_\textrm{m}}$ (a, d), ${f_\textrm{Hb}}$ (b, e) and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ (c, f) extracted from a hyperspectral image of a human hand (participant 2) recorded before arterial occlusion test using fitting case a (top row) and case z (bottom row). In this fitting approach, the initial parameters presented in Table 3 were input in the IAD algorithm. Framed areas represent zoomed-in distributions of parameters on the middle finger outlined by the black rectangle.
All colour maps displayed for case z (Fig. 12(d-f)) appear more continuous and less grainy than for case z (Fig. 12(a-c)), highlighting detailed features such as wrinkles and blood vessels. At first glance, though, the colourmaps of ${f_\textrm{m}}$ might seem more heterogeneous for case z (Fig. 12(d)) than case a (Fig. 12(a)) – which is attributed to the absence of details in case z, also reflected by a higher standard deviation of ${f_\textrm{m}}$ in case z than in case a. These findings can be generalised for all in-vivo measured hyperspectral images used in this study.
To confirm our visual findings, we applied the 2D fast Fourier transform (FFT) to each parameter image of all participants included in this study and calculated the mean ratio of magnitudes for cases a and z with respect to annuli with increasing inner and outer radii. These annuli with increasing inner and outer radii were generated to cover the entire spatial frequency domain space (k-space) of a spatial Fourier transform of 2D parameter distribution maps. Figure 13(a) shows the 11 annuli with a 10 px width used in the calculation, whereas Fig. 13(b) displays the mean FFT magnitude ratio for five optical parameters inside respective annuli. The latter figure shows that the ratios generally increase as we move away from the centre of the k-space. Since the areas away from the centre and towards the edges of the k-space have a higher spatial frequency, this indicates that parameter images extracted using case a contain signal with higher spatial frequency and are thus noisier than those in case z.
Fig. 13. a) 11 annuli with varying inner and outer radii used in the analysis of FFT parameter images. b) Mean FFT magnitude ratio (case a to case z) of estimated optical parameters as a function of annulus number.
3.4.2 Second approach
In another approach described in Section 2.5, initial fitting parameters were pre-estimated rather than instinctively input in the IAD algorithm. A comparison reveals that colour maps in case a (Fig. 14(a-c)) are now less textured compared to their equivalents in Fig. 12(a-c), whereas those in case z (Fig. 12(d-f) and Fig. 14(d-f)) generally maintain the contrast and the texture of parameter images in the first approach.
Table 6. The mean values of ${{\boldsymbol f}_{\mathrm m}}$, ${{\boldsymbol f}_{{\mathrm {Hb}}}}$ and ${{\boldsymbol f}_{{\mathrm {Hb}}{{\mathrm O}_2}}}$ for two different approaches and fitting cases (a, z).
Fig. 14. Colour maps of ${f_\textrm{m}}$ (a, d), ${f_\textrm{Hb}}$ (b, e) and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ (c, f) extracted from a hyperspectral image of a human hand (participant 2) recorded before arterial occlusion test using fitting case a (top row) and case z (bottom row). In the second fitting approach, the initial parameters were pre-estimated and input in the IAD method. Framed areas represent zoomed-in distributions of parameters on the middle finger outlined by the black rectangle.
The mean values of parameters for participant 2 are higher than in the first approach: the mean values of ${f_\textrm{m}}$ are 1.66% ± 0.42% (case a) and 2.11% ± 0.98% (case z), the mean values of ${f_\textrm{Hb}}$ are 2.30% ± 0.87% (case a) and 2.25% ± 0.87% (case z) and the mean values of ${f_{\textrm{Hb}{\textrm{O}_2}}}$ equal to 2.76% ± 1.35% (case a) and 2.59% ± 1.21% (case z). The results for both approaches and fitting cases are summed up in Table 6.
Shown in Fig. 15 is the mean FFT magnitude ratio of five optical parameters computed inside the annuli showcased in Fig. 13(a). Notably, the ratios remain stationary as we move away from the centre of the k-space, demonstrating that in this approach, the quality of parameter images in terms of noise in case a is equal to case z.
Fig. 15. Mean FFT magnitude ratio (case a to case z) of estimated optical parameters as a function of annulus number.
4. Discussion
Several significant implications of this study will be discussed categorically in this section. These findings can guide the development of new skin models and fitting diffuse reflectance spectra recorded with optical imaging techniques such as HSI or diffuse reflectance spectroscopy (DRS).
Estimating skin properties from diffuse reflectance spectra requires modelling skin tissue, considering its layered anatomical structure and absorption and scattering properties originating from absorbers and scatterers present in the tissue. Although many studies employed a large variety of skin tissue models and included various optical properties, little research has been done to assess the model’s performance at fixing different model parameters. Nevertheless, this is a crucial component in the spectra fitting process as adopting reported values of optical properties from the literature can significantly alter the accuracy and precision of estimating variable model parameters. Therefore, this work provides a thorough analysis of this topic.
To begin with, a visual comparison of four simulated noiseless skin spectra (Spectrum 1, 5, 9 and 14) in Fig. 3(a-d) fitted with all free parameters (case a) and the identical spectra fitted with four fixed parameters (case q) in Fig. 4(a-d) reveals only minor differences between the two cases. On the contrary, a comparison of estimated optical parameters values shown in Fig. 3(e-h) and Fig. 4(e-h) exposes a significant difference in the standard deviation between the two fitting cases, supported by numerical data presented in Table 4 and Table 5. An almost identical fitting performance in cases a and q with respect to spectra results in an overall lower standard deviation of estimated optical parameters in case q. However, as Table 5 suggests, many estimated values in case q are far-off from the true ones. From this, two key findings can be generalised to other fitting cases. Firstly, fixing parameters leads to a less accurate but more precise (robust) estimation of variable parameters. Secondly, irrespective of the fixed parameters, an infinite number of spectra can be simulated from the optical parameters that match the actual simulated skin spectrum almost perfectly but yield significantly different values of estimated optical parameters. This result is in accordance with the findings of Bjorgan and Randeberg [24], who observed the same phenomenon for different model geometries used in their IMC model.
Furthermore, fixing parameters reduces the model’s number of degrees of freedom, usually lowering the standard deviation of estimated parameters while resulting in (visually) equivalent skin spectra. However, parameters are correlated due to skin reflectance modelling, and fixing different sets of model parameters changes the correlation between the pairs of parameters, as shown in Fig. S2 and Fig. S3. Therefore, one cannot simply conclude that fixing more parameters will result in more compliant spectra and less scattered values of estimated optical parameters.
A detailed analysis of our model’s performance based on the calculation of normalised cumulative RMSE (Eq. (9)) and MAE (Eq. (11)) values for skin spectra and optical parameters confirms the presence of an interesting pattern: some parameters such as ${f_\textrm{CO}}$, ${f_{\textrm{CO}{\textrm{O}_2}}}$ and ${f_\textrm{Ray}}$ greatly contribute to the error in the parameters but barely affect the spectrum shape and thus its error. On the other hand, a slight change in parameters such as ${f_\textrm{Hb}}$ and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ could affect the spectrum considerably. Thus, it is expected that fixing parameters such as ${f_\textrm{CO}}$ and ${f_{\textrm{COO}_{2}}} $ could systematically improve the fitting performance, and that case s is superior to most other cases where the majority of parameters are fixed (see Fig. 5(a)).
If the fitting cases are evaluated with respect to the number of fixed optical parameters instead of case-by-case comparison, another compelling conclusion can be drawn. Figure 7 demonstrates that as more optical parameters are fixed, up to a certain point, the overall error of the estimated optical parameters remains stationary, whereas the RMSE and MAE errors of spectra increase in value due to poorer spectra matching owing to fewer degrees of freedom of the model. Specifically, the fitting algorithm has fewer possibilities in fitting the spectrum than when all parameters are free, leading to worse spectra agreement. Meanwhile, the parameter space is still large enough (infinite) that the error of the estimated parameters remains unchanged. As more parameters are fixed after the shifting point, errors in both spectra and variable optical parameters gradually decrease – the parameter space is now reduced substantially, leading to much fewer remaining combinations of fitted spectra and a smaller number of local minima. This finding implies that fixing too few parameters may not help to improve the overall fitting accuracy but rather deteriorate it; yet, fixing even more parameters could result in much more robust and accurate simulated skin spectra fitting.
The analysis of two selected fitting cases, a and z, with different noise levels added to the simulated skin spectra, offers some expected results: for a particular case, the model’s general fitting performance (the error) is worst for the highest level of noise added (SNR of 30 dB), and best for the lowest noise applied (SNR of 100 dB) as seen in Fig. 9. What is more interesting is that the performance in case z is superior to performance in case a, as demonstrated by the lower standard deviation of spectra (Fig. 8(b-c)) and estimated optical parameters (Fig. 8(d-e)), confirming our previous findings in case of the noiseless simulated spectra.
In addition to the simulated spectra, experimental skin hyperspectral images were analysed to showcase the IAD algorithm’s performance on real-world data. In the biomedical optics field, the IAD method was primarily used in DRS measurements involving single- or double-integrating sphere systems to determine the optical properties of tissue-mimicking phantoms [32–36], biological tissues [1,37–41] and tissue components such as blood [42] and melanin [43]. However, similarly to widely known MCML [26] used in the light transport modelling of multi-layered tissue, IAD can be applied to other detection schemes such as hyperspectral imaging, provided that the spatially resolved information is not obtained, that laterally uniform illumination is considered, and light is normally incident upon the sample.
Figure 11 clearly shows that the algorithm performs exceptionally well in fitting cases a and z, with an almost perfect match between the measured and fitted skin spectra. A comparison of colour maps for different parameters of participant 2 (recorded before vascular occlusion) in Fig. 12 reveals a better performance of fitting case z than case a. Specifically, parameter distribution maps extracted with case z appear less textured, grainy and more continuous than in case a. As a result, the visibility of details such as wrinkles and blood vessels in case a is compromised, and image contrast deteriorates. What is more, the mean FFT magnitude ratios shown in Fig. 13(b) for different parameters increase with higher spatial frequency, confirming that parameter maps in case a contain much more noise than those in case z. From this, we can conclude that fixing multiple parameters indeed improves the parameter extraction from hyperspectral images.
In the first approach, the initial set of parameters was adopted from the literature (see Table 3). Instead, if the parameters are pre-estimated using the second approach described in 2.5, the smoothness of colour maps can be improved significantly in case a, as seen by comparing colour maps in Fig. 12(a-c) and Fig. 14(a-c). On the contrary, there is no visible improvement in terms of smoothness (or texture) in case z (colour maps in the bottom rows in Fig. 12(d-f) and Fig. 14(d-f)). All maps are equally smooth, although the values of parameters are generally higher in the second approach. We hypothesise that by pre-estimating model parameters, results are more improved in case a because it is less robust than case z, where multiple model parameters are fixed. Pre-estimated parameters are input into the fitting algorithm to perform the actual fitting process. Because these pre-estimated parameters provide a much better initial guess than the general assumptions (Table 3), the estimated skin parameters are less scattered. However, in case z, the fitting algorithm is intrinsically more robust due to multiple fixed parameters, diminishing the effect of parameter pre-estimation on the parameter estimation robustness. Our claim is supported by Fig. 15, which shows that the ratios of FFT magnitude of estimated parameters remain constant as the spatial frequency increases – meaning that there is now an equal amount of noise in parameter images of both fitting cases.
As indicated by the higher mean values of estimated parameters in the second approach (see Table 6), parameter pre-estimation forces the fitting algorithm into a different, possibly global minimum, than the first approach. Undoubtedly, this leads to the algorithm searching parameter values close to the pre-estimated values and making the parameter distributions smoother. Such parameter maps generally correspond to natural skin but might overlook sudden heterogeneities such as skin moles (nevi) and blood vessels. However, the LM fitting algorithm can find the optimal solution even if the initial guess is far-off, and as seen by comparing colourmaps in Fig. 12 to Fig. 14, minor changes in parameters values are well pronounced in the latter. In fact, as the colourmaps in Fig. 14 contain less noise than in Fig. 12, the signal-to-noise ratio (SNR) is increased, and heterogeneities should be more noticeable.
To demonstrate that pre-estimating model parameters can improve the detection of sudden minor changes in tissue parameters, we simulated a 100 × 100 px skin phantom with multiple skin nevi. Specifically, the majority of skin was generated from baseline spectra (Spectrum 2 in Fig. 2), whereas skin nevi were generated from Spectra 3, all with added 30 dB noise to simulate skin and detector inhomogeneity. Therefore, the tissue properties outside and inside the nevi were identical except for ${f_\textrm{m}}$, which was on average 2% and 5% inside and outside nevi, respectively. Three circular skin nevi were generated with different radii (1 px, 5 px and 10 px). Figure 16(a) shows the ${f_\textrm{m}}$ map extracted using fitting case a (no fixed parameters) and the first approach, whereas Fig. 16(b) displays the ${f_\textrm{m}}$ map extracted using fitting case a and the second approach with parameter pre-estimation.
Fig. 16. ${f_\textrm{m}}$ maps extracted using fitting case a (no fixed parameters) and a) the first or b) the second approach using parameter pre-estimation.
We can see that the smallest nevus is invisible in the first colour map but discernible in the second, proving that parameter pre-estimation indeed improves the detectability of sudden minor changes in the tissue. The mean values of ${f_\textrm{m}}$ outside and inside the nevi in Fig. 16(a) were 1.60% ± 0.29% and 2.69% ± 0.60%, respectively. The respective mean values of ${f_\textrm{m}}$ outside and inside the nevi in Fig. 16(b) were 2.14% ± 0.06% and 4.90% ± 0.14%. Comparing these values demonstrates that parameter pre-estimation results in more accurate and robust estimations. Specifically, the mean values in Fig. 16(b) agree more with the actual values of 2% and 5% than those in Fig. 16(a)), and mean ${f_\textrm{m}}$ values in Fig. 16(b) are less scattered than those in Fig. 16(a) – inside and outside the nevi. We also experimented with adding less noise to the spectra, but the improvement in ${f_\textrm{m}}$ homogeneity was less pronounced, while the accuracy was higher in the pre-estimated map, similar to the high noise case presented here.
All in all, results suggest that fitting hyperspectral images with multiple free parameters could be improved by pre-estimating them. At first glance, parameter pre-estimation might seem to prolong the time needed to fit hyperspectral images. However, we noticed an improved convergence of the IAD algorithm, resulting in an overall much faster fitting process: the mean number of iterations was reduced by 11, and the mean time per spectrum was lowered on average by 0.08 s.
Furthermore, the values of estimated optical parameters for the second approach presented in Table 7 agree between different fitting cases and the physiological values reported in the literature. Specifically, the respective mean values of ${f_\textrm{m}}$ are all within the acceptable range of 0.1-15% stated by Verdel et al. [16], and the respective mean values of ${f_\textrm{Hb}}$ and ${f_{\textrm{Hb}{\textrm{O}_2}}}$ also fall within the permitted range of 0.1-20% reported in [16]. As a reference, the calculated mean tissue oxygen saturation ($\textrm{St}{\textrm{O}_2}$) values coincide with ranges reported by Fredriksson et al. [14] and Bruins et al. [44]. The mean values of scattering coefficient, $as$, are generally also within the range of 29.7-68.7 [4]. It can be seen from Table 7 that the IAD algorithm can successfully differentiate between different phases of the arterial occlusion test based on the values of tissue oxygenation. In both fitting cases, the initial tissue oxygenation dramatically decreases during the arterial occlusion when blood flow is obstructed, followed by a significant increase after the blood pressure cuff is deflated due to reflow of highly oxygenated blood in the examined area.
Table 7. Mean values and standard deviations of optical parameters estimated using case a or case z (second approach) for all participants. The values in the table are presented for each phase of the vascular occlusion test separately.
Further analysis of cases a and z for both fitting approaches reveals that the cumulative RMSE of fitted spectra in the second approach is higher in case z ($3.49\cdot {10^{ - 3}}$) than case a ($1.06\cdot {10^{ - 3}}$) – a consequence of a worse agreement of fitted and measured spectra due to fewer degrees of freedom of the model in case z than case a. We can see that the mean values and standard deviations of estimated parameters presented in Table 7 are similar for both fitting cases, likely as a result of parameter pre-estimation. Another possible explanation why lower standard deviations of parameters in case z are not observed is because fixed parameters in case z are set to incorrect values, forcing variable parameters to compensate their contribution to improve the agreement of fitted and measured skin spectra.
This observation reveals the study’s primary limitation and raises an essential question on how fixing model optical parameters to incorrect values affects the accuracy of the fitting process. Although it was not considered in this study, a significant deterioration in fitting performance is observed. Typically, optical parameters such as scattering coefficient, scattering power and layer thicknesses are adopted from the literature. However, different studies report similar yet discrepant values for various optical parameters. For example, the values of the scattering coefficient of the entire skin reported by Jacques [4] ranged from 29.7 [1] to 60.1 [45] for the whole skin when Rayleigh scattering is disregarded, a twofold increase in value. Therefore, authors usually choose arbitrary values within the reported physiological values of skin optical parameters, potentially leading to significant errors in free parameter estimation. A feasible yet exhausting approach to studying this topic would be to systematically fit simulated skin spectra by fixing various optical parameters to incorrect values and examining the effect on the estimation of free parameters. However, careful attention must be paid to how incorrect values of model parameters are selected. A possible option is to uniformly select the values of parameters within the physiological values reported in the literature. Nevertheless, this could be troublesome due to the scarcity of reported values for particular parameters such as cytochrome C oxidase.
We also recognise there was only one fitting case (case a) with all free parameters and one with seven fixed parameters (case zz). We discussed that as the number of fixed parameters increases from 4 to 7 in Fig. 7, the normalised cumulative RMSE and MAE values of fitted parameters and spectra decrease. However, this is true for the average normalised cumulative RMSE and MAE values calculated with respect to the number of fixed parameters, but for 0 and 7 fixed parameters, only one point was used to calculate the average. As a result, some discrepancies occur for particular fitting cases; for example, ${f_\textrm{RMSE}}({{w_1},{w_2}} )$ in case z is lower than in case zz for some combinations of ${w_1}$ and ${w_2}$ (see Fig. 5(a)), although fewer parameters are fixed.
A question might arise why not fix all possible combinations of parameters instead of the presented 27 combinations. With the clinical utility of hyperspectral imaging to diagnose skin diseases being our primary goal, there is a need to estimate parameters such as ${f_\textrm{m}}$, ${f_\textrm{Hb}}$, ${f_{\textrm{Hb}{\textrm{O}_2}}}$ and ${f_\textrm{brub}}$ and not fix them. Parameters ${f_\textrm{CO}}$ and ${f_{\textrm{CO}{\textrm{O}_2}}}$ were mainly introduced in the model to improve the spectra fitting at the 600–700 nm slope. Because both are present in the tissue simultaneously, it is sensible to either simultaneously consider them or disregard them both from the model. Moreover, epidermal and dermal thicknesses, ${d_e}$ and ${d_d}$, were fixed in all fitting cases except case a, where no parameters were fixed. Case a served as a reference to show the accuracy and robustness of the unrestricted parameter estimation method. Fixing just ${d_\textrm d}$ does not affect the spectra shapes significantly, as it is thick enough to be considered semi-infinite (${d_\textrm d}$ = 1 cm). We recognise that ${d_\textrm e}$, on the other hand, significantly affects the spectra and should be considered individually in the research. However, a strong negative correlation between ${d_\textrm e}$ and ${f_\textrm{m}}$ (–50 to –60%) was found (Fig. S2), prohibiting robust estimation of both parameters simultaneously. Therefore, when the parameter extraction was performed with both parameters free, the resulting maps of ${d_\textrm e}$ and ${f_\textrm{m}}$ were erratic. To increase the robustness, fixed ${d_\textrm e}$ and free ${f_\textrm{m}}$ were selected. Also, a recent study by Wang et al. [46] shows that ${f_\textrm{m}}$ multiplied by ${d_\textrm e}$ was estimated more accurately than ${f_\textrm{m}}$ itself.
In this particular study, no tissue-simulating phantoms with known optical properties such as epidermal thickness were used to verify the results. There are two main reasons for that; firstly, we are interested in extracting parameters from natural skin rather than tissue-mimicking phantoms; secondly, producing tissue phantoms with anatomy and optical properties representative of natural skin is exceptionally demanding. Thus, a compromise would have to be found between the production difficulty and representativeness of the phantoms. A limited complexity of such phantoms would require our tissue model to be simplified, eventually leading to unrepresentative verification. However, we acknowledge that time and effort should be put to verify our results on tissue-simulating phantoms in future work. Doing so would help us assess the real-world accuracy (with experimental setup uncertainties introduced) of the IAD method by providing us with invaluable reference values of optical properties. Nevertheless, our recent study extracted the concentrations of absorbing silicon pigment and scattering microspheres from a set of simple hemispherical tissue phantoms using a one-layer IAD method [47]. The extracted concentrations agreed well with the actual concentrations used in phantom preparation (less than 10% relative difference).
Our model could be disregarded as a simplified representation of natural human skin. Indeed, human skin is a complex heterogeneous tissue consisting of (at least) three distinct layers [1]. The dermis, for example, comprises two structurally different layers, specifically papillary and reticular dermis, with contrasting tissue and consequentially optical properties. Therefore, another possible improvement of this study is the model expansion, namely adding the third and perhaps fourth skin layer and other variable optical parameters (e.g., water, lipids). It would undoubtedly increase the complexity of our current model by adding more parameters and thus degrees of freedom, as seen in [3,6,14,23,24]. However, we believe that the agreement between the model output and the experimental spectra will not improve significantly by adding additional layers, as demonstrated by Bjorgan and Randeberg [24].
In his review, Jacques [4] demonstrated the effect of sequentially adding optical parameters such as melanin, water and lipids to a generic tissue on the total absorption coefficient. He showed that arbitrary skin spectra could be constructed from individual components by different skin tissue models. A growing number of papers attempt to characterise biological tissue such as human skin physiologically and structurally using tissue models with numerous model parameters [3,5–25,30,47,48]. Some model parameters (or sets of parameters) such as epidermal thickness are fixed to reduce the complexity of the fitting problem and improve the speed and convergence of the fitting algorithms. While multiple studies addressed the effect of using different tissue models on the extraction of optical properties [3,16,24], no research systematically studied the effect of fixing model parameters on the remaining variable model parameters. At best, Bjorgan and Randeberg [24] evaluated the deviation of layer thickness, melanin content and scattering coefficients when a single parameter was fitted in their IMC model, when all parameters were fitted and when depth and melanin were fitted. Less extensively, Liu and Ramanujam [11] studied the effect of epidermal thickness on the deviation of the total absorption and scattering coefficients estimated by their IMC model.
Similarly, Sharma et al. [13] verified their two-layer IMC model by comparing estimated top layer thickness, reduced scattering coefficient and absorption coefficient in the top and the bottom layer with the actual values. However, the latter two studies did not explore the fixation of particular optical parameters. Therefore, with this study, we remind the biomedical optics community of the influence of the parameter’s values in tissue modelling, a problem that is too often overlooked, and clarify the effect of fixation on the fitting results.
5. Conclusion
In this study, the performance of a two-layer GPU-accelerated IAD algorithm was examined on noiseless simulated skin spectra, realistic simulated skin spectra with three different levels of noise added, and in-vivo measured hyperspectral images of human hands recorded before, during, and after arterial occlusion. The results on the simulated skin spectra suggest that fixing a large number of model optical parameters – although highly subjective to the selection of parameters – generally improves the performance of the IAD algorithm, making it more accurate and robust. Adding noise to the simulated spectra negatively affects algorithms’ performance, yet the performance is satisfactory for the experimental noise level detected in our integrated HSI system – the lower the noise, the better the fitting performance. Moreover, the IAD method successfully demonstrated its potential in extracting optical properties of skin from in-vivo measured hyperspectral images. Fixing multiple parameters generally improves parameter estimation, but results suggest that fixing parameters to incorrect values deteriorates the algorithm's overall performance. However, pre-estimating parameters before actual fitting could improve the convergence of the IAD algorithm and speed up the fitting process.
This research work provides a foundation for studying the effect of fixing optical parameters on the performance of various iterative methods such as IMC and IAD. Therefore, it should be considered in subsequent studies where diffuse reflectance spectra are analysed to extract optical properties of the skin.
Funding
Javna Agencija za Raziskovalno Dejavnost RS (J3-3083, P1-0389).
Acknowledgements
The authors acknowledge the financial support from the Slovenian Research Agency (ARRS) (P1-0389, J3-3083).
Disclosures
The authors declare that there are no conflicts of interest related to this article.
Data availability
Data is not publicly available at the moment but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: a review,” J. Innov. Opt. Health Sci. 04(01), 9–38 (2011). [CrossRef]
2. S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid media by using the adding–doubling method,” Appl. Opt. 32(4), 559 (1993). [CrossRef]
3. R. Hennessy, M. K. Markey, and J. W. Tunnell, “Impact of one-layer assumption on diffuse reflectance spectroscopy of skin,” J. Biomed. Opt 20(2), 027001 (2015). [CrossRef]
4. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef]
5. E. Häggblad, T. Lindbergh, M. G. D. Karlsson, H. Casimir-Ahn, E. G. Salerud, and T. Strömberg, “Myocardial tissue oxygenation estimated with calibrated diffuse reflectance spectroscopy during coronary artery bypass grafting,” J. Biomed. Opt. 13(5), 054030 (2008). [CrossRef]
6. T. Lindbergh, E. Häggblad, H. Ahn, E. Göran Salerud, M. Larsson, and T. Strömberg, “Improved model for myocardial diffuse reflectance spectra by including mitochondrial cytochrome aa3, methemoglobin, and inhomogenously distributed RBC,” J. Biophoton. 4(4), 268–276 (2011). [CrossRef]
7. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37(4), 779 (1998). [CrossRef]
8. S. L. Jacques, “Optical assessment of cutaneous blood volume depends on the vessel size distribution: a computer simulation study,” J. Biophoton. 3(1-2), 75–81 (2010). [CrossRef]
9. S. A. Lisenko, M. M. Kugeiko, V. A. Firago, and A. N. Sobchuk, “Analytical model of diffuse reflectance spectrum of skin tissue,” Quantum Electron. 44(1), 69–75 (2014). [CrossRef]
10. L. Vidovic and B. Majaron, “Elimination of single-beam substitution error in diffuse reflectance measurements using an integrating sphere,” J. Biomed. Opt. 19(2), 027006 (2014). [CrossRef]
11. Q. Liu and N. Ramanujam, “Sequential estimation of optical properties of a two-layered epithelial tissue model from depth-resolved ultraviolet-visible diffuse reflectance spectra,” Appl. Opt. 45(19), 4776 (2006). [CrossRef]
12. G. Mantis and G. Zonios, “Simple two-layer reflectance model for biological tissue applications,” Appl. Opt. 48(18), 3490 (2009). [CrossRef]
13. M. Sharma, R. Hennessy, M. K. Markey, and J. W. Tunnell, “Verification of a two-layer inverse Monte Carlo absorption model using multiple source-detector separation diffuse reflectance spectroscopy,” Biomed. Opt. Express 5(1), 40 (2014). [CrossRef]
14. I. Fredriksson, “Inverse Monte Carlo method in a multilayered tissue model for diffuse reflectance spectroscopy,” J. Biomed. Opt 17(4), 047004 (2012). [CrossRef]
15. I. Fredriksson, O. Burdakov, M. Larsson, and T. Strömberg, “Inverse Monte Carlo in a multilayered tissue model: merging diffuse reflectance spectroscopy and laser Doppler flowmetry,” J. Biomed. Opt 18(12), 127004 (2013). [CrossRef]
16. N. Verdel, A. Marin, M. Milanič, and B. Majaron, “Physiological and structural characterization of human skin in vivo using combined photothermal radiometry and diffuse reflectance spectroscopy,” Biomed. Opt. Express 10(2), 944 (2019). [CrossRef]
17. M. Kirillin, I. Meglinski, V. Kuzmin, E. Sergeeva, and R. Myllylä, “Simulation of optical coherence tomography images by Monte Carlo modeling based on polarization vector approach,” Opt. Express 18(21), 21714 (2010). [CrossRef]
18. I. V. Meglinskii, “Monte Carlo simulation of reflection spectra of random multilayer media strongly scattering and absorbing light,” Quantum Electron. 31(12), 1101–1107 (2001). [CrossRef]
19. I. V. Meglinsky and S. J. Matcher, “Modelling the sampling volume for skin blood oxygenation measurements,” Med. Biol. Eng. Comput. 39(1), 44–50 (2001). [CrossRef]
20. I. V. Meglinski and S. J. Matcher, “Computer simulation of the skin reflectance spectra,” Computer Methods and Programs in Biomedicine 70(2), 179–186 (2003). [CrossRef]
21. A. P. Popov, A. V. Bykov, and I. V. Meglinski, “Influence of probe pressure on diffuse reflectance spectra of human skin measured in vivo,” J. Biomed. Opt 22(11), 1 (2017). [CrossRef]
22. D. Y. Churmakov, I. V. Meglinski, S. A. Piletsky, and D. A. Greenhalgh, “Analysis of skin tissues spatial fluorescence distribution by the Monte Carlo simulation,” J. Phys. D: Appl. Phys. 36(14), 1722–1728 (2003). [CrossRef]
23. Hanna Karlsson, Anders Pettersson, Marcus Larsson, and Tomas Strömberg, “Can a one-layer optical skin model including melanin and inhomogeneously distributed blood explain spatially resolved diffuse reflectance spectra?” in (2011), Vol. 7896.
24. A. Bjorgan and L. L. Randeberg, “Exploiting scale-invariance: a top layer targeted inverse model for hyperspectral images of wounds,” Biomed. Opt. Express 11(9), 5070 (2020). [CrossRef]
25. H. Jonasson, I. Fredriksson, S. Bergstrand, C. J. Östgren, M. Larsson, and T. Strömberg, “In vivo characterization of light scattering properties of human skin in the 475- to 850-nm wavelength range in a Swedish cohort,” J. Biomed. Opt. 23(12), 1 (2018). [CrossRef]
26. L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Computer Methods and Programs in Biomedicine 47(2), 131–146 (1995). [CrossRef]
27. Vinoin Devpaul Vincely and Karthik Vishwanath, “Lateral light losses in integrating sphere measurements: comparison of Monte-Carlo with inverse adding-doubling algorithm,” in (2020), Vol. 11231.
28. A. Sterle, “Implementacija algoritma Adding-Doubling v programskem jeziku Python,” 61 (n.d.).
29. O. Nelles, Nonlinear System Identification (Springer Berlin Heidelberg, 2001). [CrossRef]
30. L. Rogelj, U. Pavlovčič, J. Stergar, U. Simončič, and M. Milanič, “Correction of hyperspectral images with build- in 3D laser profilometry,” 16 (n.d.).
31. G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” J. Biomed. Opt 19(1), 010901 (2014). [CrossRef]
32. J. W. Pickering, S. A. Prahl, N. van Wieringen, J. F. Beek, H. J. C. M. Sterenborg, and M. J. C. van Gemert, “Double-integrating-sphere system for measuring the optical properties of tissue,” Appl. Opt. 32(4), 399–410 (1993). [CrossRef]
33. M. S. Wróbel, A. P. Popov, A. V. Bykov, M. Kinnunen, M. Jędrzejewska-Szczerska, and V. V. Tuchin, “Multi-layered tissue head phantoms for noninvasive optical diagnostics,” J. Innov. Opt. Health Sci. 08(03), 1541005 (2015). [CrossRef]
34. Paul Lemaillet, Jean-Pierre Bouchard, Jeeseong Hwang, and David W. Allen, “Double-integrating-sphere system at the National Institute of Standards and Technology in support of measurement standards for the determination of optical properties of tissue-mimicking phantoms,” J. Biomed. Opt. 20, 1–8 (2015). [CrossRef]
35. P. Lemaillet, C. C. Cooksey, J. Hwang, H. Wabnitz, D. Grosenick, L. Yang, and D. W. Allen, “Correction of an adding-doubling inversion algorithm for the measurement of the optical parameters of turbid media,” Biomed. Opt. Express 9(1), 55–71 (2018). [CrossRef]
36. D. Hu, R. Lu, Y. Huang, Y. Ying, and X. Fu, “Effects of optical variables in a single integrating sphere system on estimation of scattering properties of turbid media,” Biosystems Engineering 194, 82–98 (2020). [CrossRef]
37. E. K. Chan, B. Sorg, D. Protsenko, M. O’Neil, M. Motamedi, and A. J. Welch, “Effects of compression on soft tissue optical properties,” IEEE J. Sel. Top. Quantum Electron. 2(4), 943–950 (1996). [CrossRef]
38. J. F. Beek, P. Blokland, P. Posthumus, M. Aalders, J. W. Pickering, H. J. C. M. Sterenborg, and M. J. C. van Gemert, “In vitrodouble-integrating-sphere optical properties of tissues between 630 and 1064 nm,” Physics in Medicine and Biology 42(11), 2255–2261 (1997). [CrossRef]
39. Dhiraj K. Sardar, Felipe S. Salinas, John J. Perez, and Andrew T. Tsin, “Optical characterization of bovine retinal tissues,” J. Biomed. Opt. 9(3), 624–631 (2004). [CrossRef]
40. Dhiraj K. Sardar, Raylon M. Yow, Andrew Tsin, and Ratna Sardar, “Optical scattering, absorption, and polarization of healthy and neovascularized human retinal tissues,” J. Biomed. Opt. 10, 1–8 (2005). [CrossRef]
41. S. C. Gebhart, W. C. Lin, and A. Mahadevan-Jansen, “In vitrodetermination of normal and neoplastic human brain tissue optical properties using inverse adding-doubling,” Physics in Medicine and Biology 51(8), 2011–2027 (2006). [CrossRef]
42. D. K. Sardar and L. B. Levy, “Optical properties of whole blood,” Lasers Med. Sci. 13(2), 106–111 (1998). [CrossRef]
43. D. K. Sardar, M. L. Mayo, and R. D. Glickman, “Optical characterization of melanin,” J. Biomed. Opt. 6(4), 404 (2001). [CrossRef]
44. A. A. Bruins, D. G. P. J. Geboers, J. R. Bauer, J. H. G. M. Klaessens, R. M. Verdaasdonk, and C. Boer, “The vascular occlusion test using multispectral imaging: a validation study: the VASOIMAGE study,” J Clin Monit Comput 35(1), 113–121 (2021). [CrossRef]
45. I. S. Saidi, S. L. Jacques, and F. K. Tittel, “Mie and Rayleigh modeling of visible-light scattering in neonatal skin,” Appl. Opt. 34(31), 7410 (1995). [CrossRef]
46. C.-Y. Wang, T.-C. Kao, Y.-F. Chen, W.-W. Su, H.-J. Shen, and K.-B. Sung, “Validation of an inverse fitting method of diffuse reflectance spectroscopy to quantify multi-layered skin optical properties,” Photonics 6(2), 61 (2019). [CrossRef]
47. L. Rogelj, U. Simončič, T. Tomanič, M. Jezeršek, U. Pavlovčič, J. Stergar, and M. Milanič, “Effect of curvature correction on parameters extracted from hyperspectral images,” J. Biomed. Opt. 26(09), 096003 (2021). [CrossRef]
48. Maria Ewerlöf, E. Göran Salerud, Tomas Strömberg, and Marcus Larsson, “Estimating skin blood saturation by selecting a subset of hyperspectral imaging data,” in (2015), Vol. 9328.
