## Abstract

Detection of re-epithelialization in wound healing is important, but challenging. Hyperspectral imaging can be used for non-destructive characterization, but efficient techniques are needed to extract and interpret the information. An inverse photon transport model suitable for characterization of re-epithelialization is validated and explored in this study. It exploits scale-invariance to enable fitting of the epidermal skin layer only. Monte Carlo simulations indicate that the fitted layer transmittance and reflectance spectra are unique, and that there exists an infinite number of coupled parameter solutions. The method is used to explain the optical behavior of and detect re-epithelialization in an in vitro wound model.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

In vitro wound models are useful for investigation of wound healing in a controlled laboratory setting [1–5]. However, it is challenging to monitor the actual healing in such models without destructive histology analysis. Hyperspectral imaging is a technique providing a non-destructive, objective method for characterizing such tissues optically.

Hyperspectral imaging has a spectral and spatial resolution that has been shown to be useful for biomedical applications like wound imaging [6–14], burn wound imaging [12,15], cancer diagnostics [6,16,17] and surgical guidance [6,18]. Statistical processing techniques are often used to handle the large amounts of data [8,12,13,15,19–25]. The rich spectral content further enables use of inverse photon transport modeling [26–32]. Such modeling techniques can be used to interpret the data and relate spectral changes to changes in skin properties like blood content and blood oxygenation through constrained fitting of optical properties. This indirect way of estimating optical properties is considered to be ill-defined [6,33] since different media can have similar reflectance spectra [33,34]. Further, absorption and scattering properties are fitted to a single reflectance measurement [35]. A priori knowledge of the expected shapes of the absorption and scattering restrains the problem somewhat [35], but there is still a basic ambiguity resulting from a scale-invariance of the reflectance with respect to the absorption and scattering spectra [33,34]. Other possibilities to restrain such models can be a valuable road of study in order to obtain unique and robust estimates.

The tissue model used in this study was based on an in vitro wound model setup developed by Jansson et al. [36] and Kratz [37]. Here, samples with wounds are prepared from ex vivo human tissue and placed in a growth medium, which causes the wounds to re-epithelialize. Characterization of the re-epithelialized layer is of special interest, e.g. the presence, maturity and thickness of the layer. Photon transport modeling techniques could be used to extract this information from the reflectance spectra.

Photon transport modeling of in vitro wounds can be challenging. The spectral characteristics of the dermal part of the wound models are less defined than in normally circulated tissue due to the lack of blood absorption, and the spectra are significantly influenced by the spectral properties of the growth medium [38]. However, in wounds, both regions with and without the upper tissue layer are present. The spatial resolution of hyperspectral imaging makes multiple reflectance spectra of each tissue type available as every pixel contains data with high spectral resolution. This could potentially be used to restrain the model and target the epidermal layer specifically.

The main idea of the inverse modeling technique to be presented in this paper is that the basic scale-invariant limitations of the reflectance modeling can be exploited to enable consideration of the upper layers without having to completely model the lower layers. The availability of reflectance spectra from both wound and re-epithelialized tissue then enables quantification of skin properties in re-epithelialized layers, without the need to consider the dermal layers.

The basic method requires knowledge of an appropriate wound spectrum. The high number of available wound spectra in a hyperspectral image makes selection of one unique spectrum somewhat prohibitive. However, the same availability of multiple spectra with a high spectral resolution makes it possible to use dimensionality reduction methods to remove redundant information and represent the spectral information in a low-dimensional space [38]. Principal component analysis (PCA) is such a method, which can decompose a dataset in terms of orthonormal components (loadings) that can linearly transform each observation into new variance-maximizing coordinates (scores) [39]. This technique has been used as a pre-processing technique [40–42] and for investigation of spectra in a low-dimensional space [6,19,38]. The method is used in the current paper to reduce a discrete wound spectrum choice to a continuous choice of PCA scores that can be back-transformed to a wound-like spectrum using the inverse transform. This enables the wound spectrum selection to be an efficiently evaluated part of the model optimization. The main goal of this study is to validate a proof of concept of the presented basic inverse modeling technique including the suggested PCA extension, and explore the limitations and possibilities of this approach.

A simulation study is carried out using Monte Carlo simulations in MCML [43], which represents a gold standard for photon transport simulations. The simulation study is used to investigate and verify the assumptions that enable use of the inverse modeling method. The uniqueness and accuracy of the fitted skin parameters are explored. With the findings of the simulation study on uniqueness and parameter accuracy in place, the method is finally applied to experimental data. Hyperspectral images of an in vitro wound model sample are used as an example for the application of the technique.

The inverse modeling method and its basic assumptions are given in section 2.1, along with the PCA-based modification appropriate for wound imaging. The simulation study used to verify these assumptions and investigate the accuracy and uniqueness of the inverse modeled parameters is outlined in section 2.2. Information on the hyperspectral wound dataset used to demonstrate an application of the method is given in section 2.3. Finally, results and discussion are given in section 3.

The method represents a partial modeling approach which combines data-driven results with photon transport modeling. Model fitting is mainly constrained to the upper layer, alleviating some of the inherent undetermined nature of the inverse photon transport modeling approach. The example given here is wound healing, but the method is generic, and can be used to investigate any top layer given that reflectance spectra are available with and without modifications to or removal of the top layers. Examples include characterization of burn wounds and estimation of epidermal skin thickness in pre-term newborns.

## 2. Materials and methods

#### 2.1 Inverse modeling setup

### 2.1.1 Fundamental assumptions

Two main assumptions form the basis of the method:

- 1. The reflectance of a one-layer model can be written as a function of $\mu _a/\mu _s'$ (scale-invariance).
- 2. A multi-layer model and a one-layer representation yield identical reflectance values when a top layer is added.

### 2.1.2 Photon transport model

A diffusion model with an isotropic source function [44] is used as photon transport model. The main advantage of this model is its simple, analytic expression for the reflectance, enabling fast evaluation during optimization. In addition, it can yield $\mu _a/\mu _s'$ directly from the reflectance without iteration.

### 2.1.2.1 Theory

Photon transport in biological tissue can be modeled by the radiative transfer equation (RTE) [45]. Assuming an almost isotropic light distribution and isotropic source functions, the time-independent RTE in a one-dimensional geometry is simplified to [44,45]

where the diffusion constant is $D = \frac {1}{3(\mu _s' + \mu _a)}$ and $\phi$ is the fluence rate. A multi-layer medium is assumed, with $d_i$ describing the depth of layer interface $i$. With the source function in a layer $i$ given as [44] the solution for the fluence rate in the layer $i$ is given as [44]### 2.1.2.2 Offset correction

A correction constant is applied to the diffusion model reflectance. Comparing the diffusion model with the same boundary conditions and optical properties as a Monte Carlo simulation shows that the output reflectance from the diffusion model has an offset [46,47]. The assumption of an almost isotropic light distribution in the diffusion model leads to a less forward-directed photon flux close to the surface, as compared to other source functions such as the Delta-Eddington source function [47], and this yields a higher output reflectance contributing to the observed offset. The isotropic source function is chosen for simplicity and convenience, and it is considered as out of scope of this proof of principle to minimize this well known and systematic offset. It is however acknowledged that it introduces systematic errors in the estimated parameters.

The diffusion model and two-layer Monte Carlo spectra from model A and model C1 to be described in section 2.2.1 were compared with equal input parameters. An offset correction of 0.036 was found to minimize the average root mean squared error (RMSE) among all spectra. The offset correction is demonstrated in Fig. 2. On average, the RMSE between the model C1 Monte Carlo spectra in section 2.2.1 and corresponding diffusion model spectra were 0.037 (standard deviation 0.003) and 0.005 (standard deviation 0.001) before and after offset correction, respectively.

### 2.1.3 Inverse modeling method for skin

The main application of the method considers human skin with and without epidermis present, such as in wounds. Construction of a two-layer model from a basis reflectance was outlined in section 2.1.1. The properties of epidermis are then fitted to the reflectance from the region with the top layer. Python was used for development of the technique.

Melanin is assumed to be the main absorber in epidermis in the visible range [45,48,49]. Minor absorbers include carotene, lipids, cell nuclei and filamentous proteins [50], and are often modeled using a bulk background absorption [31,44,51,52]. A small amount of blood is sometimes included in the epidermis to account for the non-planar geometry of the papillary dermis [44,51,52]. For simplicity, neither of these are included in the current model. The epidermis is assumed to consist of a single layer with melanin as the absorber as shown in Eq. (6) and a scattering corresponding to Eq. (8), along with a defined layer thickness.

The objective function to be minimized was chosen to be the RMSE between measured reflectance and simulated reflectance, relative to the simulated reflectance. The function `minimize` from the Python package `scipy.optimize` was used to minimize the objective function, using L-BFGS-B (Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm with box constraints) as the optimization method. The parameters were bounded, and they were rescaled to values between 0 and 1 in order to make them comparable. The fitted epidermal skin parameters are listed in Table 1.

### 2.1.4 Modifications to the technique for application to wounds

A modification of the inverse model was used for hyperspectral images of wounds, as the appropriate basis spectrum for a given re-epithelialized tissue sample is not known a priori. It was desired to let selection of the wound spectrum be something which could be fitted during optimization rather than iterating through the possible choices.

A PCA transform with three components is applied to spectra from the wound region. This reduces all possible wound spectra down to a combination of three score parameters and corresponding PCA loading vectors. Three additional parameters were input during the optimization. These were used as PCA scores and inverse transformed to construct an artificial wound basis spectrum onto which an epidermis was placed. This allowed the inverse model to represent a wound spectrum that could be fitted during optimization.

The PCA scores were then fitted along with the rest of the parameters. The fitted scores were bounded within the min/max range of the scores as transformed from the original wound spectra. Reconstructed wound basis spectra were constrained to non-negative values, effectively changing the optimization method to SLSQP (Sequential least squares programming).

#### 2.2 Simulation study

### 2.2.1 Simulation setup

GPU-MCML [53,54] was used to simulate reflectance spectra from a skin-like geometry. NVIDIA GeForce GTX 670 was used for GPU parallelization. Wavelengths from $\lambda = 400$ to 850 nm were modeled with 3 nm discretization. 1 000 000 photons were used in each simulation. The total run time for a full spectrum was between 12 and 22 seconds. All simulations were run with pencil beams incident on the skin model, and the tissue was assumed to have a refraction index of 1.4. The homogeneity of the model makes the total integrated reflectance equivalent with the reflectance from the same model illuminated with an infinitely broad beam. A layer thickness of 1 meter was used in order to emulate a semi-infinite layer.

### 2.2.1.1 Absorption properties

The absorption in the top layer was modeled using an absorption model for melanin [55]

where $\mu _{\textrm {a,m,694}}$ is a parameter associated with the mean melanin content of the layer. Moderately dark skin corresponds to an absorption in the range 500-900 m$^{-1}$, while fair skin corresponds to an absorption in the range 200-300 m$^{-1}$ [44]. Assumed low and high values for melanin absorption and epidermal thickness are listed in Table 2. For the dermal layers, the absorption is modeled as### 2.2.1.2 Scattering properties

For scattering, the following model is used:

### 2.2.1.3 Model geometries

Geometries used in this paper are shown in Fig. 3: One-layer geometry, multi-layer geometry and two-layer geometry. The following simulations were run:

- • Model A (one-layer model): All combinations of deeper layer parameters in Table 2 were used, yielding in total 16 varieties.
- • Model B (multi-layer model): Deeper layer parameters in Table 2 were randomly picked in each of the layers in a three-layer model. Layer thicknesses were set to 150 µm. The simulation was run with and without an additional top layer with thickness 100 µm, melanin content 300 m$^{-1}$ and the low scattering parameters in Table 2.
- • Model C1 (systematic two-layer model): Lowest deeper layer parameters in Table 2 were used for the deeper layer. The thickness of the top layer was varied between 100, 150, 200, 250 and 300 µm, the melanin content between 150, 300 and 700 m$^{-1}$, and all combinations of the scattering values listed in Table 2 were used. This yielded in total 60 varieties, or 12 spectra per layer thickness step.
- • Model C2 (random two-layer model): The properties of the deeper layer were randomly selected among the deeper layer parameters listed in Table 2. Parameters corresponding to the lowest scattering values were selected for the top layer, while melanin in (6) and thickness were randomly selected from a uniform probability distribution bounded by the lower and upper values in Table 2. 50 model varieties were sampled.

### 2.2.2 Verification of modeling assumptions

Two main assumptions for the technique to be applicable were given in section 2.1.1.

The model A simulations (one-layer models) above yield, over the various wavelengths, a large range of possible combinations of $\mu _a$ (65 to 6714 m$^{-1}$) and $\mu _s'$ (1517 to 5237 m$^{-1}$). The assumption that the reflectance can be written as a function of $\mu _a/\mu _s'$ is checked by plotting the output reflectance values as a function of $\mu _a/\mu _s'$ across all simulations.

The model B simulations (multi-layer models) yield the reflectance from a complex multi-layer model with and without a top layer. An empirical Monte Carlo model was constructed for looking up $\mu _a/\mu _s'$ given a reflectance value through the use of a Savitzky-Golay fit and an interpolating natural cubic spline. This model was used to find the $\mu _a/\mu _s'$ ratio for a one-layer model from the complex multi-layer reflectance. A $\mu _s'$ was assumed ($\mu _s' = 1000 {\textrm{m}^{-1}} $), and a $\mu _a$ was calculated from the ratio. The assumption that a multi-layer model with an extra top layer is indistinguishable from a fitted one-layer model with the same top layer is then checked by comparing the latter model with an extra top layer to the original model with the same top layer.

### 2.2.3 Uniqueness of the inverse model solution

Diffusion model simulations were run in order to investigate whether multiple optical parameters could yield the same $R$, and find the shape of the relation among the optical properties, if any. A reflectance ($R = 0.6$) was picked, and a $\mu _a/\mu _s'$ ratio (1/200) was set in the deeper layer of a two-layer model. Top layer thicknesses ranging from 10 to 500 µm and scattering coefficients ranging from 10 to 100 000 m$^{-1}$ were put into the top layer. The absorption coefficient necessary to yield the assumed reflectance were derived for each parameter combination.

The uniqueness of the obtained skin parameter solution was then investigated. The inverse model was run repeatedly on the same simulated spectrum from model C1 with randomly generated start parameters in order to see whether it was possible to reach the same global solution.

### 2.2.4 Accuracy of the inverse modeled parameters

In addition to investigation of the uniqueness of the solution, it is desired to determine the accuracy of the inverse modeled parameters as compared to the input parameters in the Monte Carlo model. Three variations of the inverse model were run on the model C1 simulations:

- 1. Fit a single parameter, with all parameters except for one fixed. This estimates a baseline accuracy of the inverse model for each parameter.
- 2. Fit all parameters simultaneously.
- 3. Fit only thickness and melanin content, with scattering parameters fixed.

- 1. Set the scattering parameter to the true parameter, estimate melanin content and thickness. Evaluates the inverse model error when the scattering is known.
- 2. Set the scattering parameter to the results from multi-parameter optimization on the first spectrum, estimate only melanin content and thickness for the rest. This outlines the estimated parameter behavior when the scattering is unknown but known to be homogeneous.
- 3. Set the scattering parameter to the true parameter and use the PCA inverse model in section 2.1.4 to fit the rest of the parameters. Here, all spectra without epidermis were used to fit a PCA transform, and the PCA transform is used to find a best basis spectrum for the spectrum at hand during optimization.

#### 2.3 Experimental tests

**Hyperspectral acquisition** Reflectance data were acquired using a push-broom Hyspex VNIR-1600 hyperspectral camera (Norsk Elektro Optikk, Lillestrom, Norway). The images were acquired over the wavelength range 400-1000 nm, with a spectral resolution of 3.7 nm and an integration time of 7.5 ms per line of data. The camera has been radiometrically and spectrally calibrated using light sources with known characteristics. The radiometric calibration was used to apply correction factors to the images.

The reflectance data were acquired with illumination from two linear light sources (Model 2900 Tungsten Halogen, Illumination Technologies, New York). Polarizers (VLR-100 NIR, 450-1100 nm, Meadowlark Optics, Frederick, Colorado) were mounted on the camera lens and the light sources in order to avoid specular reflection. A Spectralon reflectance target (WS-1-SL, Ocean Optics, Duiven, Netherlands) was included within each image and used to convert the raw data to reflectance spectra.

**Wound model** Samples of the in vitro wound model were prepared from human abdominal skin. The project was approved by the regional ethical committee (REK-Midt-Norge), and informed consent was obtained from the donor. The sample used for demonstration in this study was prepared with a 4 mm wound using punch biopsy, and the full tissue sample was cut using 8 mm punch biopsy. The sample was incubated for 22 days in Dulbecco’s Modified Eagle Medium (Gibco, USA), with fetal calf serum (10%), penicillin (50 ug/ml), streptomycin (50 U/ml) and glutamine added. Hyperspectral images were acquired at day 1, 2 and then every other day, and the medium was changed every imaging session. The wound was exposed to air by resting the sample on a metallic grid, in order to ensure development and migration of multiple cell layers [36].

Re-epithelialization visible by visual inspection of the RGB images occurred during the last ten days. The sample was therefore investigated at day 12, 18 and 22.

A larger image subset over the wound boundary was selected at approximately the same region across these three days. PCA transforms were fitted to a subset of wound spectra at each day that by visual inspection and comparison of the spectra had no obvious re-epithelialization present. Each wound subset consisted of 7400 spectra, and 3 PCA components were used (explaining 87.9% of the variance, on average 0.012 RMSE between raw and reconstructed reflectance spectrum when running forward and inverse transforms). The modified PCA-based inverse model in section 2.1.4 was then run on each pixel in the larger subsets in order to yield skin parameters for the top layer.

## 3. Results and discussion

Simulations are presented in section 3.1. The modeling assumptions that form the basis of the inverse modeling technique are investigated using Monte Carlo simulations in section 3.1.1. The uniqueness and accuracy of inverse modeled parameters as compared to Monte Carlo simulations are given in 3.1.2 and 3.1.3. The simulation results are then summarized and discussed in section 3.1.4. The technique is used to estimate re-epithelialized layer thickness from hyperspectral images of wounds in section 3.1.5, and the performance of the technique is discussed in light of the findings from the simulation study.

#### 3.1 Simulation study

### 3.1.1 Verification of modeling assumptions

Two assumptions were given in section 2.1.1. These were that the reflectance of a one-layer model can be written as a function of $\mu _a/\mu _s'$, and that a multi-layer geometry with a top layer has a reflectance indistinguishable from a one-layer representation with the same top layer.

The validity of the first assumption is confirmed in Fig. 4. The figure shows reflectance spectra acquired across a wide range of one-layer models, and a corresponding plot over the same reflectance values as a function of $\mu _a/\mu _s'$. The latter clearly demonstrates that there is a one-to-one correspondence between the $\mu _a/\mu _s'$ ratio and the reflectance. Other studies also confirm this fact [34].

The second assumption was that a one-layer representation of a multi-layer model yields identical reflectance to the multi-layer model when a top layer is added to either. This is demonstrated in Fig. 5. Here, the reflectance from one-layer models constructed from each layer of the multi-layer model are used to verify that none of the upper layers completely shield the deeper layers. Further, the reflectance from the multilayer model with an epidermis on top is compared to a one-layer representation with the same epidermis on top. As they are indistinguishable, the second assumption is verified. The various optical properties at the different wavelengths represent a wide range of layer combinations that all demonstrate the validity of the second assumption. This is also mostly given by the fact that placing a layer on top of some existing model will not modify the existing parts of the model. Here, the RMSE between the two example spectra was 0.0010. Further testing the same assumption on 20 randomly generated multi-layer skin models yielded an RMSE of 0.0015.

The consequence of the two assumptions is that a top layer can be investigated without having to fully model the deeper layers, if measurements with and without the top layer are available.

### 3.1.2 Uniqueness of the inverse modeled solution

The combined results above means that the properties of the deeper layer are scale-invariant. This is a consequence of the output reflectance having no units [34]. A top layer is independent from the deeper layers, and a similar scale-invariant relation must exist between $\mu _a$, $\mu _s'$ and $d_1$ in order to produce a unit-less reflectance at the output. The epidermal skin parameters will therefore likely be coupled.

Parameter couplings between $d$ and $\mu _a$ for several $\mu _s'$ that yield the exact same reflectance are shown in Fig. 6. This demonstrates that there exists an infinite number of parameter sets that all yield the exact same reflectance, and sketches out the hyperplane on which the valid optical parameters reside.

To check potential skin parameter coupling, the model was fitted at random start parameters in Fig. 7. A clear relation between the estimated $d_1$ and the melanin content is evident. Fixing the scattering to the true parameters yields the same solution regardless of the start parameters. The relation is less clear between the estimated $d_1$ and the scattering parameters, but there are indications of a noisy quadratic relation similar to the melanin content relation. The scattering parameters were not found to influence the reflectance as much as the other parameters within the varied range, which can explain the noisiness of the relation. The lack of influence can be observed in Fig. 6, by the relations here being significantly changed only for larger scattering coefficients.

The RMSE shows that each of these solutions are identical with respect to the simulated reflectance. Fitting all unknown parameters at the same time is therefore not expected to yield a unique solution.

The found parameter non-uniqueness can be argued from the scale-invariance between $\mu _a$, $\mu _s'$ and $d_1$. The absorption coefficient $\mu _a$ is here modeled as a varying parameter multiplied by a fixed wavelength-dependency. Since the wavelength-dependency is fixed, it can be argued that the scale-invariance is translated into the parameter, and that this has a scale-invariant relation with $d_1$ and $\mu _s'$. The scattering model can be re-written to $\mu _s = A\left [f (\lambda /500)^{-b_{\textrm {Mie}}} + (1-f)(\lambda /500)^{-4}\right ]$, where $A = \mu _{\textrm {s,Mie,500}}' + \mu _{\textrm {s,Ray,500}}'$ and $f = \frac {\mu _{\textrm {s,Mie,500}}'}{A}$. Since $f$ is a dimensionless number and the wavelength-dependencies are fixed, the scale-invariant relation can be translated into $A$ and to the scattering parameters $\mu _{\textrm {s,Mie,500}}'$ and $\mu _{\textrm {s,Ray,500}}'$. A coupling between the skin parameters is therefore expected.

All of these solutions are valid with respect to the stated problem, as all yield the same fitted reflectance. It can be observed that they all apparently characterize a unique diffuse layer despite the large variation in skin parameters. Each of the valid parameter sets in Fig. 7 were used to model a single-layer model with finite thickness, and Monte Carlo simulations were used to obtain reflectance and transmittance spectra. The spectra are plotted in Fig. 8. These are more or less identical. Deviations can be attributed to albedo-dependent inaccuracies in identical diffusion model simulations that would lead to variations in a Monte Carlo simulation.

### 3.1.3 Accuracy of inverse modeled parameters

### 3.1.3.1 Systematic variation in input parameters (model C1)

The error deviation results across the different parameters are shown in Fig. 9 for three cases: Fit of the particular parameter only, fit of all parameters simultaneously, and fit of melanin content and layer thickness only. Averages over the corresponding reflectance errors between fitted and original spectra are shown in Fig. 10.

The case where all parameters are fitted simultaneously is first considered. Each fitted parameter deviates over a large range (highest relative error, low/high input parameter: thickness 69%/35%, melanin content 61%/40%, Mie scattering 116%/38%, Rayleigh scattering 93%/45%). This is to be expected due to the non-uniqueness of the solution. The variation in simulated reflectance likely trigger various local minima along the scale-invariance.

Fitting a single parameter is more interesting. The deviation range for thickness and melanin content is more limited for this case, due to the uniqueness of the solution (highest relative error, low/high input parameter: thickness 17%/14%, melanin content 14%/15%). The deviation of the scattering still matches the order of magnitude of the input, however (Mie scattering 122%/46%, Rayleigh scattering 61%/29%). Changing the epidermal scattering parameters do not change the reflectance as much as the absorption and layer thickness. The scattering can therefore not be reliably estimated, even if it is the only fitted parameter. The mismatch between the diffusion model and the Monte Carlo spectra at the true parameter leads to a bias in all parameters. Varying the particular parameter only does not bridge the mismatch entirely, as seen in Fig. 10 by the higher reflectance errors as compared to the cases where multiple parameters are fitted.

The parameter $b_{\textrm {Mie}}$ was not varied in the simulations. Recovered values (input 0.22) ranged from 0 to 3 in full parameter fit, and from 0 to 2 when it was the only parameter fitted.

Last, fixing scattering and fitting the rest of the parameters is considered. Fitting the thickness and melanin at the same time apparently yields a less biased estimate than when considering them apart. Further, the error is at most 80 m$^{-1}$ for larger melanin contents (highest relative error, low/high input: 12%/12%), while the thickness can be estimated down to an error below 14 µm (11%/5%).

### 3.1.3.2 Random selection of input parameters (model C2)

Parameter estimation results over random skin models are shown in Fig. 11. Using the true scattering yields a reasonable estimate of thickness and melanin content. The RMSEs are 17 µm (relative RMSE: 7%) and 57 m$^{-1}$ (8%) for the thickness and melanin contents, respectively. This is in line with the deviations found for thickness (below 14 µm) and melanin contents (below 80 m$^{-1}$) when fitting these simultaneously on the systematic variation earlier.

Exact knowledge of the correct scattering is challenging. It is expected that it must be guessed in the application at hand. A case where the scattering is estimated from a single reflectance spectrum and fixed for the rest of the spectra is shown in the same figure above. Both melanin content and thickness estimates become biased, but they remain correlated with the true parameters and retain variations that are similar to the case where the true scattering is used. Although the scattering is not known and the thickness estimates are incorrect, this can then still be used to detect relative thickness changes.

The method is to be applied to hyperspectral images of wounds, where the basis spectrum is not known a priori. The modification of the technique, as outlined in section 2.1.4, was used to fit basis spectra along with the rest of the epidermal skin parameters. The resulting RMSEs were 16 µm (relative RMSE: 6%) and 50 m$^{-1}$ (13%) for the thickness and melanin parameters, respectively, similar to RMSEs in the case where the basis spectra are known exactly.

### 3.1.4 Summary and general discussion of the simulation results

The inverse modeling method presented in this paper could be a valuable tool for characterizing hyperspectral images of re-epithelialized tissues in wounds.

The main inverse modeling assumptions have been verified. The reflectance of a one-layer model can be written as a function of $\mu _a/\mu _s'$. This means that any combination of $\mu _a$ and $\mu _s'$ yields the same reflectance as long as they obey the given ratio. Further, an arbitrary multi-layer model can be represented by a one-layer model. Adding a top layer to either of these models yields indistinguishable reflectance spectra. Thus, the top layer can be fitted and investigated without considering the deeper layers.

It has been shown that no unique solutions exist for the top layer. The solutions are coupled, however, and yield unique $R(\lambda )$ and $T(\lambda )$ through the upper layer that are common for all parameter sets. The fitted, diffuse layer is therefore unique, though the lack of a known thickness means that the optical properties are undeterminable. The wavelength dependency in $R$ and $T$ can be valuable for drawing conclusions about the nature of the layer.

The solutions for some of the skin parameters were found to be robust to changes in the start parameters during fitting when at least one parameter was fixated. This indicates that unique solutions may be possible to find in such cases. The reflectance was not found to be very sensitive to changes in the scattering parameter within the expected range. The scattering parameter is therefore a first choice for fixation. Fair estimates of the absorption parameter and layer thickness can be found when given the correct scattering, and the main expected levels can be discriminated. The method has been shown to yield reasonable relative estimates when the scattering is homogeneous, but unknown. The parameter value will then vary around a mean level within a small deviation, and be correlated with the true value. Such estimates are useful for determining whether a given location has a top layer thickness greater or less than the thickness of some other location. In practice, the scattering parameters can be fixed to e.g. the low parameter values outlined in section 2.2.1. Another possibility is to let all parameters be fitted simultaneously for a single spectrum in order to estimate a best fit for the wavelength-dependency, and then fixate the scattering parameters to these parameters for the rest of the spectra.

The method is thus useful for in vitro wounds in two ways. First by demonstrating whether the optical properties of various tissue regions can be explained by wound optical properties with an epidermal layer on top. Second by evaluating relative layer thicknesses at different positions, and further use these to explain spatial variations by layer thickness differences.

Melanin and thickness parameter fits have been found to have relative errors from 5 to 12%. Inverse methods in other studies that estimate epidermal thickness and melanin content are reported to have errors in the range of 9% for epidermal thickness and 8-15% for melanin content [58], or 6-8%, 16-20% for epidermal thickness and below 0.5% for melanin content [31], subject to modeling details. Relative errors of the current model are thus in the same range as methods reported in the literature.

A weakness of the method is that the basis spectrum representing the deeper layers must be known. While known exactly for the simulations, wounds have inhomogeneities that result in no clear basis candidate. Taking a mean spectrum over the wound was not found to yield correct wavelength-dependencies. Iterating over all possible wound spectra and selecting the best candidate was found to yield better wavelength-dependencies and lower RMSEs, but this is problematic for a larger number of pixels. Using a PCA transform to represent the possible wound spectra was found to be a viable alternative that could be fitted during optimization. This alternative has been shown to yield similar parameter RMSEs to the case where the basis spectrum is known exactly. This thus represents a suitable modification to the technique for hyperspectral images of wounds.

The layer uniqueness results show that a more direct approach technically could be taken in obtaining the reflectance and transmittance of the diffuse top layer, using a similar approach as the adding-doubling technique [59]. Such a technique would obtain reflectance and transmittance directly. A separate characterization using a one-layer model with finite thickness would be necessary for parameter estimation. The method in the current study obtains the parameters directly as a part of the procedure. Obtaining reflectance and transmittance would have to be done as a second step. Which method would be better to use would then depend on the application and the desired end result of the technique. A variant of adding-doubling would more clearly show that no unique parameters exist. It would not assume anything on the form of the optical properties of the top layer during fitting, which could be valuable as a more independent result. On the other hand, the form assumptions are necessary for enabling application of the PCA modification of the technique.

For this study, an inverse diffusion model with an offset correction obtained from Monte Carlo simulations was used. Its performance would thus be similar to an inverse Monte Carlo model with some minor inaccuracies. This allows the technique to be evaluated in terms of the basic idea rather than being overshadowed by systematic errors, while making the model suitable for hyperspectral applications. Similar corrections include a background absorption included by Svaasand et al. [44] and blood volume fraction scaling done by Randeberg et al. [46]. The offset correction is not expected to work outside the parameter range for which it was fitted, and is thus not appropriate for any unknown spectrum. More elaborate correction schemes or better model approximations are required. The model could be replaced by an empirical Monte Carlo model, or a diffusion approximation more appropriate for the absorption/scattering ratios in human tissue. Examples include the $\delta$-Eddington/$\delta$-P1 approximation [47,60]. The latter will not eliminate the offset between the model and the Monte Carlo spectra entirely [47]. In addition, it should be noted that correction factors developed for simulated reflectance spectra in an integrating sphere geometry will not directly apply to reflectance spectra from hyperspectral images. Model replacement is out of scope for the current study, where the main aim is to present and demonstrate a proof of concept. Refining the core model will be a part of future studies.

The simulations have thus verified the applicability of the technique, identified limitations and indicated what it can be used for. The technique can then be applied to experimental data. Thickness estimation of re-epithelialized areas in hyperspectral images of wounds is used as an example application.

### 3.1.5 Experimental results

In the following, hyperspectral images of an in vitro wound model sample were used to demonstrate the inverse modeling technique. The PCA-based modification in section 2.1.4 was used to represent the wound spectra using PCA during fitting. Layer thickness results over a hyperspectral image subset at days 12, 18 and 22 during the wound healing process are shown in Fig. 12. Model fits for selected spectra from day 18 are shown in Fig. 13.

The first main conclusion to be drawn from these results is that the spectral properties of the edge of the wound are explained by a gradually increasing re-epithelialization. This is modeled as a diffuse, epidermal-like layer placed on top of reflectance representing wound tissue. The layer has been shown in the simulation study to be unique. The fitted model then works as an explanatory model. The model shows that these regions have re-epithelialized, and that the optical properties here are no more than the optical properties of the wound with a typical epidermis on top. The main strength of the technique is that this can be shown without having to consider the optical properties of the dermal layers. This is a major advantage of the method as the optical properties of in vitro wound models are largely unknown. Minor changes that are challenging to identify in the RGB images can be found by the technique, as demonstrated by a thin epidermis apparently being present at the wound edge of day 12 in Fig. 12.

Depth profiles along lines placed at the approximately same position across days are shown in Fig. 14. The estimated layer thicknesses here provide a relative estimate of the re-epithelialization thickness, given that the scattering properties are homogeneous, as shown by the simulation study. In vitro wounds which are exposed to air and incubated in the medium used in this study are expected to have migration of multiple epidermal cell layers which quickly stratify into more mature epithelium [36]. This migration occurs from the edge of the wound and towards the center of the wound, with a tip of non-cornified epidermis extending on the wound side of a more mature neo-epidermis [61]. Such behavior is consistent with the migration and increases in thickness represented by the depth profiles. An epidermal thickness of magnitude 50 um is expected [3], however, indicating that higher absorption or scattering magnitudes than the ones currently used in the model are in order. Histologies were not available for the wound model samples in this study, and repetition of the experiment is necessary for proper attribution of the reflectance changes to corresponding changes in epidermal layer composition. However, the epidermal layer presence indicated by the inverse model results is in agreement with statistical characterization of these data [38].

Variations in absorption properties are not expected to be decouplable from thickness variations for real measurements. The fitted inverse model is able to match the reality in the simulations, which gives a clear minimum of the RMSE during optimization. More complex geometries or changes to the assumed optical properties broadens the minimum for real measurements, due to existence of multiple slightly sub-optimal solutions to the problem. Here, a melanin content range from 150 m$^{-1}$ to above 700 m$^{-1}$ and corresponding layer thicknesses yielded identical solutions. A clear minimum could only be found for high scattering levels, but in this case, this led to compensation by unrealistically high melanin contents.

The only absorber included in epidermis was melanin. Inclusion of a background absorption is expected to perturb the fitted parameter results, and could reduce the required melanin absorption and make the minimum mentioned above more clear. This was not investigated further, however, and tuning of such modeling details should wait until confirmation of the epidermal composition by histology. This will therefore be investigated in future work.

The simulation study indicates that at least one parameter should be fixed. All parameters were fitted simultaneously here, however, and no parameters were fixed, since the optimization seemed to produce stable estimates of both absorption and scattering properties. Only minor instabilities are evident in the day 12 profile in Fig. 14. This then shows that fitting a multi-parameter model to some reflectance might apparently give stable, unique results - but only by accident. Care must be taken, since the end result is dependent on the start parameters. Fixing at least one of the parameters is necessary for trustworthy results, as shown by the simulation study. Yet, as the parameters were stable in current case, these are the same results that would be obtained if e.g. the scattering parameters were fixed. Estimated parameters are then expected to correlate with the true parameters, as shown by the simulation study.

An infinitely wide beam illuminating a spatially invariant slab is effectively assumed in the simulations. The width of the beam is sufficiently broad with respect to the extent of the wound model sample, but the illuminated geometry is not homogeneous. Edge effects will be present in sharp transitions of tissue types, and lead to escaped photons from one type of tissue into the other. This will lead to an under- or over-estimation of thickness or absorption properties in some parts [32], but this is unlikely to have a significant effect within the slowly varying parts of the tissue. More investigation could be made in future studies into adjusting the model to the measurement geometry.

PCA was used to find a low-dimensional representation of the wound spectra that could be fitted during optimization. The PCA inverse transform with the selected number of components was found to be able to appropriately reproduce a given spectrum, and fitting the PCA scores during optimization gave reasonable results in the epidermal parameters. The simulation results indicate no significant decrease in estimation accuracy for the simulations. However, correctness for measurements should be investigated further in future studies, as the method might need tuning in e.g. number of components, or a different decomposition method than PCA might be more appropriate. The method is promising, however, and combines a data-driven, statistical approach to information extraction with physics-based photon transport modeling.

The method was tailored towards wounds, as the basis spectrum is available and can be used to fit spectra with a top layer. With adaption it might be possible to use the technique to estimate relative variations in the epidermal skin thickness of pre-term newborns and characterize burn wounds. Further, the technique is suitable for characterizing strongly absorptive inclusions in scattering media.

A spectrum from a single pixel was on average fitted in 0.66 seconds on a single CPU core, and 0.14 seconds when naively parallelizing the fitting of different pixels across 8 CPU cores (Intel Core i7-3840 QM, 2.80 GHz, 8 cores). The small subset of 50 x 40 pixels considered here would take 4 minutes and 40 seconds to fit using naive multiprocessing. The method currently runs a full, separate optimization of every pixel, which might not be needed. Future work will include adaption of GPU parallelization to reduce the running times. The current method, albeit slow, represents a proof of concept against which optimized solutions may be compared for correctness, and represents a first step towards a more scalable algorithm.

## 4. Conclusion

A technique for estimating the skin parameters of the re-epithelialized layer in wounds has been developed. The method has been found to characterize a unique diffuse layer defined by a unique reflectance and transmittance spectrum. There exists an infinite number of valid skin parameters that might characterize this layer. Fixing e.g. scattering parameters, however, can yield good relative estimates of layer thickness. The method has been used to characterize a larger area over the boundary of an in vitro wound model sample, showing the usefulness of the approach in characterizing the re-epithelialized layer. Here, a PCA modification to find the optimal wound basis spectrum has also been demonstrated, and represents a successful combination of data-driven techniques with physical photon transport modeling.

## Acknowledgments

Thanks to Ingvild Haneberg and Matija Milanic for acquisition of in vitro skin model data. Thanks to Brita Pukstad for collaboration in the in vitro wound model experiment. Thanks to Terje Schjelderup for help with the proofreading.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **T. A. Eikebrokk, B. S. Vassmyr, K. Ausen, C. Gravastrand, and B. Pukstad, “Cytotoxicity and effect on wound re-epithelialization after topical administration of tranexamic acid,” BJS Open **3**(6), 840–851 (2019). [CrossRef]

**2. **S. Lönnqvist, P. Emanuelsson, and G. Kratz, “Influence of acidic ph on keratinocyte function and re-epithelialisation of human in vitro wounds,” J. Plast. Surg. Hand Surg. **49**(6), 346–352 (2015). [CrossRef]

**3. **S. Lönnqvist, J. Rakar, K. Briheim, and G. Kratz, “Biodegradable gelatin microcarriers facilitate re-epithelialization of human cutaneous wounds - an in vitro study in human skin,” PLoS One **10**(6), e0128093 (2015). [CrossRef]

**4. **G. Kratz and C. C. Compton, “Tissue expression of transforming growth factor-β1 and transforming growth factor-αduring wound healing in human skin explants,” Wound Rep. Reg. **5**(3), 222–228 (1997).

**5. **E. Nyman, F. Huss, T. Nyman, J. Junker, and G. Kratz, “Hyaluronic acid, an important factor in the wound healing properties of amniotic fluid: In vitro studies of re-epithelialisation in human skin wounds,” J. Plast. Surg. Hand Surg. **47**(2), 89–92 (2013). [CrossRef]

**6. **G. Lu and B. Fei, “Medical hyperspectral imaging: a review,” J. Biomed. Opt. **19**(1), 010901 (2014). [CrossRef]

**7. **L. Khaodhiar, T. Dinh, K. T. Schomacker, S. V. Panasyuk, J. E. Freeman, R. Lew, T. Vo, A. A. Panasyuk, C. Lima, J. M. Giurini, T. E. Lyons, and A. Veves, “The use of medical hyperspectral technology to evaluate microcirculatory changes in diabetic foot ulcers and to predict clinical outcomes,” Diabetes Care **30**(4), 903–910 (2007). [CrossRef]

**8. **M. Denstedt, B. S. Pukstad, L. A. Paluchowski, J. E. Hernandez-Palacios, and L. L. Randeberg, “Hyperspectral imaging as a diagnostic tool for chronic skin ulcers,” Proc. SPIE **8565**, 85650N (2013). [CrossRef]

**9. **A. Nouvong, B. Hoogwerf, E. Mohler, B. Davis, A. Tajaddini, and E. Medenilla, “Evaluation of diabetic foot ulcer healing with hyperspectral imaging of oxyhemoglobin and deoxyhemoglobin,” Diabetes Care **32**(11), 2056–2061 (2009). [CrossRef]

**10. **D. Yudovsky, A. Nouvong, and L. Pilon, “Hyperspectral imaging in diabetic foot wound care,” J. Diabetes Sci. Technol. **4**(5), 1099–1113 (2010). [CrossRef]

**11. **A. Holmer, J. Marotz, P. Wahl, M. Dau, and P. W. Kämmerer, “Hyperspectral imaging in perfusion and wound diagnostic - methods and algorithms for the determination of tissue parameters,” Biomed. Tech. **63**(5), 547–556 (2018). [CrossRef]

**12. **M. A. Calin, T. Coman, S. V. Parasca, N. Bercaru, and S. R. S. D. Manea, “Hyperspectral imaging-based wound analysis using mixture-tuned matched filtering classification method,” J. Biomed. Opt. **20**(4), 046004 (2015). [CrossRef]

**13. **M. A. Calin, S. V. Parascal, D. Manea, and R. Savastru, “Hyperspectral imaging combined with machine learning classifiers for diabetic leg ulcer assessment - a case study,” in MIUA 2019: Medical Image Understanding and Analysis, vol. 1065, (2019), pp. 74–85.

**14. **G. Daeschlein, I. Langner, T. Wild, S. von Podewils, C. Sicher, T. Kiefer, and M. Jünger, “Hyperspectral imaging as a novel diagnostic tool in microcirculation of wounds,” Clin. Hemorheol. Microcirc. **67**(3-4), 467–474 (2017). [CrossRef]

**15. **L. A. Paluchowski, H. B. Nordgaard, A. Bjorgan, S. A. Berget, and L. L. Randeberg, “Can spectral-spatial image segmentation be used to discriminate burn wounds?” J. Biomed. Opt. **21**(10), 101413 (2016). [CrossRef]

**16. **M. Halicek, G. Lu, J. V. Little, X. Wang, M. Patel, C. C. Griffith, M. W. El-Deiry, A. Y. Chen, and B. Fei, “Deep convolutional neural networks for classifying head and neck cancer using hyperspectral imaging,” J. Biomed. Opt. **22**(6), 060503 (2017). [CrossRef]

**17. **H. Akbari, L. Halig, D. M. Schuster, B. Fei, A. Osunkoya, V. Master, P. Nieh, and G. Chen, “Hyperspectral imaging and quantitative analysis for prostate cancer detection,” J. Biomed. Opt. **17**(7), 0760051 (2012). [CrossRef]

**18. **J. Shapey, Y. Xie, E. Nabavi, R. Bradford, S. R. Saeed, S. Ourselin, and T. Vercauteren, “Intraoperative multispectral and hyperspectral label-free imaging: A systematic review of in vivo clinical studies,” J. Biophotonics **12**(9), e201800455 (2019). [CrossRef]

**19. **E. L. Larsen, L. L. Randeberg, E. Olstad, O. A. Haugen, A. Aksnes, and L. O. Svaasand, “Hyperspectral imaging of atherosclerotic plaques in vitro,” J. Biomed. Opt. **16**(2), 026011 (2011). [CrossRef]

**20. **A. O. H. Gerstner, W. Laffers, F. Bootz, D. L. Farkas, R. Martin, J. Bendix, and B. Thies, “Hyperspectral imaging of mucosal surfaces in patients,” J. Biophotonics **5**(3), 255–262 (2012). [CrossRef]

**21. **M. Wahabzada, M. Besser, M. Khosravani, M. T. Kuska, K. Kersting, A.-K. Mahlein, and E. Stürmer, “Monitoring wound healing in a 3d wound model by hyperspectral imaging and efficient clustering,” PLoS One **12**(12), e0186425 (2017). [CrossRef]

**22. **Y. Khouj, J. Dawson, J. Coad, and L. Vona-Davis, “Hyperspectral imaging and k-means classification for histologic evaluation of ductal carcinoma in situ,” Front. Oncol. **8**, 17 (2018). [CrossRef]

**23. **E. J. M. Baltussen, E. N. D. Kok, S. G. B. de Koning, J. Sanders, A. G. J. Aalbers, N. F. M. Kok, G. L. Beets, C. C. Flohil, S. C. Bruin, K. F. D. Kulhmann, H. J. C. M. Sterenborg, and T. J. M. Ruers, “Hyperspectral imaging for tissue classification, a way toward smart laparoscopic colorectal surgery,” J. Biomed. Opt. **24**(01), 1 (2019). [CrossRef]

**24. **B. Regeling, W. Laffers, A. O. H. Gerstner, S. Westermann, N. A. Müller, K. Schmidt, J. Bendix, and B. Thies, “Development of an image pre-processor for operational hyperspectral laryngeal cancer detection,” J. Biophotonics **9**(3), 235–245 (2016). [CrossRef]

**25. **A. Grigoroiu, J. Yoon, and S. E. Bohndiek, “Deep learning applied to hyperspectral endoscopy for online spectral classification,” Sci. Rep. **10**(1), 3947 (2020). [CrossRef]

**26. **D. Yudovsky, A. Nouvong, K. Schomacker, and L. Pilon, “Monitoring temporal development and healing of diabetic foot ulceration using hyperspectral imaging,” J. Biophotonics **4**(7-8), 565–576 (2011). [CrossRef]

**27. **S. Vyas, A. Banerjee, and P. Burlina, “Estimating physiological skin parameters from hyperspectral signatures,” J. Biomed. Opt. **18**(5), 057008 (2013). [CrossRef]

**28. **W. Feng, R. Shi, C. Zhang, T. Yu, and D. Zhu, “Lookup-table-based inverse model for mapping oxygen concentration of cutaneous microvessels using hyperspectral imaging,” Opt. Express **25**(4), 3481–3495 (2017). [CrossRef]

**29. **E. Zherebtsov, V. Dremin, A. Popov, A. Doronin, D. Kurakina, M. Kirillin, I. Meglinski, and A. Bykov, “Hyperspectral imaging of human skin aided by artificial neural networks,” Biomed. Opt. Express **10**(7), 3545–3559 (2019). [CrossRef]

**30. **M. Kewin, A. Rajaram, D. Milej, A. Abdalmalak, L. Morrison, M. Diop, and K. S. Lawrence, “Evaluation of hyperspectral nirs for quantitative measurements of tissue oxygen saturation by comparison to time-resolved nirs,” Biomed. Opt. Express **10**(9), 4789–4802 (2019). [CrossRef]

**31. **D. Yudovsky and L. Pilon, “Rapid and accurate estimation of blood saturation, melanin content, and epidermis thickness from spectral diffuse reflectance,” Appl. Opt. **49**(10), 1707–1719 (2010). [CrossRef]

**32. **A. Bjorgan, M. Milanic, and L. L. Randeberg, “Estimation of skin optical parameters for real-time hyperspectral imaging applications,” J. Biomed. Opt. **19**(6), 066003 (2014). [CrossRef]

**33. **M. S. Patterson, B. C. Wilson, and D. R. Wyman, “The propagation of optical radiation in tissue. ii: Optical properties of tissues and resulting fluence distributions,” Lasers Med. Sci. **6**(4), 379–390 (1991). [CrossRef]

**34. **G. Zonios and A. Dimou, “Modeling diffuse reflectance from semi-infinite turbid media: application to the study of skin optical properties,” Opt. Express **14**(19), 8661–8674 (2006). [CrossRef]

**35. **A. Kim and B. C. Wilson, * Optical-Thermal Response of Laser-Irradiated Tissue* (Springer, 2011), chap. Measurement of ex vivo and in vivo tissue optical properties: Methods and theories, pp. 267–319, 2nd ed.

**36. **K. Jansson, G. Kratz, and A. Haegerstrand, “Characterization of a new in vitro model for studies of reepithelialization in human partial thickness wounds,” In Vitro Cell. Dev. Biol.: Anim. **32**(9), 534–540 (1996). [CrossRef]

**37. **G. Kratz, “Modeling of wound healing processes in human skin using tissue culture,” Microsc. Res. Tech. **42**(5), 345–350 (1998). [CrossRef]

**38. **A. Bjorgan, B. Pukstad, and L. L. Randeberg, “Hyperspectral characterization of re-epithelialization in an in vitro wound model,” J. Biophotonics (2020).

**39. **G. James, D. Witten, T. Hastie, and R. Tibshirani, * An introduction to statistical learning* (Springer, 2013), 1st ed.

**40. **G. Shaw and D. Manolakis, “Signal processing for hyperspectral image exploitation,” IEEE Signal Process. Mag. **19**(1), 12–16 (2002). [CrossRef]

**41. **J. Khodr and R. Younes, “Dimensionality reduction on hyperspectral images: A comparative review based on artificial datas,” in * 2011 4th International Congress on Image and Signal Processing* (IEEE, 2011), pp. 1875–1883.

**42. **M. D. Farrell and R. M. Mersereau, “On the impact of pca dimension reduction for hyperspectral detection of difficult targets,” IEEE Geosci. Remote Sens. Lett. **2**(2), 192–195 (2005). [CrossRef]

**43. **L. Wang, S. L. Jacques, and L. Zheng, “Mcml monte carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. **47**(2), 131–146 (1995). [CrossRef]

**44. **L. Svaasand, L. Norvang, E. Fiskerstrand, E. Stopps, M. Berns, and J. Nelson, “Tissue parameters determining the visual appearance of normal skin and port-wine stains,” Laser. Med. Sci. **10**(1), 55–65 (1995). [CrossRef]

**45. **L. V. Wang and H. Wu, * Biomedical Optics, Principles and Imaging* (John Wiley & Sons, 2007).

**46. **L. L. Randeberg, A. Winnem, R. Haaverstad, O. A. Haugen, and L. O. Svaasand, “Performance of diffusion theory vs. monte carlo methods,” Proc. SPIE **5862**, 586200 (2005). [CrossRef]

**47. **T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. **39**(34), 6453–6465 (2000). [CrossRef]

**48. **R. R. Anderson and J. A. Parrish, “The optics of human skin,” J. Invest. Dermatol. **77**(1), 13–19 (1981). [CrossRef]

**49. **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]

**50. **T. Lister, P. A. Wright, and P. H. Chappell, “Optical properties of human skin,” J. Biomed. Opt. **17**(9), 0909011 (2012). [CrossRef]

**51. **N. Verdel, A. Marin, M. Milanic, 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–960 (2019). [CrossRef]

**52. **L. Vidovic, M. Milanic, and B. Majaron, “Objective characterization of bruise evolution using photothermal depth profiling and monte carlo modeling,” J. Biomed. Opt. **20**(1), 017001 (2015). [CrossRef]

**53. **E. Alerstam, W. C. Y. Lo, T. D. Han, J. Rose, S. Andersson-Engels, and L. Lilge, “Next-generation acceleration and code optimization for light transport in turbid media using gpus,” Biomed. Opt. Express **1**(2), 658–675 (2010). [CrossRef]

**54. **E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed monte carlo simulation of photon migration,” J. Biomed. Opt. **13**(6), 060504 (2008). [CrossRef]

**55. **T. Spott, L. O. Svaasand, R. E. Anderson, and P. F. Schmedling, “Application of optical diffusion theory to transcutaneous bilirubinometry,” Proc. SPIE **3195**, 234–245 (1998). [CrossRef]

**56. **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]

**57. **S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. **58**(11), R37–R61 (2013). [CrossRef]

**58. **P. Naglic, L. Vidovic, M. Milanic, L. L. Randeberg, and B. Majaron, “Suitability of diffusion approximation for an inverse analysis of diffuse reflectance spectra from human skin in vivo,” OSA Continuum **2**(3), 905–922 (2019). [CrossRef]

**59. **S. A. Prahl, *Optical-Thermal Response of Laser Irradiated Tissue* (1995), chap. 5, The adding-doubling method, pp. 101–129.

**60. **S. A. Carp, S. A. Prahl, and V. Venugopalan, “Radiative transport in the delta-p1 approximation: accuracy of fluence rate and optical penetration depth predictions in turbid semi-infinite media,” J. Biomed. Opt. **9**(3), 632–647 (2004). [CrossRef]

**61. **G. D. Glinos, S. H. Verne, A. S. Aldahan, L. Liang, K. Nouri, S. Elliot, M. Glassberg, D. C. DeBuc, T. Koru-Sengul, M. Tomic-Canic, and I. Pastar, “Optical coherence tomography for assessment of epithelialization in a human ex vivo wound model,” Wound Rep. Reg. **25**(6), 1017–1026 (2017). [CrossRef]