Open Access
Add to BibSonomyAbstract
Hyperspectral images of tissue contain extensive and complex information relevant for clinical applications. In this work, wavelet decomposition is explored for feature extraction from such data. Wavelet methods are simple and computationally effective, and can be implemented in real-time. The aim of this study was to correlate results from wavelet decomposition in the spectral domain with physical parameters (tissue oxygenation, blood and melanin content). Wavelet decomposition was tested on Monte Carlo simulations, measurements of a tissue phantom and hyperspectral data from a human volunteer during an occlusion experiment. Reflectance spectra were decomposed, and the coefficients were correlated to tissue parameters. This approach was used to identify wavelet components that can be utilized to map levels of blood, melanin and oxygen saturation. The results show a significant correlation (p <0.02) between the chosen tissue parameters and the selected wavelet components. The tissue parameters could be mapped using a subset of the calculated components due to redundancy in spectral information. Vessel structures are well visualized. Wavelet analysis appears as a promising tool for extraction of spectral features in skin. Future studies will aim at developing quantitative mapping of optical properties based on wavelet decomposition.
© 2014 Optical Society of America
1. Introduction
Hyperspectral imaging, which combines spectroscopy with optical imaging, has been proposed and employed as a diagnostic tool for clinical examination of tissue [1–3]. Hyperspectral imaging is essentially spatially resolved diffuse reflectance spectroscopy. Diffuse reflectance spectra in the visible to near infra-red region contain a vast amount of information, originating from the absorption and scattering properties of the tissue constituents. Such information may be valuable in characterization and diagnosis of tissue in applications such as e.g. monitoring of wound healing, cutaneous burns and radio therapy induced skin changes. Although the properties of the individual chromophores are well known, it is challenging to extract accurate information on chromophore concentrations and distributions from measured spectra due to the scattering nature of tissue. Hyperspectral data and reflectance spectra are currently processed by simple wavelength ratios, statistical image analysis or by using analytical or numerical methods such as diffusion theory or Monte Carlo techniques.
The use of band ratio metrics can provide some information about the concentration of a chromophore (e.g. oxygenated hemoglobin [2]) or a property of interest (e.g. erythema index [1]). This is a fast and simple approach which only uses the information from a few wavelength bands, but does not take advantage of the full potential in hyperspectral imaging. The method is prone to error if other or more chromophores than expected are present.
Statistical classification is another method [3–7]. This approach makes use of the full spectral information, which reduces the problem of cross-talk. However, due to the composition of tissue, reflectance spectra from qualitatively different areas may still be statistically similar. A global classification scheme might fail to catch more subtle features when the overall dynamics of the spectrum is high. Numerical methods are often considered as the gold standard in biomedical optics. However, such methods might be computationally demanding and thus time consuming. This makes them less applicable for real-time diagnostics. Due to all these considerations, there is still a need for investigating new techniques to identify physically informed methods that are suitable for real-time implementation.
The inherent challenge in characterizing tissue spectra is that important information appears with varying scale and spectral locations. Wavelet analysis is designed for this exact situation. The transform provides a means to analyze signals on different levels of detail, without losing information about locations of the isolated features. Wavelet analysis has previously proven efficient in dimensionality reduction and classification in agricultural hyperspectral imaging [8]. The method is fast and not computationally demanding, which makes it applicable for real-time implementation.
In this work, the use of wavelet decomposition is explored for feature extraction from hyper-spectral images of tissue. The main aim of the study was to correlate the results from wavelet decomposition in the spectral domain with physical parameters such as tissue oxygenation, blood and melanin content. Correlations between wavelet components and skin properties are investigated through Monte Carlo simulations of tissue spectra, measurements of a tissue phantom and hyperspectral data acquired during an occlusion experiment.
2. Theory and methods
2.1. Hyperspectral imaging
Hyperspectral imaging provides data with high spatial and spectral resolution. In every pixel of a hyperspectral image, light intensities are recorded for numerous wavelengths, forming a 3-D data cube spanning two spatial dimensions and one spectral. In this study, a pushbroom-scanning hyperspectral camera (VNIR 1600, Norsk Elektro Optikk AS [9]) was used to collect in vivo reflectance data from human skin. The camera covers the wavelength range 400 – 1000 nm and has a spectral resolution of 3.6 nm (160 bands). The width of the imaged area (FOV) was approximately 10 cm, with a spatial resolution of 1600 pixels. Each image pixel had an approximate area of 0.0035 (mm)2 at the surface of the sample. The scanning speed was about 0.5 cm per second with an integration time of 30 ms per line.
The scene was illuminated by a halogen lamp (Model 2900 Tungsten Halogen, Illumination Technologies), projecting two lines of light, overlapping at the focal distance. To reduce specular reflectance, linear polarizers (VersaLight polarizer, VLR-100 NIR, 450–1100 nm) were placed in front of the camera lens and each of the light sources (cross-positioned). A gray reflectance standard (Labsphere Inc., SRT-50-050 Reflectance Target) was placed in the scene, and flat field calibration was performed in the software pre-processing of data, transforming light intensities from radiance to reflectance. Prior to further analysis, the image noise was reduced using the Minimum Noise Fraction algorithm (MNF) [10].
2.2. Tissue optics
The main chromophores in skin are hemoglobin and melanin. The characteristics of the respective absorption spectra can be used to identify or quantify the presence and relative concentrations of skin chromophores in the analysis of collected hyperspectral images. There is a complex and non-linear relation between the absorption and the resulting diffuse reflectance. Features observed as peaks in an absorption spectrum tend to appear as minima in the reflectance. However, if several chromophores with overlapping spectra, layered structures and multiple scattering is present, the picture becomes more complicated. The relation between optical parameters and observed reflectance can in such situations be explored by the use of numerical methods e.g. Monte Carlo simulations. Figure 1 show the absorption spectra of the most important chromophores, and how they relate to the observed reflectance.
Fig. 1 (a): Absorption coefficients of oxygenated and deoxygenated hemoglobin [14] and melanin [11]. (b): Absorption coefficients of water [15] and mammalian fat [16]. Note that the two figures have different scales. Sample reflectance generated using a Monte Carlo model is shown in (c). See also Fig. 7 and 11 for further examples of diffuse reflectance spectra. The scattering models are plotted in (d).
Oxygenated hemoglobin (HbO2) exhibits absorption maxima at 415, 542 and 576 nm, whereas deoxygenated hemoglobin (Hb) has peaks at 550 and 760 nm (Fig. 1(a)).
Melanin absorption is substantial in pigmented skin, and decreases exponentially with increasing wavelength in the VIS-NIR range [11], as can be seen in Fig. 1(a).
Water has a significant absorption peak at 970 nm, which can be detected in the NIR region of reflectance spectra from tissue. Fat (lipids) absorbs in the blue/green region and in the NIR part of the spectrum, and is notable at 930 nm. The absorption spectra of water and fat are shown in Fig. 1(b).
The other important optical process that takes place in light-tissue interaction, is scattering. Scattering occurs where there are differences in the refractive index, and is caused by e.g. lipids, water and collagen. These interactions can be described as a combination of Mie and Rayleigh scattering and is commonly modeled by a scattering coefficient μs, and the anisotropy factor g (average cosine of scattering angle), which are both wavelength dependent [17]. Figure 1 shows the scattering properties employed in this study. A further description of these models is given in section 2.6.
2.3. Wavelet analysis
The reflectance spectrum in each individual pixel can be regarded as a continuous function of wavelength, R(λ), sampled at regular intervals, . These signals are characterized by localized features on varying scales, which makes them suitable for Multi Resolution Analysis (MRA) (see e.g. [18, Ch. 5, p. 190ff]), such as the Discrete Wavelet Transform (DWT).
The DWT of {rn} at scale j and location k, is calculated through discrete convolution with a scaled wavelet function ψj, and can be expressed as
where k is the translation variable and An efficient implementation of the DWT, is the dyadic filter bank algorithm [19], in which the wavelet basis is represented by a set of orthogonal high-pass and low-pass filters. In the first step of the algorithm, the original signal r is decomposed into an approximation, a1, and a detail signal, d1, through convolution with the low-pass filter h and the high-pass filter g, respectively, followed by down-sampling by a factor of 2. This procedure is then cascaded and successively applied to the approximation signals, producing ever coarser approximations and details in each iteration. This is illustrated in Fig. 2. The maximum number of iterations is given by the length of the analyzed signal and the filter length used.Fig. 2 High and low-pass filtering (filters g and h, respectively) followed by down-sampling, produces approximation and detail signals. In the first step of the algorithm, the original signal is decomposed, then the process is iterated and repeatedly applied to the new approximation.
Figure 3 shows the decomposition of a sample reflectance spectrum. The signals shown have been reconstructed from the corresponding approximation and detail coefficients through inverse DWT [19].
Fig. 3 Reflectance signal r decomposed in 6 iterations using the high-pass filter g and low-pass filter h associated with the Symlet wavelet [20]. Reconstructed approximation and detail components are shown. The coefficients of all levels of detail and the coarsest approximation are kept in the implemented algorithm (marked by the red frame).
In the process of choosing an appropriate wavelet family, a number of properties need to be considered. First of all, the wavelet must be designed for discrete analysis, i.e. the filter functions have to be discretized to a finite number of points. If energy preservation is required, which is a prerequisite for exact reconstruction, the wavelet basis functions have to be orthogonal. The support width determines how well features can be localized in time-space (or wavelength in this case), so compact support is a desirable property.
In order to avoid introducing artifacts, the wavelet should be smooth (high number of derivatives) and preferably symmetric (linear phase). Orthogonal wavelets cannot be symmetric, but may be optimized to be as close to symmetry as possible. Another alternative is to construct biorthogonal wavelets, which can be symmetric, but with the loss of energy preservation. Biorthogonal wavelets provide an invertible transform, and are optimized to minimize the deviation from orthonormality.
The number of vanishing moments can be thought of as a measure of a wavelet’s ability to represent a complex function. The more vanishing moments, the sparser the set of coefficients needed to represent a signal. A high number of vanishing moments, requires a larger support width, so there is a trade-off between these two properties.
Finally, if a wavelet already resembles the signal to be decomposed, or parts of it, it is likely to capture features efficiently, i.e. most of the energy ends up in a few large wavelet coefficients. To be able to capture the spectral features of hemoglobin the wavelet should resemble the spectral characteristics of the hemoglobin spectrum.
Several wavelet families satisfy the pure mathematical requirements, e.g. Daubechies wavelets, Symlets, Coiflets and Biorthogonal wavelets. Several of these filters generate interesting decompositions that may be used to analyze tissue signals. However, when the decomposition components are examined, the more symmetric Symlets and Coiflets tend to render more high contrast images. In other words, these filters separate characteristics more efficiently for the current data sets. The results of using either of these two wavelet types are very similar for short filter lengths. The Symlet [20] filters of length 8 were chosen in this study. The filters are orthogonal with 4 vanishing moments, near linear-phase and provide a reasonable trade-off between information detail and efficiency. This choice fulfills the required mathematical properties and resemble the blood absorption features of a typical skin reflectance spectrum (see Fig. 1). The wavelet function is shown in Fig. 4.
Fig. 4 The symlet-4 wavelet function. Its resemblance to the blood absorption features of skin reflectance ensures that important features in the reflectance are captured efficiently.
Six iterations of the above described algorithm provides one approximation and 6 sets of detail coefficients. This was done in every pixel of the collected images, resulting in a total of 199 component images.
2.4. In vivo imaging
Hyperspectral image data were collected in vivo from the volar part of the left forearm of a 40 year old Caucasian female (fair skin, Fitzpatrick skin type I/II). Blood flow was occluded using a strap around the upper arm. The arm was imaged before, during (after 5 minutes of occlusion) and three minutes after occlusion.
Applying external pressure primarily causes occlusion of the superficial veins, and not as much of the deeper arteries. This leads to a net increase of blood in the limb. Eventually, it also induces a lower oxygen saturation level, as the drainage of deoxygenated venous blood is obstructed and the inflow of arterial oxygen rich blood is limited by the high pressure in the limb.
The data were converted to reflectance and underwent noise removal using MNF-filtration [10]. Reflectance spectra were then extracted from three qualitatively different regions of the imaged limb (normal skin, mole, blood vessel). All spectra were averaged over 10 by 10 pixels.
The three sets of data from the different localizations can be compared to investigate e.g. how increase/decrease in blood content as well as oxygenation affects the spectral properties in different regions of the imaged area.
Oxygenation and blood volume fraction were determined for a region of interest from each image using an inverse diffusion model [21]. This inverse model yields relative differences of the derived properties down to a shallow penetration depth (blue/green wavelength interval) and a deep penetration depth (NIR wavelength interval).
DWT was run on the hyperspectral images, and the extracted DWT components were compared to DWT components obtained from simulations and phantom data. For visualization, individual DWT coefficients were displayed with a fixed intensity range in order to be able to compare coefficients from different images.
2.5. Phantom data
Tissue phantom data from a previous study [22] were used to investigate the effect of depth and oxygenation on the DWT coefficients. The phantom setup is shown in Fig. 5. Blood vessels at various depths in skin tissue are emulated by three capillary tubes filled with whole blood, submerged with a tilt in a 10 % intralipid/water solution. Two of the tubes were filled with oxygenated blood (close to 99% oxygenated), while the third tube was deoxygenated to below 30% (measured to 22%) using Sodium dithionite (Na2S2O4).
Fig. 5 Phantom design. The phantom design was previously published by Randeberg et al. [22]. Capillary tubes filled with blood were placed in a bath filled with a 10% intralipid/water solution.
2.6. Monte Carlo simulations
Monte Carlo models were used to simulate single reflectance spectra and larger hyperspectral images in order to investigate the link between DWT components and physical skin properties. Two different models were used for this: A 1D (one-dimensional) model with a layered, laterally homogeneous geometry and a 3D model with more complex geometry. The 1D model was used for correlation of the DWT components to the input parameters, while the 3D model was used to test the methods in a non-homogeneous, known situation.
2.6.1. 1D Monte Carlo modeling
Simulated spectra were obtained using a GPU (Graphics Processing Unit) accelerated Monte Carlo code [23].
The skin model was based on a skin model used in a previous study [21]. The geometry was composed of 3 laterally infinite layers. The uppermost layer (epidermis) had a thickness of 60 μm and a default melanin content of 3%. The second layer (dermis) was 1.5 mm thick with a default blood concentration of 1% and oxygenation of 80%. The third, semi-infinite layer (subcutis) contained 40% fat and the background absorption given by Salomatina et al. [24],
Epidermis and dermis had a fixed background absorption of 25 m−1.The reduced scattering in epidermis and dermis was based on an expression given by Jacques [25],
The parameter values were set to fRayleigh = 0.48, a = 36400 m−1 and bMie = 0.22, based on data obtained by Bashkatov et al. [12,26] (scattering alternative I). Simulations were also run for one other set of scattering parameters in order to test the robustness of the developed approach. These parameters were obtained from data presented by Anderson et al. [13]: fRayleigh = 0.22, a = 45600 m−1 and bMie = 1.184 (scattering alternative II).The reduced scattering coefficient in subcutis was modeled by an expression given by Naglic et al. [27]:
Simulations were run for every wavelength from 415 to 992 nm in 3.6 nm intervals, resulting in 160 data points for each spectrum. Each simulation had 3 000 000 photons per wavelength.
To investigate the effect of varying dermal blood fraction, epidermal melanin content and oxygenation, three sets of 5 spectra were simulated. The simulation parameters are specified in table 1.
Table 1. Parameter combinations of the simulated spectra. All spectra were simulated using both scattering models.
All simulated spectra were subsequently decomposed using wavelets, as described in section 2.3. Correlation coefficients and p-values between input parameters (levels of blood content, melanin and oxygen saturation) and output DWT components were calculated. The correlation coefficients were used as guidance in choosing DWT components, and only images which could intuitively be interpreted (e.g. high blood vessel contrast) were chosen for further analysis.
2.6.2. 3D Monte Carlo modeling
3D Monte Carlo simulations were run in order to test the developed approach on a hyperspectral dataset with known tissue geometry. A GPU-accelerated version of the model presented by Milanic et al. [28] was used. The geometry was implemented as discrete voxels assigned tissue optical properties. The light source was implemented as a homogeneous plane wave. Lateral mirror boundary conditions were enforced, while subcutis was treated as a semi-infinite layer. Simulations were run for 160 wavelengths from 415 to 992 nm with a step size of 3.6 nm. Noise removal was run on the simulated data using the MNF transform. Optical properties in the homogeneous parts of the models were set to the same as the 1D simulations unless otherwise specified.
Three different geometries were simulated: a mole, large subcutaneous vessels and smaller vessels going from superficial layers to deeper layers. Vessels going across dermis from superficial- to deeper layers is unphysiological, but were modeled in this way in order to match the phantom setup.
The mole was modeled as a disc of 0.5 mm radius with 5% melanin, located in the center of an otherwise homogeneous geometry with an epidermal melanin content of 1%. The voxel size was set to 50 μm × 50 μm. The thickness of epidermis was set to 50 μm. Simulations were run with 50 000 000 photons per wavelength.
For the vessels, the voxel size was decreased to 10 μm × 10 μm in order to reduce the influence of voxelization on the results. The depth and melanin content of the epidermis were set to 60 μm and 3%. Reflectance values were meaned over 5 × 5 pixels in order to yield a pixel size corresponding to the pixel size in the mole simulation. Scattering and absorption in the blood were modeled using properties obtained from Friebel et al. [29]. The refraction index was set to the values reported by Li et al. [30],
Simulations were run with 100 000 000 photons per wavelength. Geometries are shown in Fig. 6.3. Results and discussion
The DWT algorithm rearranges the information in the hyperspectral images and emphasizes features of the analyzed spectra appearing at different scales and wavelengths. By visual inspection, useful images can be recognized, highlighting some desired parameters or structures. As a pure mathematical method, however, the bridge between the derived results and tissue properties is not readily apparent.
In this work, simulations have been used to investigate the link between tissue properties and DWT coefficients. Simulated reflectance spectra with known parameters were decomposed, and the coefficients were correlated to the tissue parameters. This was used to identify which components can be utilized to map levels of blood, melanin and the oxygen saturation. These coefficients were tested on measured hyperspectral data, simulated data and tissue phantoms.
3.1. Simulated reflectance spectra
Diffuse reflectance spectra simulated using Monte Carlo and scattering alternative I are shown in Fig. 7. The effect of varying the dermal blood volume fraction (BVF), the melanin content, and oxygen saturation is shown in Figs. (a), (b) and (c), respectively.
Fig. 7 Diffuse reflectance spectra, simulated using a three-layer Monte Carlo model and scattering alternative I. The spectra in (a) represent tissues with different dermal blood content. The effects of varying melanin content in epidermis are shown in (b). In (c), the effect of varying oxygen saturation levels in dermal blood can be seen.
Variations in BVF generally result in a down-shift of the reflectance spectra. Changes in oxygenation, with a fixed BVF, has an impact mainly on the spectral shapes. For this set of simulations, this is most pronounced for the oxyhemoglobin peaks at 542 and 576 nm. Changes in melanin content affects the steepness of the curves and decreases the visibility of the spectral features of hemoglobin in addition to reducing the overall reflectance. Melanin is known to show cross talk with deoxy-hemoglobin absorption, and these chromophores might be difficult to distinguish, especially in the red range of the spectrum [31].
The same spectra were also generated using scattering alternative II. Due to redundancy in the spectral properties displayed by these additional spectra, only one of these is compared to one of the spectra simulated using scattering alternative I in Fig. 8. The reflectance level is shifted, but the same features are still present.
Fig. 8 Comparison between diffuse reflectance spectra simulated using scattering alternative I and II. Oxygen saturation was 99%.
The relation between simulated and measured spectra is discussed in section 3.4.
3.1.1. Correlation to wavelet components
In table 2, the selected DWT coefficients have been listed along with correlation coefficients and p-values for all three input parameters and scattering levels. The corresponding approximate center wavelength of each coefficient is stated in the same table. For coefficient n, at detail level m, the center wavelength is calculated as the mean of the nth out of m equal parts of the full spectrum. The correlation results for scattering alternative II show the same trends as scattering alternative I, though with a slightly lower correlation for the oxygenation and higher p-values.
Table 2. Correlation between input parameters and a few selected DWT coefficients. Correlation coefficients R and probabilities p have been calculated for varying input parameters and output DWT coefficients of the simulated spectra as described in section 2.6. λ is the approximate center wavelength for the DWT coefficient.
3.2. Simulated hyperspectral images
3D Monte Carlo simulations were run for a set of skin geometries. DWT coefficients (approximation 9, detail 2.17 and a combination of approximation 4 and 5) are shown in Fig. 9.
Fig. 9 DWT coefficients obtained from the 3D Monte Carlo simulations. Results obtained from the mole are shown in figures 9(a) to 9(c). Results obtained from the deep vessels are shown in Figs. 9(d) to 9(f). The results from the vessels going from the superficial layers to the deep layers are shown in Figs. 9(g) to 9(i). See Fig. 6 for the vessel geometry. Note that some of the coefficients have different intensity ranges in order to enhance features.
The combination of approximation 4 and 5 is affected mainly by melanin. This combined coefficient is featureless unless melanin is present. However, it seems to be affected by super-ficial deoxygenated vessels, as seen in Fig. 9(i). Approximation 9 is sensitive to the presence of vessels. Detail 2.17 shows a high sensitivity to low oxygenation in superficial vessels. The detail shows a larger footprint than the actual vessel. This might be an artifact caused by scattering. This detail is also highly affected by photon noise. Detail 2.17 is negatively affected by superficial blood content as the coefficient is lower for the superficial vessel than the surrounding tissue although the oxygenation in the vessel is lower than the surrounding tissue (0.7 and 0.8, respectively).
3.3. Phantom data
DWT coefficients extracted from phantom data are shown in Fig. 10. The images were cropped so that only results from the submerged parts of the capillary tubes are shown. Approximation 9 in Fig. 10(a) has a different scaling than approximation 9 in Fig. 12 for visualization purposes. It should be noted that the simple phantom does not completely represent the properties of skin. No epidermis, and thus no pigmentation is included. The blood was contained in glass tubes with a different refractive index, which might affect the results. Despite this, the phantom does illustrate the effects of depth and oxygenation quite well. Detail 2.17 in Fig. 10(b) displays a higher sensitivity to oxyhemoglobin than deoxyhemoglobin. The values of the detail are lower above the capillary tubes than in the surrounding intralipid, even though the blood volume fraction is higher. The combination of approximation 5 and 4 shows some sensitivy to deoxyhemoglobin, as it did for the 3D Monte Carlo simulations, this illustrated the known crosstalk between de-oxyhemoglobin and melanin [31]. The high-contrast line around line 780 and 230 is due to an artifact from the imaging setup.
3.4. In vivo image data
The wavelet methods were tested on in vivo image data where oxygenation and blood volume were manipulated using cuff occlusion. A collection of measured spectra are shown in Fig. 11. The spectra were extracted from three qualitatively different regions of the imaged limb (normal skin, mole, blood vessel). All shown spectra are averaged over 10 by 10 pixels.
Oxygenation and BVF results obtained from normal skin using the inverse diffusion model are shown in table 3.
Table 3. Calculated oxygen saturation levels and blood volume fractions for plain skin before, during and after occlusion.
These results show that the oxygenation decreases during occlusion and rises after. The oxygenation has not yet reached the pre-occlusion value 3 minutes after the cuff was released. The BVF show an increase during occlusion, followed by a decrease after releasing the strap around the arm. The BVF is not normalized completely after 3 minutes. This is also evidenced by inspection of the skin reflectance spectra. It is at its lowest during occlusion and slightly increased after occlusion, indicating a higher blood volume fraction. In the spectra, this can be seen as a reduced reflectance for all wavelengths. However, the effect is more pronounced at shorter wavelengths indicating a pooling of blood in the smaller vessels due to the occlusion. The change in oxygenation during occlusion can be seen as a transition from the double minima at 542 nm and 576 nm caused by oxyhemoglobin absorption into a broader minimum at 550 nm caused by deoxyhemoglobin. However, the changes seen at these wavelengths will reflect changes in the oxygenation in the superficial vessels due to a short penetration depth at these wavelengths. At 760 nm the penetration depth is longer. The minimum at 760 nm due to de-oxyhemoglobin becomes more prominent during occlusion, showing reduced oxygenation also in deeper layers. This agrees well with the calculated values.
Yu et al. [32] performed a similar occlusion experiment. Here, the limb was subjected to arterial cuff occlusion for 3 minutes. The oxygenation was shown to decrease during occlusion, as was the measured blood flow. The oxygenation was shown to increase after the release, however, before reverting to the original oxygenation. The same is seen in a similar study [33] for a longer occlusion period (5 minutes). Another study [34] shows a higher blood volume fraction and lower oxygenation during the occlusion, and a lower blood volume fraction and higher oxygenation after.
The data in our study show a slightly different trend, as the oxygenation and BVF has not completely returned to the starting value after 3 minutes. This might be due to individual differences, differences in the way the occlusion is imposed, and in the measurement configuration and geometry. However, the effects of the occlusion imposed in this study can be interpreted based on the measured reflectance spectra, which then serves as a basis for the interpretation of the DWT coefficients.
All 9 measured spectra generally indicate the presence of lipids (at 930 nm) and water (at 970 nm). These features are not present or less pronounced in the simulated spectra in Fig. 7. Water was not explicitly included in the simulations, which explains the lack of the water peak. The subject had a very fair skin with an observed reflectance of almost 60% for the longer wavelengths. This indicates a very low melanin content and/or a high scattering coefficient. The assumed skin thickness in the simulations might also affect the visibility of the lipid peak, as it is assumed to be due to subcutaneous fat.
3.5. Blood vessel visualization
Approximation coefficient 9 was identified as an image with high contrast between blood vessels and surrounding tissue through visual inspection of the DWT output. See Fig. 12. It shows very strong negative correlation to all tested parameters (blood, melanin and oxygenation), as can be seen in table 2. Regardless of the cross-talk in this coefficient, it provides images with good contrast for blood vessels that can be used for further processing to extract vessel structures. This information is of interest for example in diagnosing arthritis in finger joints.
The vessels become more pronounced during and after the cuff occlusion, in accordance with the increased blood volume evidenced by the reflectance spectra and inverse modeling results.
Detail 2, coefficient 17 (detail 2.17) is shown in Fig. 13. This coefficient is more complicated to interpret and shows some conflicting trends that have to be further investigated. The image in Fig. 13 shows good contrast between vessels and normal tissue. The visualized vessels are different from the most apparent vessels seen in Fig. 12. The central wavelength for approximation 9 is 960 nm, while the central wavelength for detail 2.17 is 627 nm. Due to the higher absorption and scattering in this wavelength range, the penetration depth is shorter. Thus, detail 2.17 would likely be sensitive to more superficial changes in the tissue, than approximation 9. This is promising and in accordance with the 3D Monte Carlo results.
The 3D Monte Carlo confirm the differences in penetration depth. Both detail 2.17 and approximation 9 show sensitivity to superficial vessels in Fig. 9(h) and 9(g). Detail 2.17 is more clearly affected. Only approximation 9 shows some features for the deep vessels in Fig. 9(d). Here, detail 2.17 is almost featureless (see Fig. 9(e)). The phantom data in Fig. 10 show the same trend. It is thus expected that the vessels seen with strong contrast in Fig.12 are deep, while the vessels shown in Fig. 13 are more superficial.
The 3D Monte Carlo simulations and phantom data show that oxygenated blood leads to lower coefficient values as compared to deoxygenated blood. This is in accordance with the trend predicted by the 1D Monte Carlo simulations. Table 2 shows that detail 2.17 has a statistically significant negative correlation to oxygenation (R = −1.0 for scattering alternative I, p < 0.02).
However, these results show a discrepancy when compared to the 1D Monte Carlo simulations. During occlusion, where oxygenation decreases, the detail 2.17 coefficient decreases (see Fig. 13(b)), although the negative correlation to oxygenation indicates that the coefficient should increase. There is no significant correlation with blood content in this case.
The 3D Monte Carlo simulations show the same trend. The superficial parts of the vessel in Fig. 9(h), and the deoxygenated vein (oxy = 0.70) exhibit lower coefficient values than the rest of the tissue (oxy = 0.80). The same can be seen for the phantom data. Both the blood filled tubes have lower coefficient values than the rest of the phantom. The observed trends will thus be further investigated using more advanced statistical tools than correlation analysis to identify the reason for the observed discrepancy.
It might be that detail 2.17 is sensitive to a very high blood content superficially or throughout the skin, leading to decreased coefficient values. If this is the case, it will not be in accordance with the 1D Monte Carlo simulations. This might be due to the chosen simulation geometry and parameter variation. A weakness in the 1D Monte Carlo simulations in this study is that the variations of one specific input parameter was not subjected to the co-variation of the other parameters. This will be investigated further in future work.
Detail 2.17 and approximation 9 were visualized using a fixed intensity range in order to compare the variation throughout the occlusion experiment. However, for the purposes of vessel visualization, it will be more appropriate to scale each coefficient image dynamically, as the intensity level can be subject to variations in melanin and scattering.
3.6. Melanin visualization
While individual DWT coefficients can be used to map the properties of the imaged tissue, they also display a high redundancy, as shown in table 2. This is to be expected as the wavelet method creates 199 new images out of the original 160 bands. Many of the observed features seen along the spectral axis are due to variations in the same chromophores, and the coefficients will thus be highly correlated. Differences in penetration depth across the included spectral range can make it viable to combine these coefficients into new coefficients. One example is the enhancement of differences in melanin.
The negative correlation of oxyhemoglobin for approximation 5 and the positive correlation for approximation 4, together with the consistently negative correlation for melanin, suggests that these could be combined to a oxygenation-independent measure of the melanin content. This assumption has been investigated using measured data or 3D Monte Carlo simulations. As can be seen by studying the moles in Figs. 14(a) and 14(b), the correlation to blood is somewhat different in the two coefficients. High levels of melanin are indicated by lower values in coefficient 4 than in coefficient 5. Inspection of the difference between the two images indicate that the correlation to melanin can be emphasized using these coefficients. Due to the similar correlation to blood content and oxygenation in these two coefficients, effects from blood can be suppressed by subtracting one from the other. In Fig. 14(c), coefficient 4 has been subtracted from coefficient 5, highlighting the moles containing melanin.
Fig. 14 DWT approximation coefficients 4 and 5, and the difference between them, which emphasizes melanin content. Image acquired before occlusion.
The combination of DWT coefficients was tested on simulated data with known geometries. Results for 3D Monte Carlo simulations have already been shown, suggesting that blood related properties are suppressed while melanin related properties are enhanced. In Fig. 15, the combined DWT coefficients was plotted as a function of input parameter for different input parameters. This shows that the new ”coefficient” maps the melanin content dynamically. This suggests that the combination could be used for melanin mapping, e.g. as a melanin index. Changes in scattering levels shows that the absolute measure obtained from this combination of coefficients could be subject to changes in scattering, though it can serve as an relative measure within skin regions having the same levels of scattering.
Fig. 15 DWT approximation coefficient 05 minus coefficient 04, correlation to melanin, blood and oxygenation, as applied on simulated 1D Monte Carlo spectra.
In Fig. 16(a), this melanin mapping has been applied to a synthesized image constructed from the simulated spectra to illustrate how varying melanin levels are clearly separated, while differences in blood content and oxygenation are suppressed. Figure 16(c) is a scaled version of 14(c), where a blood vessel turns up with negative values, presumably due to high oxygenation and/or blood content.
It is interesting to compare this melanin map to a previously established melanin index, introduced by Dawson et al [35]. Applying Dawson’s melanin index to the synthesized image, seems to map melanin well, but also falsely indicates high blood fractions as high melanin content. For comparison, Dawson’s index has also been calculated for the measured data, as can be seen in Fig. 16(d). The cross-talk with hemoglobin is visible here as well. Based on this simple comparison it seems like wavelets can be used to reduce the crosstalk between deoxyhemoglobin and melanin.
3.7. Noise
The wavelet decomposition methods have been tested on noise-reduced reflectance spectra. These spectra were obtained either through a high photon count in the Monte Carlo simulations or by employing spectral noise removal on the results. It is likely that the method will be affected by noise in the data. As an example, detail 2.17 shows sensitivity to noise in the Monte Carlo simulations in Figs. 9(e) and 9(h), while the approximations remain largely unaffected.
All data in this study were subject to noise removal, and noise removal is thus considered to be a requirement for consistent results.
3.8. Summary and future work
The above examples show the potential of applying the wavelet transform to hyperspectral images of tissue to extract features.
Wavelets can be used to extract information from hyperspectral images and emphasize features at different scales and wavelengths. However, to be able to fully utilize the potential of the method it is important to relate the output from this mathematical decomposition of spectral data to physiological parameters. The analysis performed in this study suggests that selected co-efficients can be combined to yield estimates of melanin content, blood volume or oxygenation. Other multispectral or spectral methods can achieve similar results [21], but wavelet methods would have an advantage due to simplicity and computational efficiency. Although this paper does not show absolute chromophore concentration values, it shows the relation between the parameters. These relations will be further investigated in future work, with the aim of more quantitative results.
This study is intended as proof of concept, and the main goal was to show a relation between the mathematically derived coefficients and the true physiological parameters. At the present this method can be used for treatment monitoring and follow-up over time, where relative changes in the physiological parameters are the most important. The method will be applied to extract information about the vasculature and blood content in patients with arthritis [36].
Future work will include building a larger database of simulated skin spectra in order to elucidate the more subtle relationships between the true physiological state and the derived co-efficients, and more advanced statistical evaluation of the method. There is obviously a relation between the depth of the observed structures and the output from the wavelet algorithm. This will be further investigated. Other coefficients or other combinations of coefficients than the ones presented in this paper can also prove to be more efficient. There is still a large number of DWT components that have not been fully investigated.
Other authors have reported the use of wavelets for explicit oxygenation determination from tissue reflectance spectra [37] and several wavelet based feature extraction schemes have been suggested for classification in remote sensing [38,39]. The present study aims to generalize the application of wavelet analysis for feature extraction and visualization in hyperspectral imaging of tissue. The presented results can potentially be used for classification and diagnostics in clinical hyperspectral imaging.
4. Conclusion
DWT analysis has been used to extract features of interest from diffuse reflectance spectra, collected with a hyperspectral imaging system. The output from the wavelet algorithm was correlated with blood content, blood oxygenation and melanin content. Through examination of coefficient images and correlation analysis of simulated spectra, melanin maps and high contrast images suitable for blood vessel visualization was constructed. The DWT technique is general but can be optimized for specific features and systems, which makes it a promising tool for visualization and characterization of tissue.
Acknowledgments
The project has been funded by NTNU and through the Iacobus project, supported by the European Commission’s 7th RTD Framework Programme Collaborative Project No. 305760.
Acquisition of in vivo hyperspectral data was done by Lukasz Paluchowski. The authors would also like to acknowledge the work of Tore Berg and Tore Landsem in the mechanics workshop at the Department of Electronics and Telecommunications, NTNU, for building the bed-side scanner which was used in the clinical imaging.
A special thanks to professor Yurii Lyubarskii at the Department of Mathematical Sciences, NTNU, for sharing his insight into wavelet theory.
References and links
1. G. N. Stamatas, B. Z. Zmudzka, N. Kollias, and J. Z. Beer, “In vivo measurement of skin erythema and pigmentation: New means of implementation of diffuse reflectance spectroscopy with a commercial instrument,” Brit. J. Dermatol. 159, 683–690 (2008). [CrossRef]
2. D. Yudovsky, A. Nouvong, and L. Pilon, “Hyperspectral imaging in diabetic foot wound care,” J. Diabetes Sci. Technol. 4, 1099–1113 (2010). [CrossRef] [PubMed]
3. M. Denstedt, B. S. Pukstad, L. Paluchowski, J. E. Hernandez-Palacios, and L. L. Randeberg, “Hyperspectral imaging as a diagnostic tool for chronic skin ulcers,” Proc. SPIE 8565, 85650N (2013). [CrossRef]
4. K. Rajpoot and N. M. Rajpoot, “Hyperspectral colon tissue cell classification,” in “SPIE Medical Imaging,” (2010).
5. G. Camps-Valls and L. Bruzzone, “Kernel-based methods for hyperspectral image classification,” IEEE T. Geosci. Remote 43, 1351–1362 (2005). [CrossRef]
6. M. Fauvel, J. A. Benediktsson, J. Chanussot, and J. R. Sveinsson, “Spectral and spatial classification of hyper-spectral data using svms and morphological profiles,” IEEE Trans. Geosci. Remote 46, 3804–3814 (2008). [CrossRef]
7. S. Prigent, X. Descombes, D. Zugaj, and J. Zerubia, “Spectral analysis and unsupervised svm classification for skin hyper-pigmentation classification,” in “2nd Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing, WHISPERS 2010 - Workshop Program,” (2010).
8. L. M. Bruce, C. H. Koger, and J. Li, “Dimensionality reduction of hyperspectral data using discrete wavelet transform feature extraction,” IEEE Trans. Geosci. Remote 40, 2331–2338 (2002). [CrossRef]
9. T. Skauli, P. E. Goa, I. Baarstad, and T. Løke, “A compact combined hyperspectral and polarimetric imager,” Proc. SPIE 6395, 639505 (2006). [CrossRef]
10. A. A. Green, M. Berman, P. Switzer, and M. Craig, “Transformation for ordering multispectral data in terms of image quality with implications for noise removal,” IEEE Trans. Geosci. Remote 26, 65–74 (1988). [CrossRef]
11. L. O. Svaasand, L. T. Norvang, E. J. Fiskerstrand, E. K. S. Stopps, M. W. Berns, and J. S. Nelson, “Tissue parameters determining the visual appearance of normal skin and port-wine stains,” Laser Med. Sci. 10, 55–65 (1995). [CrossRef]
12. A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “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). [CrossRef]
13. R. R. Anderson and J. A. Parrish, The Science of Photomedicine (Plenum Press, 1982), chap. 6: Optical properties of human skin, pp. 147–194. [CrossRef]
14. T. Spott, L. O. Svaasand, R. E. Anderson, and P. F. Schmedling, “Application of optical diffusion theory to transcutaneous bilirubinometry,” Proc. SPIE 3195, 234–245 (1997). [CrossRef]
15. D. SegelsteinThe complex refractive index of water, Master'sthesis, University of Missouri-Kansas City (1981).
16. R. L. P. van Veen, H. J. C. M. Sterenborg, A. Pifferi, A. Torricelli, and R. Cubeddu, “Determination of vis- nir absorptioncoefficients of mammalian fat, with time- and spatially resolved diffusereflectance and transmission spectroscopy,”(OSA, 2004).
17. R. Anderson and J. Parrish, “The optics of human skin,” J. Invest. Dermatol. 77, 13–19 (1981). [CrossRef] [PubMed]
18. A. Boggess and F. J. Narcowich, A First Course in Wavelets with Fourier Analysis, 2nd ed. (Wiley, 2009).
19. S. G. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Trans. Pattern Anal. 11, 674–693 (1989). [CrossRef]
20. I. Daubechies, Ten Lectures on Wavelets, no. 61 in CBMS/NSF Series in Applied Math. (1992). [CrossRef]
21. A. Bjorgan, M. Milanic, and L. L. Randeberg, “Estimation of skin optical parameters for real-time hyperspectral imaging applications,” J. Biomed. Opt. 19, 066003 (2014). [CrossRef] [PubMed]
22. L. L. Randeberg, E. L. P. Larsen, and L. O. Svaasand, “Hyperspectral imaging of blood perfusion and chromophore distribution in skin,” Proc. SPIE 7161, 71610C (2009). [CrossRef]
23. E. Alerstam, W. C. Y. Lo, T. D. Han, J. Rose, S. Andersson-Engels, and L. Lilge, “Next-generation acceleration and code optimization for light transport in turbid media using gpus,” Biomed. Opt. Express 1, 658–675 (2010). [CrossRef]
24. E. Salomatina, B. Jian, J. Novak, and A. N. Yaroslavsky, “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11, 064026 (2006). [CrossRef]
25. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58, R37–R59 (2013). [CrossRef] [PubMed]
26. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: A review,” J. Innov. Opt. Health Sci. 4, 9–38 (2011). [CrossRef]
27. P. Naglic, L. Vidovic, M. Milanic, L. L. Randeberg, and B. Majaron, “Applicability of diffusion approximation in analysis of diffuse reflectance spectra from healthy human skin,” Proc. SPIE 9032, 90320N (2013). [CrossRef]
28. M. Milanic and B. Majaron, “Three-dimensional monte carlo model of pulsed-laser treatment of cutaneous vascular lesions,” J. Biomed. Opt. 16, 128002 (2011). [CrossRef] [PubMed]
29. M. Friebel, A. Roggan, G. Mueller, and M. Meinke, “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, 034021 (2006). [CrossRef]
30. H. Li, L. Lin, and S. Xie, “Refractive index of human whole blood with different types in the visible and near-infrared ranges,” Proc. SPIE3914 (2000). [CrossRef]
31. G. N. Stamatas and N. Kollias, “Blood stasis contributions to the perception of skin pigmentation,” J. Biomed. Opt. 9, 315–322 (2004). [CrossRef] [PubMed]
32. G. Yu, T. Durduran, G. Lech, C. Zhou, B. Chance, E. R. Mohler, and A. G. Yodh, “Time-dependent blood flow and oxygenation in human skeletal muscles measured with noninvasive near-infrared diffuse optical spectroscopies,” J. Biomed. Opt. 10, 024027 (2005). [CrossRef] [PubMed]
33. R. A. De Blasi, S. Fantini, M. A. Franceschini, M. Ferrari, and E. Gratton, “Cerebral and muscle oxygen saturation measurement by frequency-domain near-infra-red spectrometer,” Med. Biol. Eng. Comput. 33, 228–230 (1995). [CrossRef] [PubMed]
34. G. Zonios, J. Bykowski, and N. Kollias, “Skin melanin, hemoglobin, and light scattering properties can be quantitatively assessed in vivo using diffuse reflectance spectroscopy,” J. Invest. Dermatol. 117, 1452–1457 (2001). [CrossRef]
35. J. B. Dawson, D. J. Barker, D. J. Ellis, J. A. Cotterill, E. Grassam, G. W. Fisher, and J. W. Feather, “A theoretical and experimental study of light absorption and scattering by in vivo skin,” Phys. Med. Biol. 25, 695–709 (1980). [CrossRef] [PubMed]
36. IACOBUS project web page, http://www.iacobus-fp7.eu.
37. P. C. Chen and W. C. Lin, “Spectral-profile-based algorithm for hemoglobin oxygen saturation determination from diffuse reflectance spectra,” Biomed. Opt. Express 2, 1082–1096 (2011). [CrossRef] [PubMed]
38. T. R. West, S. Prasad, and L. M. Bruce, “Wavelet packet tree pruning metrics for hyperspectral feature extraction,” in Geoscience and Remote Sensing Symposium, 2008. IGARSS 2008 (IEEE, 2008), pp. II–946–II–949.
39. S. Prasad, W. Li, J. E. Fowler, and L. M. Bruce, “Information fusion in the redundant-wavelet-transform domain for noise-robust hyperspectral classification,” IEEE. Trans. Geosci. Remote 50, 3474–3486 (2012). [CrossRef]
