1.

Introduction

Hyperspectral imaging has been recently adopted for diagnostic imaging of human skin. Hyperspectral imaging combines high spatial and spectral resolution in one modality, giving images with full spectral resolution in every pixel.^1,2 This makes it a promising tool for tissue characterization and optical diagnostics.^1,3,4 Flexible wide-field imaging options give the possibility of rapid scanning of both larger areas and samples and smaller details, e.g., full body scans and close up imaging with microscopic resolution.⁵ The physical size of the equipment makes it possible to place it on a small trolley, and it can thus be used for bedside scanning.⁶ The hyperspectral camera used in this study is a push broom, line-scanning device with a capture speed of 30 ms per line of data.

The amount of collected data is an obstacle for fast hyperspectral data processing while scanning. With typical files being in the order of several gigabytes, processing time is prohibitive for advanced analysis. Time is a crucial parameter for the technique to be clinically relevant. Processing speed should preferably obey an external real-time deadline limit defined by the acqusition speed of the hyperspectral camera. Commodity graphics processing units (GPUs) are well known to satisfy heavy computing requirements.^7,8 GPUs are normally used to process graphics on personal computers and are relatively inexpensive. Such graphics cards can potentially boost processing speed due to the inherent parallelizability of hyperspectral data processing.

It is valuable to estimate the concentrations of various tissue components (blood, melanin, and water) for diagnostic purposes. Spectral unmixing algorithms are frequently used to obtain concentration maps of materials in remote sensing.⁹ Applying the same methods to analyze the hyperspectral images of skin is challenging due to the turbid, layered nature of human skin. One aim of this study is to develop a method for spectral unmixing of the hyperspectral images of human skin.

Forward light transport models, such as Monte Carlo¹⁰ or the diffusion approximation,^11–13 may be used to simulate light transport in media such as human skin. These type of models have previously been used to extract optical properties of tissue.^13–15 This has been typically done by comparing measured, diffuse reflectance spectra to simulations, and tweaking the input parameters of the simulations to obtain the best possible fit. To our knowledge, these methods have never been applied to hyperspectral images in real-time during data collection. Some reports exist on the application of inverse photon transport models on hyperspectral images after image acquisition.^16,17

This paper presents a deterministic hyperspectral inverse modeling approach based on diffusion theory and spectral unmixing. These analytical methods are suitable for GPU implementation. The analysis chain fulfills the real-time requirements of the hyperspectral imaging process. The optical properties are delivered line by line during data collection, and no additional computational time is needed beyond the time it takes to scan the sample.

The accuracy of the results obtained using the inverse model was evaluated against numerical simulations. However, as this is a proof of a new concept, real-time processing concerns were given priority over simulation accuracy. This issue will be addressed in further development of the concept.

2.

Materials and Methods

2.2.

Optical Modelling

2.2.1.

Inverse model

Skin model. The skin model applied in the inverse model is a two-layered model which was previously described by Spott et al.,¹³ Randeberg et al.,^20,21 and Svaasand et al.¹¹ In short, it consists of two planar layers. The first layer represents an epidermal layer of finite thickness containing melanin and a small amount of blood. Blood is added to the epidermal layer to correct for the depth variations of the papillae. The second layer is a homogeneous, semi-infinite dermal layer containing blood and other chromophores evenly distributed throughout the layer.

The applied melanin absorption model was adapted from Spott et al.¹³ and is given as

Eq. (1)

$${\mu }_{a,m}(\lambda )={\mu }_{a,m,694}·{\left(\frac{\lambda }{694\mathrm{nm}}\right)}^{-3.46}.$$

The melanin content is denoted by the melanin absorption at 694 nm (${\mu }_{a,m,694}$) in units of ${\mathrm{m}}^{-1}$. The absorption in blood is assumed to arise from deoxygenated and oxygenated hemoglobins.²² It depends on the blood volume fraction (BVF) and oxygen saturation (oxy). A constant background absorption of $25{\mathrm{m}}^{-1}$ was added to the model.¹¹ Absorption spectra of hemoglobin, fat, and melanin are shown in Fig. 1.

Fig. 1

The absorption coefficients of blood,²² water,²³ fat,²⁴ and the melanin model¹³ for sun-protected north European skin (${\mu }_{a,m,694}=350{\mathrm{m}}^{-1}$).

In the model, the absorption spectra are multiplied by their volume fractions and summed to yield the total wavelength-dependent dermal absorption coefficient, ${\mu }_{a,d}(\lambda )$.

The reduced scattering coefficient in the tissue is modeled by an expression given by Jacques²⁵

Eq. (2)

$${\mu }_s^{\prime }(\lambda )=100{\mathrm{m}}^{-1}·[a_{\mathrm{Mie}}{\left(\frac{\lambda }{500\mathrm{nm}}\right)}^{-b_{\mathrm{Mie}}}+a_{\mathrm{Ray}}{\left(\frac{\lambda }{500\mathrm{nm}}\right)}^{-4}].$$

Bashkatov et al.²⁶ give $a_{\mathrm{Mie}}=18.780$, $a_{\mathrm{Ray}}=17.6$, and $b_{\mathrm{Mie}}=0.22$ as values for the coefficients in the expression. This is based on data obtained from whole skin ex vivo samples, and is used for both dermis and epidermis. The anisotropy factor $g$ is modeled as

Eq. (3)

$$g(\lambda )=0.62+\lambda ·29·{10}^{-5},$$

reported by Van Gemert et al.²⁷ The index of refraction used for both skin layers was 1.4.

Light transport model

The light transport model is based on a diffusion model with isotropic source functions, derived by Svaasand et al.¹¹ This model is used to obtain diffuse reflectance as a function of wavelength and skin optical parameters. The diffusion model has fast, analytic solutions well suitable for a real-time computing environment. In the diffusion approximation, it is assumed that scattering dominates over absorption. This assumption is less valid below 600 nm for skin due to the high absorption coefficients for hemoglobin and melanin (see Fig. 1). It has been shown that the approximation still results in reflectances close to Monte Carlo simulations for the 450–800 nm spectral range.²⁰

The error between the Monte Carlo model and the diffusion approximation was reported by Randeberg et al.²⁰ to be minimized by applying a constant scaling factor to the absorption coefficients in the diffusion model. This scaling factor has not been explicitly applied. However, a constant deviation is present in all extracted parameters, as found by Randeberg et al.²⁰ The extracted parameters are thus expected to characterize the relative variations of properties between different tissues well.

A Monte Carlo-based inverse model will not fulfill real-time requirements, even if they are GPU-accelerated.⁷ Typical computing times for 1600 pixels and 160 wavelengths, based on our own experiences, are in the order of days.

The derivation of the diffusion model assumes a laterally homogeneous skin model within the range of a few mean transport lengths. A valid application of the diffusion model to single hyperspectral pixels would require the lateral broadening of the light to be smaller than the spatial width of the hyperspectral pixel. Since this is not the case, photons will be scattered from pixel to pixel. The properties extracted from a single pixel will therefore be influenced by the surrounding pixels. This effect has been further investigated in this study.

Fitting algorithm

The diffusion model is applied to calculate the wavelength-dependent diffuse reflectance. By comparing the calculated reflectance against a measured reflectance, the skin parameters may be found iteratively. The iteration strategy is designed with GPU implementation in mind, requiring independence across wavelengths and pixels and determinism in the set of instructions used. The iteration strategy is therefore based around the independent estimation of the epidermal and dermal absorption coefficients using Newton–Rhapson’s method²⁸ and the spectral unmixing of these.

The number of layers is kept down to two layers in order to keep the number of fitted parameters down to a minimum, thus simplifying the reflectance expression. It is known that the skin is a nonhomogeneous organ and is approximated more accurately by three or more skin layers. The two-layered skin model is, therefore, applied on separate wavelength ranges where the penetration depth can be assumed to be more or less uniform across the given wavelength range. It is assumed that this approach will yield the mean properties down to the given penetration depth. This is illustrated in Fig. 2 and has been done previously by Randeberg et al.²¹ This approach is also similar to work done by Tseng et al.²⁹ and Saager et al.³⁰ The two chosen fitting intervals are 510–590 and 690–820 nm. Light penetrates more superficially at the shorter wavelengths due to the high absorption of hemoglobin and melanin, while light penetrates more deeply at the longer wavelengths (see Fig. 1). The exact limits have been chosen through thorough testing on real data and close observation of the chromophore absorption spectra. Typical mean optical penetration depths through dermis are estimated to around 200 μm for 510–590 nm and above 500 μm for 690–820 nm.

Fig. 2

Illustration of the two-layered skin model applied to a three-layered situation. The two-layered skin model approximates the properties of multiple, inhomogeneous layers to a single, homogeneous and semi-infinite (s.i.) layer where the derived properties are distributed evenly throughout the layer. Different parts of the diffuse reflectance spectrum contain information about different penetration depths.

The fixed parameters within the chosen skin model are:

• • scattering coefficients [Eq. (2)]

• • thickness of the layers (100 μm for epidermis).

The parameters used in Eq. (2) are assumed to have a slow spatial variation for normal skin tissues and can be assumed fixed, although they may be different for different types of tissues. The thickness of the epidermis varies by person and location, and thickness under- or overestimation results in either a lowered or increased melanin content.

The fitted parameters within the chosen skin model, as applied to a given wavelength interval, are:

• • melanin absorption in epidermis

• • BVF in dermis

• • oxygen saturation of the blood in dermis

• • other chromophores (in this paper: water, fat, and constant baseline absorption).

The basic fitting procedure is outlined in Fig. 3. The dermal absorption coefficient ${\mu }_{a,d}$ is estimated and unmixed at two wavelength intervals following determination of melanin in epidermis. The melanin extraction method is outlined in Fig. 4. The melanin absorption coefficient ${\mu }_{a,m,694}$ is assumed to be $100{\mathrm{m}}^{-1}$ as a starting value.

Fig. 3

Outline of the inverse model algorithm, where the absorption in dermis is corrected for melanin in epidermis.

Fig. 4

Outline of the melanin estimation algorithm. The method is run on the wavelength interval from 730 to 830 nm, and is run twice.

The dermal absorption coefficients are derived for the wavelength interval 730–830 nm. Melanin and hemoglobin are then fitted to the absorption coefficients. Next, the hemoglobin parameters are fixed and the epidermal absorption coefficients are derived. Melanin is then fitted to the epidermal absorption coefficients to get the corresponding melanin content. Using the new melanin estimate, the method is run a second time to get the final estimate of the melanin content.

The wavelength-dependent dermal absorption coefficient is then found across the entire spectral range. Due to a linear relationship between the chromophore absorption coefficients and the total dermal absorption, the chromophore contributions to the dermal absorption spectrum are found using a spectral unmixing algorithm.

The chromophore spectra for this unmixing procedure are known a priori and only estimation of chromophore volume fractions is desired for the spectral unmixing part. For a single pixel, the derived dermal absorption can be written as

Eq. (4)

$${\overrightarrow{\mu }}_{a,d}=A\overrightarrow{x},$$

where the matrix $A$ consists of the chromophore absorption spectra and $\overrightarrow{x}$ are the concentrations $\{x_0,x_1,x_2,\dots \}$ of each material to be found under a non-negativity constraint.

The sequential coordinate-wise algorithm for non-negative least-squares problems (SCA)³¹ is used to estimate the solutions to Eq. (4), and thus unmix the absorption spectra. This algorithm has no proven convergence guarantees, but is reportedly fast. It requires only the matrix $H=A^{T}A$ and optimizes with respect to one variable at a time, suitable for a memory-effective GPU parallelization in a demonstration prototype.

The unmixing of the dermal absorption coefficient determines the fitted chromophores (hemoglobin, water, fat, and the constant baseline absorption) in two separate wavelength intervals. In addition, melanin is included in the fitting to rectify potential melanin underestimation in epidermis.

2.2.2.

Numerical simulations

Reflectance spectra simulated using a four-layered Monte Carlo model were used to test the reliability of the inverse diffusion model. In the Monte Carlo model, the dermis was subdivided into two layers, the first representing papillary dermis, and the second reticular dermis. The fourth layer represented a blood-less subcutaneous layer containing 40% fat and 60% background absorption. The implemented background absorption was taken from Salomatina et al.³²

Eq. (5)

$${\mu }_{a,b}(\lambda )=100{\mathrm{m}}^{-1}·(0.82+16.82e^{-(\lambda -400\mathrm{nm})/80.5\mathrm{nm}}).$$

The reduced scattering coefficient in subcutis was modeled by an expression given by Naglic et al.,³³ comparable to the scattering coefficients reported by Salomatina et al.,³²

Eq. (6)

$${\mu }_s^{\prime }(\lambda )=1500{\mathrm{m}}^{-1}·(16.34+303.8e^{-\lambda /180.3\mathrm{nm}}).$$

The refraction index $n$ was set to 1.4 for all layers. Dermis was modeled using blood and the constant background absorption presented earlier in Sec. 2.2.1. The reduced scattering coefficient for dermis was calculated using Eq. (2). This will serve to ease the applicability of the inverse diffusion model since the same optical properties are modeled. The main purpose of the simulations is to:

• • Test the accuracy of the optical parameters given by the inverse diffusion model.

• • Test the applicability of the one-dimensional (1-D) diffusion model to describe a 3-D situation.

A 1-D Monte Carlo model⁷ was used to simulate reflectance spectra. The melanin in epidermis and the BVF in the superficial dermal layer were varied. The inverse model was tested for its ability to estimate the changes in these parameters and the stability of the other estimated parameters (oxygenation and BVF in the deeper dermal layer). The skin model is shown in Fig. 5. The melanin absorptions corresponded to lightly and more pigmented Caucasian skin. Diffuse reflectance spectra in the wavelength range between 400 and 848 nm with a step size of 2 nm were simulated.

Fig. 5

Geometry of the one-dimensional (1-D) Monte Carlo forward model used to simulate skin reflectance. Parameter values within brackets are varied, the other parameters are fixed.

The inverse diffusion model used to derive optical properties from the 1-D Monte Carlo spectra involved the use of an ordinary non-negative least-squares algorithm³⁴ for the unmixing of the absorption spectra, instead of SCA unmixing as presented earlier. This was done in order to evaluate the performance of SCA against a more ordinary non-negative least-squares algorithm.

The iteration strategy has been evaluated against results obtained using MATLAB’s (Version 8.1.0.604, The MathWorks Inc., Natick, Massachusetts, USA) lsqcurvefit routine, using a Trust-region reflective optimization approach with upper and lower boundary constraints. This was run on the wavelength range 690–820 nm for all 1-D Monte Carlo spectra.

We also performed 3-D Monte Carlo simulations to test the applicability of the 1-D diffusion model on single pixels in a mock hyperspectral image. An implementation of 3-D Monte Carlo developed by Milanic and Majaron³⁵ was used to obtain reflectance images of a 1-mm diameter junctional nevus, and an intradermal vessel at a depth of 0.25 mm and with a diameter of 0.2 mm. Skin geometries are shown in Figs. 6 and 7. The simulated wavelength range for the diffuse reflectance was 42 wavelengths from 687 to 835 nm with a step size of 3.6 nm. This wavelength discretization was used in order to simulate a hyperspectral dataset.

Fig. 6

Geometry of the skin model containing a superficial vessel. The vessel was modeled as a circular tube. A horizontal cross section through the center of the vessel is shown on top. A vertical cross section is shown below. The image scale of the vertical cross section differs from the horizontal cross section.

Fig. 7

Geometry of the skin model containing a mole. The mole was modeled as a circular area with increased melanin content. A horizontal cross section through the center of the mole is shown on top. A vertical cross section is shown below. The image scale of the vertical cross section differs from the horizontal cross section.

The blood vessel was simulated using the blood absorption and scattering values from Friebel et al.³⁶ and a blood oxygenation of 95%. The refraction index of the vessel was set to

Eq. (7)

$$n(\lambda )=1.357+\frac{6.9·{10}^3}{{\lambda }^2}+\frac{7.6·{10}^8}{{\lambda }^4},$$

as reported by Li et al.³⁷

The inverse diffusion model used for the 3-D Monte Carlo signals involved SCA unmixing for the unmixing of the absorption spectra, as presented earlier in Sec. 2.2.1.

3.

Results

3.1.

Timing Results

The specific timing results for each GPU operation as applied on a 1600 samples $\times$ 160 bands hyperspectral data line are shown in Fig. 8. The final results, i.e., melanin and the blood parameters from different wavelength intervals, are delivered within 3.5 ms. This computation time is well within the real-time deadline limit imposed by the hyperspectral system, meeting the fast computing requirement, and leaving GPU time for other future processing operations.

Fig. 8

Total computational times for the inverse modelling of one hyperspectral line of data. Three wavelength intervals were used in the unmixing of the dermal absorption, the line of data had $1600\text{samples}\times 160\text{bands}$ of data.

3.2.

Simulation Results

3.2.1.

One-dimensional Monte Carlo modeling

The GPU inverse diffusion model was evaluated on diffuse reflectance spectra obtained by a multilayered, 1-D Monte Carlo model. This was compared against an ordinary inverse diffusion model based on multivariate objective fitting (due to potential confusion, henceforth referred to as “the objective fitting”). The BVFs extracted from the wavelength range 510–590 nm are plotted in Fig. 10. The means and standard deviations of the other parameters are displayed in Table 1. A full list of extracted parameters is shown in Table 2. Some Monte Carlo spectra and their inverse diffusion model fits are shown in Fig. 9.

Table 1

Mean and standard deviation for extracted parameters from the multilayered Monte Carlo forward model. Statistics were computed across the dataset generated by varying the blood volume fraction (BVF) in the superficial dermal layer from 0.5% to 3.0% in steps of 0.5%. The extracted parameters are the melanin absorption at 694 nm (μa,m,694), the oxygenation in the 510–590 and 690–820 nm wavelength intervals (oxy510−590, oxy690−820) and the BVF in the 690–820 nm wavelength interval (BVF690−820). Inv. DM is the inverse diffusion model, obj. fit is the inverse diffusion model using a multivariate objective fitting scheme. Input oxygenation for the superficial and deep dermal layer was 50% and 80%, respectively. BVF in the deep dermal layer was 3%.

Method	Input ${\mu }_{a,m,694}$ (${\mathrm{m}}^{-1}$)	Output ${\mu }_{a,m,694}$ (${\mathrm{m}}^{-1}$)	${\mathrm{oxy}}_{510-590}$	${\mathrm{oxy}}_{690-820}$	${\mathrm{BVF}}_{690-820}$
Inv. DM	227	$264\pm 7$	$0.65\pm 0.01$	$0.77\pm 0.01$	$0.026\pm 0.002$
Inv. DM	683	$692\pm 11$	$0.67\pm 0.01$	$0.76\pm 0.03$	$0.026\pm 0.002$
Obj. Fit	227	$230\pm 39$	$0.66\pm 0.02$	$0.76\pm 0.02$	$0.025\pm 0.002$
Obj. Fit	683	$682\pm 31$	$0.70\pm 0.03$	$0.77\pm 0.02$	$0.028\pm 0.002$

Table 2

Extracted parameters from the one-dimensional (1-D) Monte Carlo forward model for the different inverse models. The number below the parameter name denotes from which wavelength range the parameter was extracted (510–590 or 690–820 nm). μa,m is the epidermal melanin absorption coefficient at 694 nm. Input oxygenation for the superficial and deep dermal layer was 50% and 80%, respectively. Blood volume fraction (BVF) in the deep dermal layer was 3%.

Input parameters		Inv. DM					Objective fit, DM
${\mu }_{a,m}$	${\mathrm{BVF}}_1$	${\mu }_{a,m}$	${\mathrm{oxy}}_1$	${\mathrm{BVF}}_1$	${\mathrm{oxy}}_2$	${\mathrm{BVF}}_2$	${\mu }_{a,m,1}$	${\mathrm{oxy}}_1$	${\mathrm{BVF}}_1$	${\mu }_{a,m,2}$	${\mathrm{oxy}}_2$	${\mathrm{BVF}}_2$
227	0.005	254	0.66	0.017	0.77	0.023	248	0.63	0.018	260	0.78	0.023
227	0.0075	263	0.65	0.019	0.79	0.024	235	0.63	0.020	261	0.78	0.024
227	0.01	260	0.64	0.022	0.78	0.025	226	0.63	0.022	259	0.78	0.024
227	0.0125	266	0.63	0.024	0.79	0.024	218	0.64	0.025	265	0.78	0.024
227	0.015	252	0.64	0.027	0.76	0.025	206	0.66	0.027	259	0.76	0.025
227	0.0175	271	0.65	0.028	0.79	0.025	201	0.67	0.030	266	0.78	0.025
227	0.02	262	0.64	0.031	0.76	0.027	187	0.66	0.033	260	0.76	0.026
227	0.0225	264	0.65	0.033	0.75	0.027	179	0.68	0.036	257	0.74	0.027
227	0.025	273	0.65	0.035	0.76	0.028	172	0.68	0.039	261	0.75	0.027
227	0.0275	267	0.66	0.038	0.75	0.029	160	0.69	0.042	260	0.74	0.027
227	0.03	275	0.66	0.040	0.76	0.028	151	0.70	0.045	265	0.75	0.027
683	0.005	697	0.70	0.019	0.80	0.022	702	0.67	0.019	710	0.81	0.025
683	0.0075	710	0.67	0.021	0.80	0.026	690	0.67	0.022	706	0.79	0.026
683	0.01	709	0.67	0.024	0.79	0.026	686	0.68	0.024	706	0.79	0.026
683	0.0125	689	0.67	0.028	0.76	0.025	673	0.69	0.027	703	0.78	0.027
683	0.015	701	0.64	0.030	0.78	0.026	669	0.68	0.030	705	0.78	0.027
683	0.0175	685	0.65	0.034	0.78	0.027	659	0.70	0.033	711	0.78	0.027
683	0.02	679	0.67	0.037	0.73	0.023	648	0.72	0.036	705	0.76	0.028
683	0.0225	680	0.66	0.040	0.74	0.027	639	0.72	0.039	702	0.76	0.029
683	0.025	687	0.67	0.042	0.73	0.026	633	0.73	0.042	705	0.75	0.029
683	0.0275	687	0.67	0.045	0.73	0.027	622	0.73	0.046	704	0.75	0.030
683	0.03	693	0.68	0.048	0.74	0.028	621	0.75	0.049	701	0.74	0.030

Fig. 9

Monte Carlo simulations and their corresponding diffusion model fits. Spectrum 1 corresponds to a melanin absorption of $227{\mathrm{m}}^{-1}$ while both spectra 2 and 3 have a melanin absorption of $683{\mathrm{m}}^{-1}$. Blood volume fractions (BVFs) in the upper dermal layer are respectively 1%, 1%, and 3%. Dots in the graph are not indicative of fitting points and are only used for discrimination. Fitted wavelength ranges were 510–590 and 690–820 nm.

The melanin estimates from the diffusion inverse model and objective fitting at the near infrared (NIR) wavelength range are overestimated by 16% and 15% for the lower melanin content and 1% and 3% for the higher melanin content, respectively. The inverse diffusion model has random variations in the extracted melanin parameter, while the objective fit is more or less stable. The melanin extracted from the blue–green wavelength interval for the objective fit has systematic deviations with increasing BVF in the superficial dermal layer.

The oxygenations extracted from the NIR wavelength range agree with each other across the two methods, as do the BVFs. The extracted oxygenation in the NIR wavelength range decreases by 2–3 percentage points with respect to the increased BVF in the superficial dermal layer. The oxygenation in the blue–green wavelength range does not change with respect to the BVF for the inverse diffusion model. For the objective fit, the oxygenation in the blue–green wavelength range increases for increased BVF in the upper dermal layer. The BVFs extracted from the NIR wavelength range using either method follow a small increase with increasing blood volume in the upper dermal layer.

The oxygenation in the NIR wavelength range extracted using the inverse diffusion model is lower compared to the input oxygenation in the lower dermal layer. The oxygenation in the blue–green wavelength range is higher compared to the input oxygenation in the upper dermal layer.

The BVFs in Fig. 10 follow the increase in the BVF in the upper dermal layer. The values are overshot compared to the BVF extracted from the NIR wavelength range.

Fig. 10

Output BVF from the wavelength range 510–590 nm. ${\mu }_{a,m}$ is short for ${\mu }_{a,m,694}$.

3.2.2.

Three-dimensional Monte Carlo modeling

Results of applying the hyperspectral inverse model to a 3-D Monte Carlo model of a mole are shown in Fig. 11. The 3-D vessel results are shown in Fig. 12.

Fig. 11

Melanin content (a) and oxygenation for 690–820 nm (b) extracted from the 3-D Monte Carlo model of a mole. Blood volume fraction was as in Fig. 12 outside the vessel.

Fig. 12

Blood volume fraction for 690–820 nm (a) and melanin content (b) from the 3-D Monte Carlo model of a vessel. Oxygenation for 690–820 nm was constant throughout the image, 100%.

The determined melanin values outside the mole are around $250{\mathrm{m}}^{-1}$. The determined melanin values inside the mole approach $650{\mathrm{m}}^{-1}$. Melanin input values were 225 and $1135{\mathrm{m}}^{-1}$, respectively. Increased melanin values are seen outside of the mole boundary. The mean optical penetration depth can be calculated to be approximately 0.97 mm in dermis and 0.2 mm in the epidermis of the mole area. This was calculated using the following definition of the mean optical penetration depth,⁴⁰

$$\delta =\sqrt{\frac{1}{3[(1-g){\mu }_s+{\mu }_a]{\mu }_a}},$$

and the modeled absorption and scattering coefficients at 687 nm. The BVF is about 3.7%. The oxygen saturation stays constant at 100%. Both of these parameters are extracted from the 690–820 nm wavelength interval. Some random, disconnected pixels inside the mole have a slightly lowered oxygenation saturation, to around 97%. This is likely due to increased noise levels in the mole due to higher absorption.

The determined melanin values for the vessel image stay constant to about $250{\mathrm{m}}^{-1}$. The BVFs outside the vessel are around 3.7%, while the area above the vessel approaches 5.0%. The modeled BVF in dermis was 1.0%, while the vessel was assumed to be a small tube with pure blood absorption, located at a depth of 0.25 mm. The increased BVF ranges over an area corresponding to a diameter of 0.5 mm. The oxygenation stays constant to 100% throughout the vessel. Both parameters were extracted from the 690–820 nm wavelength interval.

3.3.

Hyperspectral Image Inverse Modeling Results

Red, green, and blue (RGB) images and results of the inverse model are shown in Figs. 13 and 14, before and after 5 min of cuff-induced occlusion of the arm, respectively. Individual spectral fits are shown in Fig. 15. The images have been subsetted in order to ignore nonskin regions with a high signal-to-noise ratio and to analyze only the parts of the image that are well illuminated and in focus. The approximate width and height of the imaged subsetted area are 60 and 120 mm.

Fig. 13

Extracted skin properties from the volar side of the forearm of a Caucasian woman. These properties are the blood volume fraction in (a) 510–590 nm and (b) 690–820 nm, the oxygenation in (c) 510–590 nm and (d) 690–820 nm, (e) the melanin content and (f) the RGB image.

Fig. 14

Extracted skin properties from the volar side of the forearm of a Caucasian woman after 5 min of cuff-induced occlusion. These properties are the blood volume fraction in (a) 510–590 nm and (b) 690–820 nm, the oxygenation in (c) 510–590 nm and (d) 690–820 nm, (e) the melanin content and (f) the RGB image.

Fig. 15

Diffusion model fits of pixels located close to the mole. Spectrum 1 is before occlusion, spectrum 2 is after 5 min of occlusion. Note that the dots are used for differentiation and are not indicative of fitting points. Fitted wavelength ranges were 510–590 and 690-820 nm.

The total time used to process the subsetted ($1800\text{lines}\times 900\text{samples}$) image was 7 s.

Water and fat were fitted for the absorption spectra, but are not shown here.

Focus and illumination problems in the experimental data are propagated into an overestimation of the BVF and melanin in the upper left corner and lower right corner of the image. These artifacts in the data are due to the curvature of the arm. The melanin content was estimated to be about $600{\mathrm{m}}^{-1}$ in the central parts of the mole. The blood parameters extracted from either wavelength range show structures reminiscent of vessels.

After occlusion, the estimated blood content is increased both in normal tissue and in the blood vessels. The oxygenation is decreased at both wavelength ranges. The melanin is more or less unaffected, although the melanin in the mole is decreased compared to unaffected skin.

4.

Discussion

The aim of this paper is to present a proof of concept for a real-time inverse modeling method for hyperspectral images of skin.

The steps in the inverse modeling approach have been chosen in order to make a real-time implementation using GPU parallelization viable. A two-layered model is applied to multilayered tissue in each point of the hyperspectral image. The two-layered approach is verified through 1-D Monte Carlo simulations, while the point-based approach is verified through the 3-D Monte Carlo simulations. In addition, the model is applied to a set of basic hyperspectral images to show the potential of the technique.

Large individual variations in optical properties can be expected from measurements on in vivo tissue. The optical properties will vary as a function of temperature and hydration.^41–44 Fixed spectra of optical properties are applied in the inverse models, which have been obtained ex vivo under specific conditions and after preparation techniques, which may affect the coefficients.^43,45 As a result, there will always be a large uncertainty in the obtained tissue properties, no matter how accurate the inverse model may be when applied on numerical simulations. This has to be taken into account when setting the accuracy requirements and evaluating these kind of methods.

In general, our developed inverse diffusion model provides a mean of the tissue parameters. The resulting mean is a mean of all tissue parameters reached by the light. The model characterizes the changes of this mean. The output parameters from the simulations are stable despite the change in the BVF of the upper dermal layer. The BVF extracted using the blue–green wavelength range characterizes this change well.

The isotropic diffusion model is known to result in higher reflectance values than the corresponding Monte Carlo reflectance for the same set of input parameters and high absorption.^20,46 This will explain the slightly increased melanin values produced by the inverse diffusion model and the increased BVF extracted from the blue–green wavelength range. The difference gets higher with increased input melanin absorption. The results obtained using our methods and the objective multivariate fit are comparable for the NIR wavelength range. Mainly the BVFs are comparable for the blue–green wavelength range, while the objective fit deviates for the melanin and oxygenation. This seems to be due to cross talk, possibly between hemoglobin and melanin.

The situation for the 3-D Monte Carlo model is comparable to the results for the 1-D case. The melanin outside of the mole has the same trend as for the 1-D situation. The determined melanin content inside the mole is lowered compared to the input value, although there is no apparent misfitting. Using the estimated penetration depths as an estimate of the mean transport length, it is clear that some of the light entering the mole will exit normal skin and vice versa. This results in a washed out appearance of the mole and higher reflectance values inside the mole region. This gives the surrounding pixels the appearance of a higher melanin content. The diffuse reflectance of the mole has a slightly higher diameter (1.5 mm) than the actual mole (1 mm). This is to be expected in a highly scattering tissue.

The BVF extracted from 3D Monte Carlo is increased compared to the BVF extracted from the 1-D Monte Carlo simulation. This is not due to the difference between 1-D and 3-D Monte Carlo modelings, as the extracted properties outside both mole and vessel should be the same as for the 1-D case due to similar assumptions. The source of this difference might be the SCA method, which was used for unmixing of the absorption. SCA (sequential coordinate-wise algorithm for non-negative least-squares problems) was chosen for the unmixing of the derived absorption coefficients in the hyperspectral images due to its suitability for GPU implementation. This method was not used in unmixing of the absorptions in the 1-D case, where an ordinary non-negative least-squares algorithm was applied. The use of SCA apparently also results in increased oxygenation. The oxygenation should be slightly lower than the input oxygenation of 95%. This can be traced to a cross talk between the constant baseline absorption and the blood absorptions, where SCA encounters more challenges than an ordinary non-negative least-squares algorithm.

The changes are still characterized well. The properties do not change significantly across the mole, and an increased BVF is seen for the vessel area. Many scattering events give a washed-out appearance, which makes complete reconstruction difficult. We are still able to extract useful information from the model.

The same behavior is seen for the experimental data. The properties have a washed-out, diffuse look due to lateral broadening. The characterized properties are, on the other hand, more realistic than the properties extracted from the 3-D Monte Carlo model. The oxygenation saturation for the NIR wavelength range does not approach 100% and the BVFs are not unrealistically high. It seems that the cross talk between the constant baseline absorption and the blood parameters is not seen here.

With the results from the simulated mole in mind, we can expect the melanin content to be underestimated inside the real mole. This is seen. The other extracted parameters seem to be affected by this. They were not affected for the 3-D Monte Carlo. This seems to be mainly due to a change in scattering parameters. Compared to normal skin, scattering properties are altered due to the structural difference of a mole. The method does not rectify variations in scattering, though scattering variations may be observed in the more extreme cases by quantifying the misfit between the different wavelength ranges.

The trends of the extracted parameters are as expected from the experiences with the simulations. The BVFs are increased where blood vessels are localized. The vessels extracted from the NIR wavelength range are likely to carry venous blood since the depths of these are less than the depths of the large arteries. The oxygenation here is lowered compared to the oxygenation in the rest of the tissue, which is assumed to be correct. The extracted oxygenations are decreased after occlusion of the arm, which is expected.

In general, the oxygenation extracted from the blue–green wavelength interval is lower than the oxygenation extracted from the NIR wavelength interval. This agrees with results found by Tseng et al.,²⁹ and is likely due to differences in vascularization down to the different penetration depths.²¹

The melanin is slightly underestimated for lower oxygenations and at the locations of larger blood vessels. The former is mainly due to cross talk with deoxyhemoglobin. The latter may be due to both cross talk and changes in the boundary condition assumptions (e.g., changes in refraction indices).

SCA was used in the unmixing process. While having no convergence guarantees, it still minimizes the differences between the fitted and derived absorption coefficients. This is evident from the displayed spectra fits in Fig. 15. Still, it should not be trusted for the unmixing of larger wavelength intervals and many chromophores. The experiences with the Monte Carlo simulations showed SCA to be less trustworthy than a more ordinary non-negative least-squares algorithm. SCA has, however, the advantage of being more suitable for GPU implementation. Future work will involve improving the unmixing algorithms. This can be done either by improving SCA, adapting other non-negative least-squares GPU implementations,⁴⁷ or by implementing a non-negative least squares algorithm optimized for problems such as the unmixing of hyperspectral imagery.⁴⁸

Illumination problems in the image lead to some artifacts in the extracted parameters. This issue in the imaging technique is currently being addressed in another study which aims to use 3-D modeling to obtain a digital elevation model and flatten the image.

The developed model is a proof of concept where we have shown the possibility of characterizing spatially resolved tissue properties in real-time. The model will still be further developed to obtain more complexity and accuracy. The isotropic source functions used in the model will be exchanged by more accurate Delta-Eddington source functions.⁴⁶ No special assumptions were made for the blood vessel distribution in the skin model. It is known that assuming an average blood vessel diameter will affect the apparent absorption levels in blood.⁴⁹ The lack of such a correction of the blood absorption may have affected the extracted blood oxygenation and BVF, and is something which will be implemented in the future.

The initial results obtained from the hyperspectral images show promise in the characterization of tissue properties. The method can be used to identify interesting areas during image scan, such as areas where the optical properties are changed (i.e., wounds, moles, bruises, or other skin lesions). Classification and statistical methods can be run on the estimated tissue properties to automatize this identification process. The interesting areas may then be more closely investigated.

5.

Conclusion

An inverse photon transport model with real-time performance has been developed for a hyperspectral image scanning system using GPU parallelization. The model fulfills the real-time analysis constraints set by the hyperspectral setup, leaving a lot of computational time for additional image processing.

Simulations have shown that the inverse model has an ability to characterize changes in optical properties. Running the inverse model on hyperspectral images of skin shows promising results.

Future work will involve improvement of the unmixing algorithms and complexity of the models.