## 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.^{5} 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.^{6} 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.^{9} 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^{10} or the diffusion approximation,^{11}12.^{–}^{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}14.^{–}^{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.1.

### Background

The main objective of the study was to develop an algorithm for extraction of skin properties from hyperspectral images in real-time during scanning. This has been done using GPU parallelization. GPU hardware is essentially an SIMD (single instruction and multiple data) vector processor.^{18} For maximum parallelization, this requires the same sequence of instructions to be independently applied to multiple data. Since the hyperspectral data are discretized in pixels and wavelengths, doing individual processing on pixels and wavelengths is the natural way to parallelize the GPU processing.

A photon transport model is iterated with respect to skin parameters, such as blood volume fraction and melanin concentration. In order to apply the same sequence of GPU instructions for all data, the iteration strategies will be limited to rather simple iteration strategies. To fulfill the independency requirement, each pixel in the hyperspectral image is treated as an independent diffuse reflectance spectrum measured on a laterally infinite layered medium. It is possible to achieve pixel interdependency and nondeterminism, but at the cost of reducing GPU processing optimality. This is not applied in this preliminary study.

The GPU code was designed for a computer with an Intel Core i7 CPU (8 cores), 6 GB RAM and an NVIDIA GeForce GTX 670 GPU. Debian GNU/Linux (jessie) was used as the operating system. The code was integrated into a hyperspectral streaming framework developed by the Norwegian Defence Research Establishment.^{19}

## 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.,^{13} Randeberg et al.,^{20}^{,}^{21} and Svaasand et al.^{11} 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.^{13} and is given as

## (1)

$${\mu}_{a,m}(\lambda )={\mu}_{a,m,694}\xb7{\left(\frac{\lambda}{694\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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.^{22} It depends on the blood volume fraction (BVF) and oxygen saturation (oxy). A constant background absorption of $25\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$ was added to the model.^{11} Absorption spectra of hemoglobin, fat, and melanin are shown in Fig. 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^{25}

## (2)

$${\mu}_{s}^{\prime}(\lambda )=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}\xb7[{a}_{\mathrm{Mie}}{\left(\frac{\lambda}{500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}}\right)}^{-{b}_{\mathrm{Mie}}}+{a}_{\mathrm{Ray}}{\left(\frac{\lambda}{500\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}}\right)}^{-4}].$$Bashkatov et al.^{26} 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

^{27}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.^{11} 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.^{20}

The error between the Monte Carlo model and the diffusion approximation was reported by Randeberg et al.^{20} 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.^{20} 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.^{7} 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^{28} 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.^{21} This approach is also similar to work done by Tseng et al.^{29} and Saager et al.^{30} 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.

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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$ as a starting value.

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

The sequential coordinate-wise algorithm for non-negative least-squares problems (SCA)^{31} 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.^{32}

## (5)

$${\mu}_{a,b}(\lambda )=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}\xb7(0.82+16.82{e}^{-(\lambda -400\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm})/80.5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}}).$$The reduced scattering coefficient in subcutis was modeled by an expression given by Naglic et al.,^{33} comparable to the scattering coefficients reported by Salomatina et al.,^{32}

## (6)

$${\mu}_{s}^{\prime}(\lambda )=1500{\mathrm{m}}^{-1}\xb7(16.34+303.8{e}^{-\lambda /180.3\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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^{7} 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.

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^{34} 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^{35} 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.

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

## (7)

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

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.

## 2.3.

### Experimental

The developed inverse diffusion model was also tested on measured hyperspectral data.

Hyperspectral images of skin were collected using a push-broom HySpex VNIR-1600 camera (Norsk Elektro Optikk, Lillestrøm, Norway).^{38} A healthy, female volunteer (Caucasian, 39 years old) with fair skin (Fitzpatrick skin type I/II) had the volar side of her forearm imaged. Two images were obtained. The first image was obtained as a baseline. The second image was obtained after 5 min of occlusion induced by applying a blood pressure cuff. This was done in order to modulate the oxygenation and blood content in the skin. The images had a size of $1600\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{pixels}\text{\hspace{0.17em}}(\text{samples})\times 160\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{wavelengths}\text{\hspace{0.17em}}(\text{bands})\times \phantom{\rule{0ex}{0ex}}\text{a varying number of lines}$. Hyperspectral data lines were scanned at a speed of 30 ms per line of data, which was also chosen to be the real-time deadline limit for the processing.

The lens had a focal length of 30 cm. The pixel field of view was approximately 0.4 mrad.^{38} Pixel size on the skin surface after magnification using a 30 cm lens was approximately $60\times 60\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{m}$.

Two linear light sources were used for illumination (Model 2900 Tungsten Halogen, Illumination Technologies, New York). Polarizers were mounted on the camera lens and the light sources (VLR-100 NIR, 450–1100 nm, Meadowlark Optics, Frederick, Colorado) in order to avoid specular reflection.

The images were converted into reflectance and corrected for uneven illumination across the field of view using a Spectralon reflectance target (SRT-50-050 Reflectance Target, $12.7\times 12.7\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{cm}$, ACAL Bfi Nordic AB, Uppsala).^{6} Spectral variations in the specified intensity of the reflectance standard were taken into account in the conversion. The images were denoised using the maximum noise fraction transform (MNF).^{39}

## 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.

## 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 μa,m,694 (m−1) | Output μa,m,694 (m−1) | oxy510−590 | oxy690−820 | BVF690−820 |
---|---|---|---|---|---|

Inv. DM | 227 | 264±7 | 0.65±0.01 | 0.77±0.01 | 0.026±0.002 |

Inv. DM | 683 | 692±11 | 0.67±0.01 | 0.76±0.03 | 0.026±0.002 |

Obj. Fit | 227 | 230±39 | 0.66±0.02 | 0.76±0.02 | 0.025±0.002 |

Obj. Fit | 683 | 682±31 | 0.70±0.03 | 0.77±0.02 | 0.028±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 | ||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

μa,m | BVF1 | μa,m | oxy1 | BVF1 | oxy2 | BVF2 | μa,m,1 | oxy1 | BVF1 | μa,m,2 | oxy2 | BVF2 |

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 |

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.

## 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.

The determined melanin values outside the mole are around $250\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$. The determined melanin values inside the mole approach $650\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$. Melanin input values were 225 and $1135\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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,^{40}

The determined melanin values for the vessel image stay constant to about $250\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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.

The total time used to process the subsetted ($1800\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{lines}\times \phantom{\rule{0ex}{0ex}}900\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\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}42.43.^{–}^{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.,^{29} and is likely due to differences in vascularization down to the different penetration depths.^{21}

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,^{47} or by implementing a non-negative least squares algorithm optimized for problems such as the unmixing of hyperspectral imagery.^{48}

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.^{46} 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.^{49} 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.

## Acknowledgments

Thanks to Norsk Elektro Optikk and the Norwegian Defence Research Establishment for exchange of code. The Norwegian Research School in Medical Imaging has provided funding for this project through the MedIm Bridging Grant. Thanks to Lukasz Paluchowski for help with the experimental work.

## References

G. LuB. Fei, “Medical hyperspectral imaging: a review,” J. Biomed. Opt. 19(1), 010901 (2014).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.19.1.010901Google Scholar

T. Skauliet al., “A compact combined hyperspectral and polarimetric imager,” Proc. SPIE 6395, 639505 (2006).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.693484Google Scholar

L. L. RandebergE. L. P. LarsenL. O. Svaasand, “Characterization of vascular structures and skin bruises using hyperspectral imaging, image analysis and diffusion theory,” J. Biophotonics 3(1–2), 53–65 (2010).JBOIBX1864-063Xhttp://dx.doi.org/10.1002/jbio.200910059Google Scholar

E. L. Larsenet al., “Hyperspectral imaging of atherosclerotic plaques in vitro,” J. Biomed. Opt. 16(2), 026011 (2011).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3540657Google Scholar

J. Hernandez-Palacioset al., “Hyperspectral characterization of fluorophore diffusion in human skin using a scmos based hyperspectral camera,” Proc. SPIE 8087, 808717 (2011).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.889816Google Scholar

M. Denstedtet al., “Hyperspectral imaging as a diagnostic tool for chronic skin ulcers,” Proc. SPIE 8565, 85650N (2013).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.2001087Google Scholar

E. Alerstamet al., “Next-generation acceleration and code optimization for light transport in turbid media using gpus,” Biomed. Opt. Express 1(2), 658–675 (2010).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.1.000658Google Scholar

Y. Tarabalkaet al., “Real-time anomaly detection in hyperspectral iimage using multivariate normal mixture models and gpu processing,” J. Real-Time Image Proc. 4(3), 287–300 (2009).1861-8200http://dx.doi.org/10.1007/s11554-008-0105-xGoogle Scholar

N. KeshavaJ. Mustard, “Spectral unmixing,” IEEE Signal Proc Mag. 19, 44–57 (2002).ISPRE61053-5888http://dx.doi.org/10.1109/79.974727Google Scholar

L. WangS. L. JacquesL. Zheng, “Mcml monte carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47(2), 131–146 (1995).CMPBEK0169-2607http://dx.doi.org/10.1016/0169-2607(95)01640-FGoogle Scholar

L. Svaasandet al., “Tissue parameters determining the visual appearance of normal skin and port-wine stains,” Laser Med. Sci. 10, 55–65 (1995).LMSCEZ1435-604Xhttp://dx.doi.org/10.1007/BF02133165Google Scholar

R. C. Haskellet al., “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. A 11, 2727–2741 (1994).JOAOD60740-3232http://dx.doi.org/10.1364/JOSAA.11.002727Google Scholar

T. Spottet al., “Application of optical diffusion theory to transcutaneous bilirubinometry,” Proc. SPIE 3195, 234–245 (1998).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.297907Google Scholar

L. L. Randeberget al., “A novel approach to age determination of traumatic injuries by reflectance spectroscopy,” Laser Surg. Med. 38(4), 277–289 (2006).LSMEDI0196-8092http://dx.doi.org/10.1002/(ISSN)1096-9101Google Scholar

R. Zhanget al., “Determination of human skin optical properties from spectrophotometric measurements based on optimization by genetic algorithms,” J. Biomed. Opt. 10(2), 024030 (2005).JBOPFO1083-3668http://dx.doi.org/10.1117/1.1891147Google Scholar

T. Tsenget al., “Quantification of the optical properties of two-layered turbid media by simultaneously analyzing the spectral and spatial information of steady-state diffuse reflectance spectroscopy,” Biomed. Opt. Express 2(4), 901–914 (2011).BOEICL2156-7085http://dx.doi.org/10.1364/BOE.2.000901Google Scholar

H. CenR. Lu, “Quantification of the optical properties of two-layer turbid materials using a hyperspectral imaging-based spatially-resolved technique,” Appl. Opt. 48(29), 5612–5623 (2009).APOPAI0003-6935http://dx.doi.org/10.1364/AO.48.005612Google Scholar

CUDA Toolkit Documentation, CUDA C Programming Guide, 2012, http://docs.nvidia.com/cuda (8 May 2014).Google Scholar

T. Skauliet al., “An airborne real-time hyperspectral target detection system,” Proc. SPIE 7695, 76950A (2010).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.850443Google Scholar

L. L. Randeberget al., “Performance of diffusion theory vs. monte carlo methods,” Proc. SPIE 5862, 58620O (2005).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.633028Google Scholar

L. L. Randeberget al., “In vivo spectroscopy of jaundiced newborn skin reveals more than a bilirubin index,” Acta Paediatr. 94(1), 65–71 (2005).APAEEL0803-5253http://dx.doi.org/10.1080/08035250410023179Google Scholar

W. G. ZijlstraA. BuursmaO. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin, VSP, Utrecht (2000).Google Scholar

H. BuiteveldJ. M. H. HakvoortM. Donze, “The optical properties of pure water,” Proc. SPIE 2258, 174–183 (1994).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.190060Google Scholar

R. L. P. van Veenet al., “Determination of vis-nir absorption coefficients of mammalian fat, with time- and spatially resolved diffuse reflectance and transmission spectroscopy,” in Proc. Biomedical Topical Meeting, SF4, Optical Society of America, Washington DC (2004).Google Scholar

S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R59 (2013).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/58/11/R37Google Scholar

A. N. Bashkatovet al., “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys. 38, 2543–2555 (2005).JPAPBE0022-3727http://dx.doi.org/10.1088/0022-3727/38/15/004Google Scholar

M. J. C. van Gemertet al., “Skin optics,” IEEE Trans. Biomed. Eng. 36, 1146–1154 (1989).IEBEAX0018-9294http://dx.doi.org/10.1109/10.42108Google Scholar

W. H. Presset al., Numerical Recipes, Cambridge University Press, Cambridge, UK (2007).Google Scholar

S. Tsenget al., “Chromophore concentrations, absorption and scattering properties of human skin in-vivo,” Opt. Express 17, 14599–14617 (2009).OPEXFF1094-4087http://dx.doi.org/10.1364/OE.17.014599Google Scholar

R. B. Saageret al., “Method for depth-resolved quantitation of optical properties in layered media using spatially modulated quantitative spectroscopy,” J. Biomed. Opt. 16(7), 077002 (2011).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3597621Google Scholar

V. FrancV. HlavacM. Navara, “Sequential coordinate-wise algorithm for the non-negative least squares problem,” in Proc. Computer Analysis of Images and Patterns, 11th International Conference, CAIP 2005, Versailles, France, September 5-8, 2005. Lecture Notes in Computer Science, Vol. 3691, pp. 407–414, Springer, New York (2005).Google Scholar

E. Salomatinaet al., “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11(6), 064026 (2006).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2398928Google Scholar

P. Naglicet al., “Applicability of diffusion approximation in analysis of diffuse reflectance spectra from healthy human skin,” Proc. SPIE 9032, 90320N (2013).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.2044706Google Scholar

R. J. HansonC. L. Lawson, Solving Least Squares Problems, Society for Industrial and Applied Mathematics, Philadelphia (1995).Google Scholar

M. MilanicB. Majaron, “Three-dimensional monte carlo model of pulsed-laser treatment of cutaneous vascular lesions,” J. Biomed. Opt. 16(12), 128002 (2011).JBOPFO1083-3668http://dx.doi.org/10.1117/1.3659205Google Scholar

M. Friebelet al., “Determination of optical properties of human blood in the spectral range 250 to 1100 nm using Monte Carlo simulation with hematocrit-dependence effective scattering phase functions,” J. Biomed. Opt. 11(3), 034021 (2006).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2203659Google Scholar

H. LiL. LinS. Xie, “Refractive index of human whole blood with different types in the visible and near-infrared ranges,” Proc. SPIE 3914, 517–521 (2000).PSISDG0277-786Xhttp://dx.doi.org/10.1117/12.388073Google Scholar

Hyspex VNIR-1600, “Main specifications,” http://www.hyspex.no/products/hyspex/vnir1600.php (8 May 2014).Google Scholar

A. A. Greenet al., “A transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE T. Geosci. Remote 26, 65–74 (1988).IGRSD20196-2892http://dx.doi.org/10.1109/36.3001Google Scholar

L. V. WangH. Wu, Biomedical Optics, Principles and Imaging, John Wiley & Sons, Hoboken, New Jersey (2007).Google Scholar

J. Lauferet al., “Effect of temperature on the optical properties of ex vivo human dermis and subdermis,” Phys. Med. Biol. 43 2479–2489 (1998).PHMBA70031-9155http://dx.doi.org/10.1088/0031-9155/43/9/004Google Scholar

O. S. Khalilet al., “Temperature modulation of the visible and near infrared absorption and scattering coefficients of human skin,” J. Biomed. Opt. 8(2), 191–205 (2003).JBOPFO1083-3668http://dx.doi.org/10.1117/1.1559997Google Scholar

T. ListerP. A. WrightP. H. Chappell, “Optical properties of human skin,” J. Biomed. Opt. 17(9), 090901 (2012).JBOPFO1083-3668http://dx.doi.org/10.1117/1.JBO.17.9.090901Google Scholar

C. G. Rylanderet al., “Dehydration mechanism of optical clearing in tissue,” J. Biomed. Opt. 11(4), 041117 (2006).JBOPFO1083-3668http://dx.doi.org/10.1117/1.2343208Google Scholar

E. K. Chanet al., “Effects of compression on soft tissue optical properties,” IEEE J. Sel. Top. Quant. 2(4), 943–950 (1996).IJSQEN1077-260Xhttp://dx.doi.org/10.1109/2944.577320Google Scholar

T. SpottL. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. 39, 6453–6465 (2000).APOPAI0003-6935http://dx.doi.org/10.1364/AO.39.006453Google Scholar

Y. LuoR. Duraiswami, “Efficient parallel nonnegative least squares on multicore architectures,” SIAM J. Sci. Comput. 33, 2848–2863 (2011).SJOCE31064-8275http://dx.doi.org/10.1137/100799083Google Scholar

M. H. Van BenthemM. R. Keenan, “Fast algorithm for the solution of large-scale non-negativity-constrained least squares problems,” J. Chemometrics 18, 441–450 (2004).JOCHEU0886-9383http://dx.doi.org/10.1002/(ISSN)1099-128XGoogle Scholar

L. O. Svaasandet al., “Therapeutic response during pulsed laser treatment of port-wine stains: dependence on vessel diameter and depth in dermis,” Laser Med. Sci. 10, 235–243 (1995).LMSCEZ1435-604Xhttp://dx.doi.org/10.1007/BF02133615Google Scholar