Deep Learning-Based Dermoscopic Image Classification System for Robust Skin Lesion Analysis

The DICS consists of the following modules: SLS, FE, FS, and IC modules. The following subsections explain these modules clearly and in detail. The idea of robust DICS is illustrated in Figure 1.

1.png

Figure 1. Flow of DICS for the diagnosis of skin cancer

2.1 SLS module

This is the first module of DICS where the exact skin cancer region is segmented using a k-means clustering approach on RGB (Red, Green, and Blue) colour images. K-means clustering is primarily a color-based segmentation approach that uses the colour values of the pixels in the image to cluster them into different groups. This means that pixels with similar color values are grouped, which can help to isolate the skin cancer region from the rest of the image. It does not require prior knowledge or training data to segment the image, i.e., an unsupervised learning approach. The algorithm involves the initialization of the cluster centers and continues to iterate until a convergence criterion is met. The output of the k-means clustering algorithm is a set of clusters, where each cluster corresponds to a different color region in the image.

Before segmentation, the image is smoothed by an averaging filter in the colour domain with a predefined window size which will remove noise and hair in the dermoscopic images. A larger sliding window of size 21×21 reduces the intensity variation between pixels, so noises and hairs will be removed completely. Let us consider a dermoscopic image DI. The pixel value of DI is represented by DI_ij, for i, j=1, ..., n, where n is the size of the image. The output O_ij from the linear filter of size (2m+1)×(2m+1) is denoted by.

$O_{i j}=\sum_{x=-m}^m \sum_{y=-m}^m w_{x y} D I_{i+x, j+y}$ for $i, j=(m+1), . .(n-m)$          (1)

where, w_xy for x, y=-m, ..., m is the weight of the filter. In this study, a linear filter (average filter) of size 21×21 (m=10) is used to smooth the image; hence, the weight is assigned equally to 1/21 for w_xy. The application of this filter to DI replaces each pixel with the average of pixels values in a 21×21 window centered on that pixel. For simplicity, let us consider a 3×3 filter. The averaging is carried out over this window size is defined by:

$\begin{aligned} O_{i j}= & \left(w_{-1,-1} D I_{i-1, j-1}+w_{-1,0} D I_{i-1, j}+w_{-1,1} D I_{i-1, j+1}+\right. \\ & w_{0,-1} D I_{i, j-1}+w_{0,0} D I_{i, j}+w_{0,1} D I_{i, j+1}+ \\ & \left.w_{1,-1} D I_{i+1, j-1}+w_{1,0} D I_{i+1, j}+w_{1,1} D I_{i+1, j+1}\right) / 9\end{aligned}$          (2)

Figure 2 shows the preprocessed dermoscopic images.

2.png

Figure 2. (a) Input dermoscopic images (Top Row images) (b) Smoothed images (Bottom row)

The RGB colour image is converted into L*a*b* colour space to quantify the visual differences. Figure 3 shows the filtered image's L,*a, and *b components. The conversion formula is as follows:

$L=11.6 h\left(\frac{Y}{Y_w}\right)-16$          (3)

$a=500\left[h\left(\frac{X}{X_w}\right)-h\left(\frac{Y}{Y_w}\right)\right]$          (4)

$b=200\left[h\left(\frac{Y}{Y_w}\right)-h\left(\frac{Z}{Z_w}\right)\right]$          (5)

where, $h(q)=\left\{\begin{array}{ll}q^{\frac{1}{3}} & q>0.008856 \\ 7.7879 q+\frac{16}{116} & q \leq 0.008856\end{array}\right\}$ .L, a, and b represent the lightness of the colour, red minus green and green minus blue. Also, X, Y, and Z correspond to R, G, and B channels, X_w, Y_w, and Z_w are reference white tri-stimulus values. This conversion enables obtaining all colour information in only two channels; a* and b*. Hence, k-means clustering is easily applied to cluster the colour information. It partitions the given data into k clusters. As the dermoscopic images contain skin lesions (abnormal region), unaffected skin region (normal skin area), and some background information (parts other than skin), the k value is set to 3. The steps to compute k-means clustering (k=3) for skin lesion segmentation are given below:

1. Randomly choose three initial clusters $\left(m_1^{(1)}, m_2^{(1)}\right.$ and $\left.m_3^{(1)}\right)$. The superscript identifies the initial cluster number.

2. The k-means clustering proceeds with assignment and update steps. These steps are iterated until it converges in the update step when the assignment no longer changes.

(a) Assignment step: In this step, each observation (z_p) is assigned to one of the clusters (C), which has the least squared Euclidean distance. Thus, each observation (z_p) is assigned to exactly one cluster $\left(C^{(t)}\right)$. It is defined by:

$C_i^{(t)}=\left\{z_p:\left\|z_p-m_i^{(t)}\right\|^2 \leq\left\|z_p-m_j^{(t)}\right\|^2 \forall j, 1 \leq j \leq 3\right\}$          (6)

(b) Update step: This step updates the means of the observations (z) in the new clusters ( $C_i^{(t)}$ ).

$m_i^{(t+1)}=\frac{1}{\left|C_i^{(t)}\right|} \sum_{z_j \in C_i^{(t)}} z_j$          (7)

3.png

Figure 3. RGB to L*a*b Conversion (a) RGB image (b) L channel (c) *a channel (d) *b channel (e) segmented lesion

The k-means clustering algorithm is widely applied in many medical image processing applications using gray-scale images. More information about the k-means clustering algorithm can be obtained from [12]. This study is applied to colour images with the help of colour information in a* and b* channels. Figure 3 (d) shows the segmented lesion.

The colour information in a* and b* channels are combined and fed to the k-means clustering to segment the lesion effectively. The obtained cluster (segmented area) from the k-means approach is superimposed with the preprocessed image to get the original skin lesion, shown in Figure 3(d).

2.2 FE module

Three different types of features, such as CM of up to 4th order, GARCH, and texture features by LBP, are discussed in the following sub-sections.

2.2.1 Colour moments

As the colour of the skin lesion plays an essential role in dermoscopic image classification, CMs are an excellent feature to use. They are computed from the RGB colour model, which characterizes the colour distribution in the skin lesion regions. Table 1 shows the CMs used in DICS. The total number of colour features extracted in this phase is 12.

Table 1. Colour moments

2.2.2 GARCH features

The GARCH model was developed by Bollerslev [13]. As the name implies, GARCH includes a feedback mechanism to predict future variances using past variances. The robustness of a GARCH model depends on several factors, including the quality and quantity of data used to estimate the model parameters and the choice of model specification. It requires sufficient high-quality data to accurately estimate the model parameters. The proposed system uses high-quality dermoscopic images from the PH² database. It is also essential to carefully consider the assumptions and limitations of model specification and the best fits the data are chosen. The GARCH model has been applied to speech processing [14, 15], image de-noising [16], and signal classification [17]. If the intensity information in each colour channel is modeled with GARCH, then the parameters of GARCH can be used as features. The GARCH (p,q) process is defined by:

$\varepsilon_t=\sigma_t z_t$          (8)

where, σ_t is the conditional standard deviation and z_t is a random variable drawn from a Gaussian distribution with zero mean and unit variance. Also, σ_t is a process such that:

$\begin{aligned} & \sigma_t^2=\alpha_0+\sum_{i=1}^q \alpha_i \varepsilon_{t-i}^2+\sum_{j=1}^p \beta_j \sigma_{t-j}^2 \\ & \text { where } q>0 ; p \geq 0 ; \alpha_0>0 ; \alpha_i \geq 0 ; \beta_j \geq 0 ; 1 \leq i \leq q ; \\ & 1 \leq j \leq p ; \text { and } \sum_{i=1}^q \alpha_i+\sum_{j=1}^p \beta_j<1\end{aligned}$          (9)

In this study, the parameters of the GARCH model with units p(1) and q(1) are used. That is, GARCH(1,1):

$\sigma_t^2=\alpha_0+\alpha_1 \varepsilon_{t-1}^2+\beta_1 \sigma_{t-1}^2$          (10)

The constant parameter in Eq. (10) such as α₀, α₁and β₁ can be obtained using maximum likelihood estimation [17]. These parameters, along with the conditional mean constant, are considered features. To fit the GARCH model, the condition means are assumed to be constant for simplicity as GARCH processes are independent of conditional mean on the past, and it depends only on the conditional variance on the past. The total number of GARCH features extracted in this phase is 12.

2.2.3 LBP features

LBP is an efficient local descriptor in many pattern recognition approaches [18, 19]. It is invariant against any intensity changes due to illumination variances. It uses simple thresholding of neighboring pixels inside a 3x3 window with a threshold value equal to the intensity of the centre pixel. LBP features are computed by using the following Eq. (11):

$L B P\left(X_c, Y_c\right)=\sum_{n=0}^7 2^n S\left(i_n-i_c\right)$          (11)

where, (X_c, Y_c) is the centre pixel of a 3x3 window, S is the threshold function i_c and i_n represents the value of the central pixel and the value of n^th neighborhood pixels, respectively. The threshold function S is defined by:

$S= \begin{cases}0 & \text { if } i_n-i_c<0 \\ 1 & \text { if }  i_n-i_c \geq 0\end{cases}$          (12)

4.png

Figure 4. LBP Process

The application of Eq. (11) the whole image produces LBP with 256 different patterns. All these patterns are extracted for each colour channel without considering the border. The total number of LBP features extracted in this phase is 768. Figure 4 shows the working process of LBP. The value inside the parenthesis is obtained by Eq. (12).

2.3 FS module

The DICS extracts the features such as CM of up to 4^th order, parameters of the GARCH model, and texture features by LBP to classify skin lesions. The feature dimension is 792 per dermoscopic image. The extracted features by the above steps from 2.2.1 to 2.2.3 may contain redundant information, which may degrade the system's performance. To overcome this, the FS module is introduced in the DICS. A commonly used hypothesis test, the t-test, is used to decide whether the extracted features are significantly different from each other or not. The features are ranked by class separability criteria to identify the dominant feature set.

Let us consider, N={N_k1, N_k2, ..., N_kf 1≤k≤n₁} be the extracted features from n₁ several normal training samples and A={A_k1, A_k2, ..., A_kf 1≤k≤n₂} the extracted features from n₂ several abnormal training samples. Each normal and abnormal case consists of f number of features. Let the means of each feature in groups N and A are defined by:

$m_{N_i}=\frac{1}{n_1} \sum_{k=1}^{n_1} N_{k i} \quad 1 \leq i \leq f$          (13)

$m_{A_i}=\frac{1}{n_2} \sum_{k=1}^{n_2} A_{k i} 1 \leq i \leq f$          (14)

The t-test statistic value is given by:

$t_i=\frac{m_{N_i}-m_{A_i}}{\sqrt{\frac{V^2}{n_1}+\frac{V^2}{n_2}}} \quad 1 \leq i \leq f$          (15)

where, $V^2=\frac{\sum_{x=1}^{n_1}\left(N_{x i}-m_{N_i}\right)^2+\sum_{y=1}^{n_2}\left(A_{y i}-m_{A i}\right)^2}{n_1+n_2-2}$.

It is observed from Eq. (15), that the t-test provides f number of t values. The feature with a high t-value indicates that the two classes (normal and abnormal) are significantly different. In this stage, only the dominant elements are selected by t-test, and the classification is done in the next step by the selected features as Dominant Features (DF) in a predefined percentage (1%, 2%, and 3%) of total feature dimension.

2.4 IC module

This module employs a deep learning approach for dermoscopic image classification. Deep learning adds more hidden layers between input and output layers to nonlinear model relationships, whereas traditional neural networks normally have one hidden layer [20]. A simple deep-learning architecture with two hidden layers is shown in Figure 5.

Let us consider X=[x₁, x₂, ..., x_m] be the input feature vector with m features that forms the input layer of the neural network. The i^th neuron in the j^th hidden layer is denoted by, $h_i^j$ and the corresponding weight is denoted by $w_k^{i, j} \quad 1 \leq k \leq m$. The output of the first hidden neuron in the first hidden layer, i.e., $h_1^1$ is obtained by the dot product of features in the input feature vector X with their corresponding weight. It is defined by:

$z=\sum_{i=1}^m w_i^{1,1} x_i+$ bias          (16)

The obtained dot product of $h_1^1$ is fed into an activation function $h_1^1=f(z)$ to get the neuron's output. This procedure is repeated for all neurons in each hidden layer. As the step function has no useful derivatives, the tansig function is used in this study as an activation function and works better for backpropagation. The backpropagation is a descent algorithm that minimizes the error at each iteration when the error signal is backpropagated to the lower layer. The weights in this network are adjusted, so that error rate reduction usually occurs in a decent direction.

5.png

Figure 5. Neural network architecture with two hidden layers

6.png

Figure 6. Tan-Sigmoid function in the classifier input layer

The training is usually done based on the error signal with an iterative updating of the weights, which uses the mean-squared error function of the negative gradients. The error signal is the difference between the actual and the desired output values multiplied in the input layer by the sigmoid activation function $\left(a=\operatorname{tansig}(n)=\frac{2}{1+e^{-2 a}}-1\right)$.

The tansig function is a popular activation function mainly used in the hidden layers. It has several desirable properties. First, it is a smooth, continuous function that can be easily differentiated and used in gradient-based optimization algorithms. Second, a nonlinear function allows neural networks to model complex relationships between inputs and outputs. The choice of activation function, including the tansig function, depends on the nature of the data and the task at hand. The tansig function is a hyperbolic tangent function that maps its input to a range between -1 and 1. It accelerates the convergence of the backpropagation method and is well-suited for modeling nonlinear relationships between inputs of normal and abnormal dermoscopic images. Figure 6 shows the tansig function in the classifier input layer.

The classifier output is obtained from the results of j^th the hidden layer. Since the classifier is trained to produce only one output, either normal or abnormal, the weight w⁰ consists of n weights with $w_1^o, w_2^o, \ldots w_n^o$. It is defined by:

$z=\sum_{i=1}^n w_i^o h_i^j+$ bias          (17)

The obtained dot product is fed into an activation function $\hat{y}=f(z)$ to get the final classifier output. The function used in the output layer is the linear function (a=purelin(n)) which is shown in Figure 7. While training the classifier, the abnormal inputs are given with class labels as '1' and normal with class labels '0'. Thus, for a testing sample, if the classifier output is greater than 0, the testing sample is classified as abnormal else it is classified as normal.

7.png

Figure 7. Linear function in the classifier output layer

After the FS module, the deep learning network is trained using the selected features with ten hidden layer structures and validated using a ten-fold cross-validation technique.