Calibration and segmentation of skin areas in hyperspectral imaging for the needs of dermatology

Hyperspectral imaging is now one of the most developed methods of visible light imaging [1, 2]. It enables to acquire data of an object in any spectral range set in the camera. This method is entirely non-contact and non-invasive, and measurements can be carried out remotely. With these advantages, hyperspectral cameras are used in many fields of technology and medicine [3–7]. In particular, they are used in dermatology. The issue of dermatological research concerns spectral skin analysis in almost all cases. This analysis is usually carried out on the basis of data recorded by a camera and saved in dat or raw format, or another one dependent on the type of camera used [8, 9]. These data are further analysed in the software provided directly by the camera manufacturer or, less often, in another software. The latter option is rarely used as there is a need to meet the following requirements:

A large amount of data prevents simple applications from analysing such large images - a typical file size is several gigabytes, as mentioned earlier. The data in the file are saved in 32 or 16-bit format depending on its type, dat or raw. Consequently, this leads to differences in the file volume and the quality of the resultant image. Regardless of the format of data record, a sequence of hyperspectral images requires calibration. Due to the different brightness of lighting (from the definition of linear spectral characteristics), calibration of images in the range from black (0% brightness) to white (100% brightness) is necessary. This goal is achieved in different ways. In the simplest case, it is a white stripe, pattern, taken as 100% brightness, that is placed next to the imaged skin. Subsequent images are calibrated with respect to it. The calibration stage ends with suitable linear adjustment of brightness levels. Placing a white stripe pattern at the top of the image is sometimes used and recommended by Specim Company and LOT-Quantum Design Company. Automatic its analysis here is not done in a dedicated program ENVI. Therefore, it is not possible implementation of automatic calibration.

The second major area of problems related to hyperspectral imaging is image analysis. On the one hand, typical image analysis and processing in visible light can be applied. On the other hand, the acquisition of a separate image for each length of the spectrum gives much wider possibilities of analysis. Typical tasks, most commonly used in dermatology, include segmentation of skin areas. In general, morphological methods [10–15], statistical methods [16–18] and algorithms profiled to selected applications are known from hyperspectral imaging. Among the morphological methods, there are classical approaches [1, 13, 14] and those profiled to the analysis of an image sequence [12, 15]. In the statistical methods, there dominates texture analysis [17–21] used as a set of features for classification and recognition. Profiled algorithms have been applied so far to face recognition [22], analysis of skin areas [23], and others [24–27]. These methods are mainly related to segmentation of specific objects [28]. On the basis of segmented objects, their morphometric measurements or their texture analysis are carried out [18]. However, for skin areas, an important feature is their location in the image, which is vital, inter alia, in the assessment of efficacy and safety of treatments in aesthetic medicine [29–31]. A description of this type of automatic method preceded by automatic calibration is presented below.

The study used images of human healthy skin obtained with the Specim camera PFD-V10E. The measured body areas were lit using a typical lamp with flat spectral characteristics in the required range (based on HgAr emission for the VNIR spectral range). The images were obtained retrospectively during routine medical (dermatological) examinations, performed in accordance with the Declaration of Helsinki. In connection with the described algorithm, no research or experiments were carried out on humans. The acquired data were anonymised and stored in the output format, source dat (ENVI File) and raw. The frequency λ of the data obtained ranged from 397 to 1030 nm. Each image was recorded every 0.79 nm, which in total gave 800 2D images for each patient. The resolution M × N (number of rows and columns) of each image for the selected frequency was varied depending on the scan area. It was usually M × N = 899 × 1312 pixels. The number of rows N was changed most often and was in the range N∈(15,899). Regardless of the distance of the camera from the object and selected focusing parameters, one pixel covered a square area in the range of 130 μm × 130 μm. A total of 36’000 2D images were obtained for 45 hyperspectral images. These images were subjected to further analysis.

Hyperspectral image analysis method is associated with three stages:

Details of these three stages are described below.

Pre-processing

Image pre-processing concerns correct reading and interpretation of the data recorded by the camera PFD-V10E in dat and raw format. This camera records information for each line (each row) N and at the same time registers a full spectral range. In this case, λ∈(397, 1030) nm is equivalent to the adopted spectral distance with the registration of 800 lines. This process is shown in Figure 1. Depending on the type of data, raw or dat, each pixel is recorded at 32 or 16 bits of data. The exact number of rows and columns is stored in a file with the extension hdr that contains all the typical header information. This is information relating to a particular frequency of the spectral range, type of data storage, sensor type, and others. For the analysed data, the image resolution M × N was varied in the range from 15 × 1312 to 899 × 1312 pixels. The range of changes was strictly dependent on the scan area. This area was limited mainly due to image acquisition time of 12 seconds for the registration of the spectrum at the maximum resolution and in the full range. A dynamic error related to the possible displacement of the scan area during measurements was minimized by mechanical stops and skin area orientation ensuring patient’s comfort.

The images L _GRAY (m,n,k), where m-row, n-column and k–the next wavelength λ (k∈(1,K)), read from the files with the extension dat or raw were further filtered. For each image in the sequence, a median filter with a mask h sized M _h  × N _h  = 3 × 3 pixels was used. The mask size was dependent on the amount of pollution and the level of noise. In the case of recorded images, the noise and artefacts did not exceed the size of 2 pixels per one cluster. For this reason, a sufficient filter mask size was 3 × 3 pixels. In this way, the noise-free image L _M (m,n,k) was subjected to calibration.

Calibration

The acquired images L _M (m,n,k) are not calibrated. Calibration involves referring each pixel of the registered skin area to the white pattern [32–34]. For the registered cases, the pattern was a white stripe placed at the top- Figure 2. Automatic detection of the pattern position was implemented in the proposed algorithm. It concerned recognition of one of the pattern contours using information about the brightness gradient of adjacent pixels for each column, i.e.:

$L_G\left(m,n,k\right)=\left\{\begin{array}{ccc}\hfill 0\hfill & \hfill \mathit{if}\hfill & \hfill \left(L_M\left(m,n,k\right)-L_M\left(m+1,n,k\right)\right)>p_r\hfill \\ \hfill 1\hfill & \hfill \mathit{other}\hfill \end{array}\right.$

(1)

for m∈(1,N-1) where p _r is a binarization threshold determined automatically according to Otsu’s formula [35].

The pattern in the form of a white stripe was 20 × 400 mm, which was equivalent to the number of rows m _w  = 80 ± 5 for its set distance from the camera lens. The value of ±5 pixels is associated with a possible image shift or rotation. The number of columns was covered by the pattern in its entirety. Therefore, the searched pattern boundary contour was designated as:

$W_I\left(n,k\right)=\left\{\begin{array}{ccc}\hfill m\hfill & \hfill \mathit{if}\hfill & \hfill \left(L_G\left(m,n,k\right)=1\right)\wedge \left(m<m_w\right)\hfill \\ \hfill 0\hfill & \hfill \mathit{other}\hfill \end{array}\right.$

(2)

On this basis, average brightness for each column is calculated, i.e.:

$L_w\left(n,k\right)=\frac{1}{W_I\left(n,k\right)}·{\displaystyle \sum _{m=1}^{W_I\left(n,k\right)}}{L}_M\left(m,n,k\right)$

(3)

Examples of graphs of L _w (n,k) for k = 400, 401 and 402 are shown in Figure 3. The image Figure 3 a) and its zoom Figure 3 b) show the differences in average brightness values. The image must be calibrated with respect to these changes. For each value k, calibration must be performed independently. Calibration of individual images is carried out as:

$L_k\left(m,n,k\right)=\frac{L_M\left(m,n,k\right)-{\min }_{m,n}\left(L_M\left(m,n,k\right)\right)}{L_w\left(n,k\right)\underset{m,n}{-\min }\left(L_M\left(m,n,k\right)\right)}$

(4)

for: $\left(L_w\left(n,k\right)\underset{m,n}{-\min }\left(L_M\left(m,n,k\right)\right)\right)\ne 0$

In the case of pixels which exceed the value “1”, adjustment is necessary:

$L_K\left(m,n,k\right)=\left\{\begin{array}{ccc}\hfill L_K\left(m,n,k\right)\hfill & \hfill \mathit{if}\hfill & \hfill L_K\left(m,n,k\right)\le 1\hfill \\ \hfill 1\hfill & \hfill \mathit{other}\hfill \end{array}\right.$

(5)

The image L _K (m,n,k) having brightness values in the range from 0 to 1 is further subjected to the next processing steps.

Image processing

The input image L _K (m,n,k) after calibration is the basis for the segmentation process. For this purpose, a sample diagram of brightness changes in a sample ROI was made for the human skin which mainly consists of water, melanin and haemoglobin. The results for each k-th image (at different wavelengths) are shown in Figure 4 a), i.e.:

$L_{\mathit{ROI}}^{\mathit{MED}}\left(k\right)=\frac{1}{M_{\mathit{ROI}}·N_{\mathit{ROI}}}{\displaystyle \sum _{m,n\in \mathit{ROI}}}{L}_K\left(m,n,k\right)$

(6)

$L_{\mathit{ROI}}^{\mathit{MIN}}\left(k\right)=\underset{m,n\in \mathit{ROI}}{\min }{L}_K\left(m,n,k\right)$

(7)

$L_{\mathit{ROI}}^{\mathit{MAX}}\left(k\right)=\underset{m,n\in \mathit{ROI}}{\max }{L}_K\left(m,n,k\right)$

(8)

The ROI was associated with the hand area shown in Figure 2 and included the range M _ROI  × N _ROI  = 150 × 150 pixels. In this range, the values $L_{\mathit{ROI}}^{\mathit{MED}}\left(k\right),L_{\mathit{ROI}}^{\mathit{MIN}}\left(k\right),L_{\mathit{ROI}}^{\mathit{MAX}}\left(k\right)$, which are the mean, minimum and maximum values of brightness changes in the ROI respectively, were calculated. The results obtained shown in Figure 4 a) are also dependent on the individual variability of patients and the method of lighting and setting the camera angle relative to the patient. The influence of these elements on the result is revealed by the shift of curves shown in Figure 4 a) up or down, which increases or decreases the mean brightness value. Therefore, normalization performed for the entire image sequence with respect to changes for k-th images is necessary, i.e.:

$L_O\left(m,n,k\right)=\frac{L_K\left(m,n,k\right)-\underset{k}{\min }{L}_K\left(m,n,k\right)}{\underset{k}{\max }\left(L_K\left(m,n,k\right)-\underset{k}{\min }{L}_K\left(m,n,k\right)\right)}$

(9)

for $\left(L_K\left(m,n,k\right)-\underset{k}{\min }{L}_K\left(m,n,k\right)\right)\ne 0$

For the images L _O (m,n,k) modified in this way the obtained results of the mean, minimum and maximum values also change in the same sample ROI, i.e.: $L_{\mathit{ROI}}^{\mathit{MMED}}\left(k\right),L_{\mathit{ROI}}^{\mathit{MMIN}}\left(k\right),L_{\mathit{ROI}}^{\mathit{MMAX}}\left(k\right)$. The results obtained are shown in Figure 4b. The normalized images L _O (m,n,k) also enable automatic segmentation in accordance with the reference curve of melanin and haemoglobin content for each wavelength. The reference content of melanin and haemoglobin can be acquired from external sources, for example from literature data [36], or on the basis of the selected ROI. In the latter case, the result will be as follows - image L _D (m,n), i.e.:

$L_D\left(m,n\right)=\frac{1}{M·N}{\displaystyle \sum _{k=1}^K}\left|L_O\left(m,n,k\right)-L_{\mathit{ROI}}^{\mathit{MMED}}\left(k\right)\right|$

(10)

Therefore, the image L _D (m,n) contains information about the average error- Figure 5. It is calculated for individual pixels relative to the reference waveform $L_{\mathit{ROI}}^{\mathit{MMED}}\left(k\right)$. On this basis, binarization-based segmentation can be performed for the binarization threshold p _w designated manually or automatically (the afore-mentioned Otsu formula [35]). From a practical point of view, the effect of the binarization threshold selection (provided manually) on the segmentation results obtained is of interest. For this purpose, the impact of changes in the threshold p _w on the changes in the surface area A(p _w ) of the segmented object was investigated, i.e.:

$A\left(p_w\right)=\frac{1}{M·N}{\displaystyle \sum _{m=1}^M}{\displaystyle \sum _{n=1}^N}{L}_{\mathit{DB}}\left(m,n\right)$

(11)

where:

$L_{\mathit{DB}}\left(m,n\right)=\left\{\begin{array}{ccc}\hfill 1\hfill & \hfill \mathit{if}\hfill & \hfill L_D\left(m,n\right)<p_w\hfill \\ \hfill 0\hfill & \hfill \mathit{other}\hfill \end{array}\right.$

(12)

The results obtained are shown in Figure 6. The area optimal from the point of view of the segmentation results is marked in green. The term “optimal” refers to such a fragment of the curve A(p _w ) which is a flat area. For the threshold p _w  = 0.2 ± 0.9, the segmentation result is correct. This fact was proven when comparing it with the result obtained by an expert relying on manual marking (gold standard). The error is here defined as:

${\delta }_w\left(p_w\right)=\frac{A\left(p_w\right)-A_Z}{A_Z}$

(13)

where A _Z is the surface area resulting from the expert’s work.

The results obtained from equation (13) it was assumed that the results obtained by the expert are repetitive and do not differ from results obtained by other experts. Not in every case, however, it must be so. The exact considerations presented in [37–39].

The results of changes in the error δ _w for varying p _w (p _w ∈(0.1, 0.3)) are shown in Figure 7a). From the presented graph and for the case under consideration, the smallest error (δ _w  ≈ 0) is obtained for p _w  = 0.23. For automatic selection of the binarization threshold (Otsu’s formula), the error is 9%. Sample binarization results are shown in Figure 7b for the thresholds p _w ∈{0.1, 0.15. 0.2, 0.27}. Visual assessment gives the best results for p _w  = 0.23. However, it should be noted here that the adopted binarization threshold values are the acceptable standard deviation of the reference distribution of the skin spectrum from the measured pixels (formula (10)). In general, for any image containing the human skin, the following approaches are possible:

• automatic selection of the binarization threshold according to the Otsu’s formula – it enables to obtain a binary image of the object,

• manual selection of the binarization threshold p _w dependent on the acceptable tolerance of individual pixels of the image relative to the reference waveform $L_{\mathit{ROI}}^{\mathit{MED}}$ - 1% tolerance is p _w  = 0.01, 10% tolerance is p _w  = 0.1 respectively, etc.,

• automatic selection of the binarization threshold p _w depending on the location of the ‘flat’ area (Figure 6).

Depending on the desired end result, one of the above methods is selected by an operator. Figure 8 shows a sequence of images L _DB (m,n) for p _w ∈(0,1) changed with 0.2 step from the area of the hand, forearm, finger (thumb) and tattoo.

where i – subsequent image pixels, i.e. i∈(1, M∙N),$L_P^I\left(k\right)$ – reference values for the spectrum characteristic for the object I.

where: L _E(C)(m,n) – the resulting binary image after subjecting the image L _DB to conditional erosion,

L _D(C)(m,n) – the resulting binary image after subjecting the image L _DB to conditional dilation,

p _we – the constant that determines erosion effectiveness,

p _wd – the constant that determines dilation effectiveness,

p _mn (m,n) – the threshold dependent on the coordinates m, n,

s _re – the mean value of the analysed area for erosion,

s _rd – the mean value of the analysed area for dilation.

The values of p _nm ∈(0,1), whereas the values (1–p _we ) should be non-negative in the range from 0 to 2. The values (1–p _we ) and (1 + p _wd ) for p _we  = p _wd  = 0 are equal to 1, which means high intensity of conditional operations. For the other values of thresholds p _we and p _wd , for example, for p _we  = p _wd  = 1, there is a complete lack of effectiveness of the erosion operation and significant effectiveness of dilation. In the present case p _mn (m,n) = const and is independent of the location (p _mn  ≠ f(m,n)). The shape of the structural element SE in all these relationships was adopted as round sized 5 × 5 pixels due to the shape and size of the smallest objects corrected. These properties of conditional erosion and dilation enable to obtain effective correction of the surface of individual areas - in this case the tattoo, skin and background. This adjustment involves sequential implementation of conditional erosion and dilatation (in this case five times) for the images L _DB (m,n) and L _D (m,n). The binary image L _DBC (m,n), adjusted in this way, enables to obtain much better results - Table 2. The pixels not allocated to any of the areas were eliminated owing to the operations of conditional erosion and dilation. However, the skin area measurement error increased and in this case amounted to 12.9%. This error is closely dependent on the amount of segmented image areas, on the type and differences in the spectra for individual areas. Finally, the results obtained have a segmentation error of less than 13% compared to the work of the expert. Figure 10 shows segmentation errors for 45 hyperspectral image sequences, 36’000 2D images, analysed with the present method. The average value of the described method error δ _K fluctuates around 9% (the maximum value is 23%, minimum - 1%). These results are in the next section compared with the results obtained by other authors.

	Tattoo [pixels]	Skin [pixels]	Background [pixels]	Non-identified areas [pixels]	Measurement error of tattoo surface area [%]	Measurement error of skin surface area [%]	Measurement error of background surface area [%]
Algorithm	19066 (6.5%)	138365 (46.9%)	138114 (46.8%)	0 (0%)	1.8%	12.9%	1%

Hyperspectral imaging is the subject of many works on spectrum measurement and its analysis. In many studies, also simple functions related to image analysis and processing are used. Compared to conventional monochrome images, hyperspectral imaging contains much more information. This means that typical known methods of image analysis and processing must be generalized or profiled to their particular application. Among other things, the fact of having information on the complete spectrum in a specified range results in the possibility of performing more accurate segmentation than during typical segmentation. The analysis itself in the known reference works, however, usually refers to the analysis of the manually selected ROI, e.g. in [28]. There are many other areas of medicine which use similar manual or semi-automatic selection of the region of interest. These are, for example, the areas mentioned in [36] - vibrational hyperspectral imaging Filik J. [40], laparoscopic digital light processing – Olweny EO [41], blood stains at the crime scene – Edelman G. [42], prostate cancer detection – Akbari H. [43], histopathological examination of excised tissue - Vasefi F, diabetic foot ulcer - Yudovsky D [44], cancer detection - Akbari H [45], and others. In medicine, the use of hyperspectral imaging to assess the creation time of a bruise is also known – Stam B. [46]. The inaccuracy found is 2.3% for fresh bruises and 3 to 24% for bruises up to 3 days old. In conclusion, colour inhomogeneity of bruises can be used to determine their age. The experiment results in the work of Li Q [28] show that the hyperspectral based method has the potential to identify the spinal nerve more accurately than the traditional method as the new method contains both the spectral and spatial information on nerve sections. A strong resemblance of hyperspectral imaging to monochrome imaging means that typical methods such as morphological operations can be applied here. They are used in the classification of different types of artefacts visible in images [1, 13, 15] in SVM (support vector machines) [14], Gauss-Markov model [16] or wavelet analysis [18]. For example, in the work of Dicker et al. [31] spectral library was generated with 12 unique spectra that were used to classify specimens where sample preparation was varied. The work of Benediktsson J. et al. [11] presents results for the sequential use of morphological opening and closing for the increased size of the structural element. This methodology is similar to the use of conditional erosion and dilation, as in [13] and G. Rellier’s work [17]. These known methods of image analysis and processing do not enable fully automatic segmentation of the skin area, especially in conjunction with the prior automatic calibration.

The paper presents the method of calibration and segmentation of selected skin areas in hyperspectral imaging. The characteristics of the method described, whose block diagram is shown in Figure 11, include:

The presented method can be extended to:

Therefore, the discussed algorithm for image analysis and processing does not fully cover the issue. In terms of application, the techniques from spectral methods [47–49], analysis of microscopic images [50, 51] and others [52–59] can also be used. This type of analysis will be carried out by the authors in future studies in this area.