Exploiting scale-invariance: a top layer targeted inverse model for hyperspectral images of wounds

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


Introduction
In vitro wound models are useful for investigation of wound healing in a controlled laboratory setting [1][2][3][4][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][7][8][9][10][11][12][13][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][20][21][22][23][24][25]. The rich spectral content further enables use of inverse photon transport modeling [26][27][28][29][30][31][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][41][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.

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 µ a /µ s (scale-invariance).

2.
A multi-layer model and a one-layer representation yield identical reflectance values when a top layer is added.
The method considers cases where reflectance spectra are available from two regions, with and without a top layer. Given the first assumption, the reflectance from the region without the top layer is represented using a one-layer model by estimating its corresponding µ a /µ s ratio without having to consider their actual forms. From the second assumption, using these properties in the deeper layer of a two-layer model can be used to represent the reflectance from the region with the top layer. The skin properties of the top layer can then be fitted without considering the deeper layers if the optical properties of the top layer are known. Insertion of one-layer properties from the region without the top layer into the two-layer model ensures that the boundary conditions towards the top layer are correct. The basic model steps are shown in Fig. 1. A one-layer model is fitted to the reflectance R 0 (λ) from a region with the top layer missing. Due to scale-invariance, any absorption coefficient µ a and scattering coefficient µ s obeying the required ratio are viable solutions. A top layer can then be fitted to the model by reusing the optical properties in the deeper layer of a two-layer model.

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 µ a /µ s directly from the reflectance without iteration.
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 = 1 3(µ s +µ a ) and φ 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] The boundary condition at the air-tissue interface is taken as j(z = 0) = Aφ(z = 0). The property j is the photon flux. The boundary condition essentially relates the irradiance propagating back into the tissue to the irradiance propagating out of the tissue by an effective reflection coefficient [44]. A refraction index of n = 1.4 yields A = 0.17 [44]. It is required that lim z→∞ φ(z) = 0. All constants A ij can then be determined by using continuity in j(z) and φ(z) between each layer and the boundary conditions above [44]. The diffuse reflectance R d is found by [44] The last expression is obtained by considering the irradiance transmitted into the air. Analytic solution for a two-layer model can be found in Svaasand et al. [44]. The one-layer solution is trivial to obtain.
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.

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. 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). [53,54] was used to simulate reflectance spectra from a skin-like geometry. NVIDIA GeForce GTX 670 was used for GPU parallelization. Wavelengths from λ = 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.

Absorption properties
The absorption in the top layer was modeled using an absorption model for melanin [55] µ a,e = µ a,m,694 (λ/694) −3.46 , where µ 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 µ a,d = µ oxy (λ)c oxy + µ deoxy (λ)c deoxy , where µ {oxy, deoxy} (λ) are the absorption spectra for oxygenated and deoxygenated blood, respectively. Assumed low and high values are listed in Table 2. These cover blood volume fractions from 2% to 10%, and oxygenations from approx. 20% to 80%. The inverse model avoids consideration of this layer, and the main goal is to model a layer with absorption and scattering magnitudes representing human skin.
The parameters µ s,Mie,500 and b Mie describe Mie scattering. Examples of values for ex vivo human skin are 1800 m −1 and 0.22, respectively [56]. The parameter µ s,Ray,500 describes Rayleigh scattering. Example of an ex vivo value is 1700 m −1 [56]. Low and high values are listed in Table 2, and are assumed to represent the expected magnitudes in human skin. The parameter b Mie is kept constant at 0.22 [56] across all simulations. This parameter is expected to vary [57], however, changes in the coefficients were thought to represent a similar change in the wavelength-dependency that could test recovery of b Mie . Further, the reflectance was found to be less sensitive to changes in scattering in the epidermal layer, and full parameter recovery of b Mie was not expected. For simplicity, the parameter was not varied in the simulations.
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.

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 µ a (65 to 6714 m −1 ) and µ s (1517 to 5237 m −1 ). The assumption that the reflectance can be written as a function of µ a /µ s is checked by plotting the output reflectance values as a function of µ a /µ 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 µ a /µ 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 µ a /µ s ratio for a one-layer model from the complex multi-layer reflectance. A µ s was assumed (µ s = 1000m −1 ), and a µ 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.

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 µ a /µ 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.

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.
3. Fit only thickness and melanin content, with scattering parameters fixed.
Model C2 picks epidermal and dermal parameters randomly across a continuous range, and is suitable for evaluating the parameter resolving performance, evaluating at multiple basis spectra and evaluating the case where the basis spectrum is not known. The same epidermal scattering parameters are used, representing a case where the epidermal scattering is known to be homogeneous. Here, three new cases of the inverse model were run: 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.
These tests should then elucidate the accuracy of the inverse model under various conditions, from which conclusions may be drawn about its application to real measurement data.

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.

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.

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 µ a /µ 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 µ a /µ s . The latter clearly demonstrates that there is a one-to-one correspondence between the µ a /µ 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.

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 µ a , µ 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 µ a for several µ 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  Fig. 5. Demonstration using Monte Carlo simulations that a multi-layer model and its one-layer fit have identical reflectance spectra when adding additional top layers: Comparison of the reflectance spectrum from a multi-layered model and the reflectance spectra from a one-layer model constructed from each layer (top), and the reflectance from the multi-layered model with an epidermis on top compared to a single-layer approximation with the same epidermis on top (bottom). The consequence is that a one-layer model is appropriate for representing a wound spectrum when evaluating the effect of adding an epidermis to the wound spectrum. The simulations were run with 1000 000 photons per wavelength. Fig. 6. Demonstration of the coupling between d, µ a and µ s in the top layer in the two-layer diffusion model, for a single reflectance value. All these parameter combinations produce the same output reflectance (R = 0.6). 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. Fig. 7. Results from fitted inverse model at random start parameters, demonstrating the coupling of the fitted parameters. All parameter combinations produce the same reflectance spectra. Top left: Estimated layer thickness versus estimated melanin concentration. Top right: Estimated layer thickness versus scattering parameters. Bottom: Estimated layer thickness versus fit RMSE. Fitted parameters at random start parameters for thickness and melanin, with the scattering parameters fixed to the true scattering parameters, are marked as "Fixed scattering" in the plots.
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 µ a , µ s and d 1 . The absorption coefficient µ 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 µ s . The scattering model can be re-written to µ s = A f (λ/500 where A = µ s,Mie,500 + µ s,Ray,500 and f = µ 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 µ s,Mie,500 and µ 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.

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 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%).

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. Fig. 11. Modeled parameters versus estimated parameters across Monte Carlo simulations with random input parameters, fixed scattering. Two cases are shown: Scattering parameters fixed to the true scattering (blue), and scattering parameters fixed to the parameters estimated from a single reflectance spectrum (green). The former leads to reasonable estimates of thickness and melanin content, while the latter leads to estimates that are only correlated with the true parameters.
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.

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 µ a /µ s . This means that any combination of µ a and µ 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(λ) and T(λ) 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 δ-Eddington/δ-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.

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]. the RGB image with a dotted square indicating the image subset considered in the inverse model, the RGB values of the reflectance from the one-layer dermis reconstruction within the subset, and the estimated thickness of the epidermis within the subset. RGB images were constructed from the hyperspectral images at 615, 564 and 459 nm wavelength bands, and were gamma-corrected for increased contrast. The white-pink region corresponds to wound, while the brown region is intact tissue. Coordinates of the spectra plotted in Fig. 13 are marked in the day 18 image.  Fig. 12. The peaks between 690 and 750 nm are artifacts due to mismatch of the order sorting filter in the hyperspectral camera, and were not fitted by the model. Fitted parameters were the melanin content, layer thickness and scattering parameters in epidermis, and the three PCA coefficients for the wound basis in dermis. 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.

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.