Modern Techniques for Computer-Aided Melanoma Diagnosis

1. Introduction

During the last two decades we have observed an amazing development of imaging techniques targeted for bio-medical applications. New technologies, devices and equipment combined with sophisticated analysis, visualization and storage algorithms as well as programs designed for acquisition of data, signal processing, database building or data manipulation are now available commercially.

Tremendous advances were also observed in the domain of dermatology. Two types of specialized equipment have been developed for this purpose, namely the dermatoscope and the SIAscope. The dermatoscope or epiluminescence microscope has become a standard tool (Argenziano et al. 1998). Connected with a digital camera and lighting equipment providing light with specific wavelengths, it becomes a powerful diagnostic device. The SIAscope can be used to acquire images under different lighting conditions while the standard dermatoscope takes the images under single light wave-length.

Images acquired using dermatoscopes can be stored and analyzed further by digital computers. Some of such applications for clinical use are based on personal computers or laptop platforms (Burroni et al. 2004,). Newer versions can be connected even to an iPhone or hand-held device such as a palmtop computer. Several companies produce hardware solutions for dermatologic diagnosis – one of such popular solutions is proposed by DermLite.

In both techniques the digital image serves as a basis for medical analysis and diagnosis of lesions under consideration. As there is a general lack of precision in human interpretation of image content, advanced computerized techniques can assist doctors in the diagnostic process (European Consensus 2009). Several companies offer complex software solutions – let us mention here only a few, such as: Mole Expert, MoleMAX, and DDAX3 (see the websites listed under Dermatoscopic systems in the list of references).

Analysis these solutions, targeted as an aid for the practicing medical doctor, reveals they all contain data storage and manipulation software (multi-media database). Some provide also some basic tools for skin lesion images analysis. Great success in application of the dermoscopic equipment and software could be noted in Australia – one of the courtiers with highest incidence of malignant melanoma registered in the last few years. The Solarscan system has been in use in almost all dermatology clinics across Australia (compare Talbot and Bischof, Menzies et al. 2005).

In this chapter we review some of the now standard methods of image processing for melanoma diagnosis and also introduce new identification methods based on color decompositions used in video and TV image processing as well as new spatial and frequency information found in the skin texture. These new methods assume that some neighborhood properties of pixels in dermatoscopic images can be a sensitive probe of different pigmented skin atypia and the melanoma progression. The most promising methods are those decomposing different frequency scales of the texture. Such multi-resolution analysis is well suited to determine distinctive signals characterizing the class of the image. We present also classification analysis of the pigmented skin lesion images taken in white light based on the inductive learning methods proposed by Michalski (AQ). These methods are developed for a computer system supporting the decision making process for early diagnosis of melanoma. Symbolic (machine) learning methods used in our study are tested on two types of features extracted from pigmented lesion images: color/geometric features, and wavelet-based features. Classification performance with the wavelet features, although achieved with simple rules, is very high. Symbolic learning applied to our skin lesion data outperforms other classical machine learning methods.

2. Clinical diagnostic approaches for diagnosis of melanoma

Analyzing the medical literature and textbooks one can find out that there are several standard approaches for analysis and diagnosis of cutaneous lesions (Johr 2002, European Consensus 2009). Let us shortly describe the characteristic features used in each method for classification of lesions.

2.1. Menzies scale

This scale is used to differentiate melanoma lesion from a non-melanoma one. Lack of negative features and the presence of at least one positive attribute suggest that the lesion should be diagnosed as melanoma. The features used for differentiation are presented in table 1.

Negative features	Positive features
lesion axial symmetry; lesion color symmetry; the presence of one color	blue-white veil; numerous brown dots; pseudopodia; radial streaks; discoloration of a scar; black dots-globules at the periphery; many colors (5 or 6); numerous blue / gray dots; extended network

Table 1.

Menzies scale feature table.

2.2 Seven-point scale

This scale uses major and minor criteria to grade lesions in scale from 1 to 7. The presence of any major criteria adds two points, one point is for minor criteria. To diagnose a melanoma using this scale the lesion must score at least 3 points. Table 2 contains the scoring criteria.

Criteria	Points
Major; atypical net pigmentation; atypical vascular pattern; blue-white veil	2; 2; 2
Minor; Irregular streaks(radial streaks); Irregular pigmentation; irregular spots / globules; areas of regression	1; 1; 1; 1
Score	< 3 = non melanoma; ≥ 3 = suspected melanoma

Table 2.

The 7-point scale criteria and scoring table.

2.3. TDS score based on ABCD rule

TDS (Total Dermoscopy Score) – is a uniform system used for dermatoscopy assessment. Using a linear equation the ABCD rule (Nachbar et al. 1994) introduced the option to grade skin lesion depending on the degree of their malignancy. This degree is determined by the TDS value calculated from the following equation:

$$TDS=A\ast 1,3+B\ast 0,1+C\ast 0,5+D\ast 0,5$$E1

The ABCD (Asymmetry, Border, Color, Dermoscopic structures) rule is used for diagnosis of skin changes of melanocytic origin. It is used to assess lesion and brings an answer to the question whether it is a mild change, suspicious or malicious. The variables for the equation (1) are determined by visual assessment of: A - lesion shape asymmetry and color asymmetry, B – border shape and sharpness, C – presence of various colors (red, blue-gray, brown, black, white), D – presence of dermoscopic structures such as pigmentation net, regression regions, dots, globules and so on.

2.4. ABCDE rule

The ABCDE (Asymmetry, Border, Color, Diameter, Evolution) (Thomas et al. 1998) (Rigel et al. 2005) is a rule that evolved on ABCD scale used to calculate TDS value. According to this rule, the lesion is suspicious if visual assessment of the lesion is positive on any of the following features: A - lesion shape asymmetry and color asymmetry, B – border shape and sharpness, C – presence of various colors (red, blue-gray, brown, black, white), D – diameter of the lesion is greater than 6mm, E – elevation of lesion has grown over short period of time, or lesion has evolved rapidly.

3. Definitions and discussion of measures of performance in the diagnosis process

Skin cancer is a fast developing disease of modern society, reaching 20% increase of diagnosed cases every year. In Central Europe the incidence rate is 10-14 per 100,000 population and in Southern Europe the incidence is 6-10 per 100,000; in the USA 10-25 per 100,000; and in Australia, where the highest incidence is observed, 50-60 per 100,000.

Dermatoscopy is primary and commonly used method of diagnostics for nearly thirty years. This method is non-invasive and requires great deal of experience to make correct diagnosis. As described in (Menzies et al. 2005) only experts have 90% sensitivity and 59% specificity in skin lesion diagnosis, as shown in table 3. The result of strict compliance with the instructions should be similar to this document. Specificity and sensitivity are indexes calculated from the following equations:

$$Sensitivity=\frac{TruePositive}{TruePositive+FalseNegative}$$E2

$$Specificity=\frac{TrueNegative}{TrueNegative+FalsePositive}$$E3

	Sensitivity	Specificity
Experts; Dermatologists; Trainees; General practitioners	90%; 81%; 85%; 62%	59%; 60%; 36%; 63%

Table 3.

Sensitivity and specificity.

Variables like TP, TF, FN and FP are described in table 4. That is why computer aided diagnostic gains more popularity every year and more complex algorithms are being developed. Some of this algorithms use: colorimetric and geometric analysis, statistical differentiation based on knowledge databases of diagnosed cases (supplied with most video and photo dermatoscopes).

	Actually Abnormal	Actually Normal
Diagnosed as Abnormal; Diagnosed as Normal	True Positive (TP); False Negative (FN)	False Positive (FP); True Negative (TN)

Table 4.

Diagnosing possibilities.

The diagnostician has to perform difficult classification task and take the decision of telling the patient that he is in good health or not, basing on available information obtained in the process of clinical examination.

It is possible to apply computerized image processing, feature extraction and classification methods to improve diagnosis process and assist the physician in making the decision.

4. Construction of image classifiers

In process of images classification a number of classification methods (Chan et al. 1999) can be used. For the dermoscopic image differentiation, classification algorithms such as Artificial Neural Networks (Park et al. 2004) and Support Vector Machines (Wu et al. 2008) are used. Support Vector Machines (in short: SVM) are commonly used with four types of kernels: linear, polynomial, radial and sigmoid (Chang et al. 1998), each for different classification effectiveness. As for classification using Artificial Neural Networks (in short NN), the radial basis function is preferred. NN are composed of simple elements (neurons) that are inspired by biological nervous systems. These neurons are working in parallel creating a network that is trained to perform connections between elements. Before NN can be used for classification, the network need to be trained, so that when it is given certain input it will respond with desired output.

Radial basis function, known also as Gaussian function is used often due to its strong classification power. RBF consists of two layers as shown in figure 1. First layer is a hidden radial basis layer of S1 neurons.

Figure 1.

Radial Basis Function Neural Network.

The second layer is a linear basis layer of S2 neurons. In figure 1 presents a NN working with Radial Basis Functions.

Conventional neural networks are performing empirical risk minimization where network weights are determined by minimizing the error between actual and desired outputs. SVM are based on the structural risk minimization where the parameters are optimized by minimizing classification error (see http://www.kernel-machines.org/).

5. Image acquisition process and image pre-processing

While working with digital images of any kind, the image acquisition is always the first step of image analysis process. Image acquisition is a stage, at which images are collected in order to create certain data set, which is later analyzed to see if any of the gathered images share similar features or contain any of the predefined features, or meet previously defined assumption.

5.1. Image acquisition

Dermoscopic images are basically digital photographs/images of magnified skin lesion, taken with conventional camera equipped with special lens extension. The lens attached to the dermatoscope acts like a microscope magnifier with its own light source that illuminates the skin surface evenly. There are various types of dermoscopy equipment, but all of them use the same principle and allow registering skin images with x10 magnification and above. Due to light source integrated into dermatoscope lens, there happens to be problem with skin reflections. To counteract this problem, a liquid is used as a medium layer between the lens and the skin. In modern dermatoscope the liquid is not necessary, because of the polarized light source that removes the reflection problem.

Digital images acquired using photo dermatoscope are sufficiently high resolution to allow for precise analysis in terms of differential structures appearance. Dermatologist can create accurate documentation of gathered images, opening a path for computer analysis, where images are processed in order to extract information that can later be used to classify those images.

5.2. Image preprocessing

Before analysis of any image set can take place, pre-processing should be performed on all the images. This process is applied in order to make sure that all the images are consistent in desired characteristic. When working with dermatoscopic images, pre-processing can cover number of features like: image illumination equalization, color range normalization, image scale fitting, or image resolution normalization. This can be dependent on defined prerequisites and methods applied in post processing.

An example of elementary operation such as image normalization is the resolution matching. Assuming that the image size in pixels is given, and all images are in the same proportion (e.g. aspect ratio of 4:3), it is easy to find the images of smallest resolution and then scale the larger images to match the size of the smallest one. This operation allows calculating the features like lesion dimensions, lesion border length and lesions area coverage. It is possible to normalize the other parameters like color palette normalization, color saturation normalization, normalization of color components, and so on. Very common operation in preprocessing is color components normalization, known as the histogram equalization.

Image histogram is the distribution of colors values in between extreme colors used in the palette. Assuming the situation where the brightest points of the grayscale image are not white and the darkest points are not black, performing histogram equalization will redistribute all the colors of the image in a way that brightest spot of the processed image will be color and the darkest regions of the image will become black. Figure 2 illustrates the results of histogram equalization performed on a dermoscopic image.

Figure 2.

a) original image and its corresponding histogram, b) image after histogram equalization and its corresponding histogram.

Normalized image is characterized by better brightness, sharpness and color depth, thereby allowing easier separation the lesion area from the background (skin color). Applying histogram equalization operation allows for better differentiation of image detail, and thus improves the efficiency of features extraction. Histogram equalization can be performed for each of the color components separately, or on all of the components at once.

5.3. Filters

Filtering an image consists of application of a certain transformation or an algorithm to image data in order to modify certain part of it. The most common task of a filter is to separate redundant information from the relevant data. By using simple filters it is possible to sharpen image, blur image, change color, etc. Applying more complex filters is possible to enhance more important sections of the image. An example for this can be strengthening particles or detection edges.

In image processing one certain transformation is most notable, namely the binarization. It is a point transformation that is one of the basic operations used when processing any image. This operation based on gray-scale image and given threshold value will outputs a binary image, which uses only two colors (black and white) to represent data. To achieve the binary image a threshold needs to be defined. This threshold is used to determine which points of the original image will be converted to black and which to white. The binarization process is illustrated by figure 3. Having a binary image allows performing various operations like measuring the length, performing segmentation, determining the number of elements, etc. Binary image is easy to work with, as most operations are fast to perform and do not require long computing time. It is quite easy to analyze the shape of objects, calculate symmetry, or find the center of weight of an object.

Figure 3.

a) original image, b) binarized image using threshold value of 60, c) binarized image using threshold value of 100, d) binarized image using threshold value of 200.

A simple example illustrating the idea of using filters in dermoscopic images can be the removal of artifacts such as hair (figure 4). To do this, a maximum filter can be applied on a binary image. Maximum filter process each point of an image and calculate its surrounding (neighborhood) in determined distance, then the maximum value is entered in the output image. As a result, the small details of the image are removed depending on the neighborhood size. The use of this filter on the single object will reduce it in size.

Figure 4.

Hair removal. a) original image, b) binarized image, c) image with removed artifacts.

Performing binarization on dermoscopic image can sometimes create errors in output, may it be due to poor illumination of the corners, or hairs on the image. Dealing with unevenly illuminated image is another practical example of filtering unwanted data from dermoscopic image. For this purpose a specific algorithm is used and results of its application can be viewed in figure 5. This algorithm’s task is to check whether there is an object in the middle of the image, and then whether each of the corners of the image contains different object in contact with the edges of the image. If it happens that all image corners contain objects then they are removed. Condition for the proper removal of redundant objects is that on the image, there must be at least two different objects, and each of the points in the corners must be a part of one of the objects. If all conditions are met, the filter removes all points belonging to object in contact with the corners of the image. This filter can only work with images that cover entire lesion. If image contains only a fragment of the lesion, this filter will not modify the input image.

Figure 5.

Filtering of vigneting artefacts. a) original image, b) binarized image, c) image with filtered corners.

6. Numerical methods for image analysis

Several types of image analysis can be done in an automatic way. Specifically geometric features and color features/decompositions can be easily performed and are found in many commercial image analysis and dermoscope accompanying programs.

6.1. Geometric features

Many attempts have been made to automate the process of dermoscopic image feature extraction, and a number of algorithms have been developed to do this. Most common geometric features that are calculated are lesion dimensions, borders shape, lesion symmetry and color symmetry lesion area. More complex problem like various differential structures recognition still exist and is a subject of current image processing research.

Calculation the dimensions of lesion on the dermoscopic image is an easy task. The easiest way to determine the surface size is to count all points that are marked as the lesion on the binary image achieved in process of binarization. The sum of these points and its percentage cover of an image lesion image can be used as variables in classification process. These two factors can be used to calculate other factors like various color coverage of the lesion.

Border length of the lesion is another feature that is easily calculated from a binary mask of the dermoscopic lesion image. Border length can be calculates simply by counting all boundary points on binary lesion mask. It is possible to estimate how irregular the border is by comparing it to the lesion size, e.g. circular lesion would have shorter border length comparing to irregular lesion of the same size.

The symmetry of the lesion is an important factor in diagnostic process. Calculating symmetry value as a geometrical feature requires more complex operations, as it is necessary to estimate proper axis of symmetry. This task can proved to be complicated, even with the manual axis alignment. The examples images for determining the symmetry axis are illustrated in figure 6.

Figure 6.

Examples of lesion shapes.

Using the algorithm of smallest bounding box estimation prove to be effective way to determine the symmetry axis for various lesion shapes. Symmetry axis for calculated bounding rectangle is used as the symmetry axis for the lesion. As shown on the figure 7.

Calculating the value of the symmetry is very simple and is done by checking if every point of the lesion on one side of the symmetry axis has its corresponding point on the other side of this axis. Symmetry is calculated separately for vertical and horizontal axis, and then the average value is subtracted from the lesion size value. The percentage of asymmetry can be calculated from this data and later used in the classification process.

Figure 7.

Different lesion shapes and its corresponding bounding boxes and estimated symmetry axis.

6.2. Color decomposition

Skin lesions are difficult to classify because of their short color ranges, instead of real-world images (Deng et al., 2001). Malignant melanoma is a kind of skin cancer that has some characteristic color groups like: black, blue-grey, red, light brown, dark brown and skin color. These colors appear on images and can depend on cancer progress stage, lesion depth and blood vessels. These colors are also used in melanoma diagnosis scales like ABCD.

Typically in digital images we use RGB (additive) or CMYK (subtractive) color representations.

There exists however several other decompositions used in video or TV coding. These color decompositions have not been used so far for melanoma diagnosis – they were first proposed as useful features in our works (Ogorzałek et al. 2009, 2010). Typical melanoma image color decompositions are shown in two color spaces in figure 8. It can be easily seen that in the second picture the points representing colors are grouped in two areas and thus could be distinguished/classified very easily.

Figure 8.

Three different color decompositions of dermoscopic image. Every point corresponds to a region of the images. HSV decomposition (middle picture) gives the best classification possibilities.

6.3. Image segmentation

Image segmentation is performed in two steps: color quantization and spatial segmentation. In the first step image is quantized to several region classes. Quantization is done in color space only. After that a class-map is created by their belonging to the corresponding color class. Second step is spatial segmentation is performed on created class-map. Color similarity is in second step not taken into consideration. The next step is objects extraction (Liu et al. 2005). Described method extracts borders from a segmented image. By using these borders objects in the image can be extracted. Object extraction is done by a region growing algorithm. This is the most time consuming part of presented method. The result of the segmentation process is illustrated in figure 9.

Each object is represented in four color spaces: RGB, HSV, NTSC and YCbCr. Every color space is represented by few variables like: red, green and blue in RGB. There are 12 color spaces features.

7. Analysis of images based on Wavelet transformations

Our database of dermoscopic images (2272x1704, 24bit) collected with Minolta Dimage Z5 camera with epiluminescence lens consisted of 78 cases of melanoma and 80 cases of dysplastic nevi (taken at the same conditions). The selection to the categories: melanoma (label=1) and dysplastic nevus (label=0) was done on the grounds of biopsy and histopathologic examination (Żabińska-Płazak et al. 2005).

Figure 9.

a) original image, b) initial segmentation of the image, c) binary segmentation mask, d) final image segmentation.

The preparation stage for geometric/coloristic feature set (Fikrle et al. 2006) (Ogorzałek et al. 2009) was done in the following way:

All the images were scaled to 800x600 pixels, normalized in the quantization depth by a histogram correction, and removed the background (in dermoscopic images ‘background’ means skin of normal complexion). Identification of the mole spot was done by color indexing i.e. reduction of the quantization depth. Then artifacts like obscuration of the image corners or a pitch on the lesion were extracted. In the end with opening/closing operations (and binarization when required) we got an output form of the lesion image.

The preparation stage for the wavelet-based features (Patwardhan et al. 2003, 2005) (Surówka et al. 2006) was done in the following way:

Since the wavelet decomposition downscales the input image by a factor of 2 in rows or columns every iteration, the width and length had to be powers of 2, so the source images were filled with black pixels in due number of rows and columns. In order to apply any transform on the picture, the compressed image (JPG) must be changed into a 2D matrix of numbers. The usual way of extracting the numerical values of pixels consists in changing the RGB representation into a one-dimension color index.

7.1. Wavelet-based features

Skin texture analysis based on the multiresolution wavelet-based decomposition of the dermoscopic images is a new method to obtain valuable features for melanoma and dysplastic nevus classification. Wavelet transforms have been widely studied as tools for a multi-scale pattern recognition analysis. Wavelets are functions well localized both in space and frequency. The wavelet theory is closely related to the theory of digital filtering (Kadiyala et al. 1999, Mallat 1989) so the properties of the decomposition filters (the choice of the basis, degree of regularity, and the sub-bands of interest) play an important role in texture characterization.

Performance of the wavelet analysis of images depends on various factors:

• For feature extraction it is important to choose between recursive decomposition of the low-frequency branch (Mallat 1989) (the pyramidal algorithm) and, on the other hand, a more selective tree-structured analysis where the further decomposition is applied to the output of any channel (Patwardhan et al. 2003, Chang et al. 1993). The latter is advocated by observations that the significant sub-bands of the pigmented skin texture are those in the middle frequency range (Patwardhan et al. 2003, 2005).

• The decomposition should be performed over an optimal finite range of resolutions (Chang et al. 1993, Mallat 1989). This issue addresses the problem of the wavelet order.

• Different wavelet bases have diverse impact on the texture classification (Kadiyala et al. 1999). They have certain constraints and should, in principle, be orthonormal (energy preserving, non-redundant). Biorthogonal (near orthogonal) bases can produce, on the other hand, more compact and symmetric wavelets and the requirement for orthogonality is not superior (Mallat 1989).

• Images are two-dimensional signals so separable or non-separable sampling is possible. In the former case a one-dimensional decomposition algorithm is applied to the rows and to the columns of the image. The resulting signal is a tensor product of the separate 1D filter. The non-separable case, on the contrary, is based on lattice sampling. Both representations have advantages and disadvantages.

Unfortunately, up to now only a limited number of the aforementioned possibilities have been tested on skin dermoscopic textures (Patwardhan et al. 2005, Patwardhan et al. 2003, Surówka et al. 2006). Since skin textures are signals consisting mainly of middle frequency bands it has been shown that the most adequate approach to the wavelet-based multiresolution analysis is the concept of wavelet packets (Chang et al. 1993). In this approach an image is decomposed along dominant frequency channels forming a parent-child structure of a tree. A certain path of the tree is chosen according to the average energy content of the channel referenced to the highest energy of a channel on the same decomposition level. The determined subbands form a feature set for classification. Such an approach is shown in (Chang et al. 1993). The main bias of the method in (Chang et al. 1993) is the threshold value of the channel energy that triggers subsequent decompositions. Since features used for selection of melanoma and dysplastic lesions are clustered into rough sets, the structure of the tree may be affected by an arbitrary choice of the threshold. In the adaptive wavelet-based tree structure analysis (ADWAT) (Patwardhan et al. 2003,2005) decision about the further decomposition of the channel is taken by statistical analysis of average energy, maximum energy ratio and fractional energy ratio among all the features at each level of decomposition. The resulting tree structured models of melanoma and dysplastic nevus are patterns for semantic comparison with unknown skin lesion images. The sensitivity and specificity of this method is further improved by introducing Gaussian and Bell shaped fuzzy clustering of the features (Patwardhan et al. 2005). Work (Surówka et al. 2006) is based on the idea presented in (Patwardhan et al. 2003) to perform a tree-like selective wavelet analysis of the skin lesion images. Unlike in (Patwardhan et al. 2003) the feature set is built from all of the analyzed channels and the decision to include a given branch of the tree does not come from ADWAT. Instead, all the wavelet based features are input to the Ridge classification model (Merkwirth), which determines a minimum–bias optimal vector of features. Those selected features can be used in various machine-learning methods.

After the preparation stage described earlier, the 2D=1Dx1D wavelet transform was applied to the values of pixels extracted from the discussed dermatoscopic images. The class of the filter was Daubechies 3 (Daubechies 1992) and its efficient algorithm was taken from (Numerical Recipes in C 1999). The choice for the filter was made concerning simplicity, performance and possible comparisons with (Patwardhan et al. 2003).

One iteration of the wavelet algorithm produces 4 subimages which can be considered as LL, LH, HL and HH filters, where L and H denote the respective low-pass and high-pass filters. One subimage is a product of the wavelet transform acting on each raw and then on each column of the parent image. Each iteration reduces half of the rows and half of the columns from the parent image. The algorithm presented in (Numerical Recipes in C 1999) is a pyramidal one i.e. it decomposes the image along the LL band in each iteration. An example of the decomposition process with 3 iterations is shown in figure 10. Here the smooth and the detail part resulting from any iteration can be easily recognized.

Figure 10.

Multi-band images i.e. a set of images of the same mole from its different frequency segments. Here recursive decomposition of the low-pass band takes place. For the sequence of the pyramid see next figure. Note that the luminosity of this summary image is offset otherwise it would be too dark for a presentation

According to the idea of the wavelet analysis, we adapted the pyramidal algorithm to a selective analysis of any sub-band of the parent image. Altogether 3 iterations were calculated. The proper sequence of filters in that decomposition is shown below. Every iteration has the same arrangement of filters with (sub) labels denoting 1=low-low, 2=low-high, 3=high-low and 4=high-high sub-band.

Figure 11.

Sequence of the pyramid construction.

Following the original idea (Chang et al. 1993) and its extended applications (Patwardhan et al. 2003, 2005) a feature set at a given level of decomposition was built. One iteration (resulting in 4 subimages) produced 11 coefficients calculated from the energy of the pixels (e₁, e₂, e₃, e₄), the maximum energy ratio (e_i/e_max, i=any three out of four except the subimage with the maximum energy) and the fractional energy ratio (e₁/(e₂+e₃+e₄) + its two permutations). Term ‘energy’ means here the sum of the absolute values of the pixels, where M and N denote the actual dimensions of the sub-image.

Energy is calculated according to the following equation:

$$e=\frac{1}{MN}{\displaystyle \sum }_{m=1}^M{\displaystyle \sum }_{n=1}^N\left|f(m,n)\right|$$E4

Since the algorithm is applied to the image recursively, the number of coefficients amounts to (1+4+16)x11=231. Numbers 1, 4 and 16 in this calculation mean here the number of sub-images to filter in each iteration.

In the approaches presented in (Patwardhan et al. 2003, 2005) the calculated features were discriminated according to their distributions within a particular channel. Only features that generated bimodal distributions separating melanoma and dysplastic images were kept in the feature set. All other were rejected. In (Surówka et al. 2006) the features were not used to control the development of the tree structure signatures. Instead, all 231 coefficients were accepted as potentially significant discriminating signals.

Since the number of dermatoscopy images taken to the analysis was limited and the 231 coefficients affect the classification in a different way, one can additionally pre-select the coefficients by their importance. The criterion was to maximize the classification performance of 20 and then10 most significant features. In (Surówka et al. 2006) the selection was twofold and was done in a Matlab toolbox ENTOOL (Merkwirth). First one Ridge regression model with 40 penalty vectors was used to determine an optimal representation of the images in a numeric form. Ridge regression constructs a linear model$\widehat{y}=X\beta +{\beta }_0$ (X-data matrix, y-vector of categories), but instead of minimizing the sum of squared residuals${\left(y-X\beta +{\beta }_0\right)}^T\left(y-X\beta +{\beta }_0\right)$, it minimizes the regularized loss function (Tikhonov regularization):$RS{S}_{pen}={\left(y-X\beta +{\beta }_0\right)}^T\left(y-X\beta +{\beta }_0\right)+\lambda {\beta }^T\beta$.

The additional penalty $\lambda {\beta }^T\beta$shrinks the regression coefficients $\widehat{\beta }$towards zero, thereby moderately increasing bias while considerably decreasing variance of the constructed models. The penalty parameter$\lambda \ge 0$controls the amount of shrinkage and can be used to fine tune the bias-variance tradeoff. For this study, the optimal ridge penalty λ is automatically determined by Leave-One-Out Cross Validation on each training fold individually. To apply ridge regression to a binary classification problem, training outputs were coded as $y=1$ (melanoma),$y=0$ (dysplastic) and a threshold of 0.5 was applied to discriminate between both classes when doing predictions. Prior to the model construction, input variables were normalized by removing the mean and dividing by the standard deviation for each variable separately. Ten most significant features are presented below. They were obtained from a statistical analysis in ENTOOL (see Merkwirth). Indices denote the subchannels of the three different decomposition levels.

The ROC curve below presents the classification performance of 10 most significant features derived with help of 100 Ridge regression models. AUC is 95%.

average energy	$e_{\mathrm{4.1.2}},e_{\mathrm{2.1.1}},e_{\mathrm{3.3.2}},e_{\mathrm{3.3.1}}$
maximum energy ratio	$e_{\mathrm{1.2.2}}/e_{\mathrm{1.2.1}}, e_{4.2}/e_{4.1}, e_{4.3}/e_{4.1}$
fractional energy ratio	$e_{\mathrm{1.2.2}}/(e_{\mathrm{2.3.1}}+ e_{\mathrm{2.3.2}} +e_{\mathrm{2.3.4} }), e_{\mathrm{3.1.1}}/(e_{\mathrm{3.1.2}}+ e_{\mathrm{3.1.3}} +e_{\mathrm{3.1.4} })$

Table 5.

Coefficients calculated from the energy (the sum of the absolute values of the pixels): e₁, e₂, e₃, e₄, the maximum energy ratios: (e_i/e_max, i=any three out of four except the subimage with the maximum energy), and the fractional energy ratios (e₁/(e₂+e₃+e₄) + its two permutations).

Figure 12.

ROC curve presenting the classification performance of 10 most significant features derived with help of 100 Ridge regression models. AUC is 95%.

Following are the sets of images no which the system was tested:

Figure 13.

Selected images from two data sets used for testing. On the left melanoma images on the right dysplasty images.

7.2. Discussion

Classification models constructed using wavelet-based tree structure decomposition of the dermoscopy images are reported in the literature (Chang et al. 1993), (Patwardhan et al. 2003, 2005). In (Chang et al. 1993) channel decomposition is determined based on the ratio of the average energy of that channel to the highest average energy of a channel at the same level of decomposition. The decomposition is performed if this ratio is above a predefined threshold value. In the learning phase tree-like topologies of both classes are produced. In the testing phase unknown images are decomposed with the thresholds set in the learning phase to build tree structures that are classified due to the Mahanalobis distance from the known patterns. The maximum performance obtained is TPF=70% and FPF=20% for the learning phase of 15M+15D and the testing phase 10M+15D (M=melanoma, D=dysplastic nevus). Arbitrary selection of the maximum energy thresholds is the main drawback of this method.

In the ADWAT method (Patwardhan et al. 2003), like in this work, different combinations of channel energy ratios over a given decomposition level are used as additional features and a statistical analysis of the feature data (not solely the maximum energy) is used to find the threshold values. Based on histograms of the features only linearly separable and bimodally distributed features between image classes are further decomposed. The thresholds are selected to optimally partition the two image classes. In the testing phase unknown skin lesion images are semantically compared with the tree structure of the melanoma and the dysplastic nevus class.

Results from (Patwardhan et al. 2005) for the ADWAT method read: TPF=86.66%, FPF=11.11%, and the method itself is more robust than (Chang et al. 1993). Such results are cited for 15M+15D in the learning phase and 15M+45D in the testing phase.

In (Patwardhan et al. 2005) an extended ADWAT method is reported. In this approach three classes are defined: melanoma, dysplastic nevus and suspicion lesion. Unknown images are classified as melanoma only if the tree structure completely matches the pattern assigned to the melanoma image class. Incomplete tree structures are assigned to the suspicious lesion class and their trans-illumination images are analyzed in the second step. The trans-illumination images are taken in 580nm and 680nm light by the Nevoscope tool in an illumination mode that provides lesion depth and structural information. For the learning/testing scheme of 15M+15D and 15M+45D this method yields TPF=93.33% and FPF=8.88%.

When additionally fuzzy membership partitions of the melanoma and dysplastic lesion images by the Bell or Gaussian membership function are applied, better results are obtained: TPF=100% and FPF=4.44% (FPF=17.77%) for the Gauss (Bell) partition functions.

8. Machine learning of wavelet-based features of melanoma

We experimentally apply to the dermoscopic images the main global machine learning methods, namely neural networks, and support vector machines. They serve as a generalization engine in the classification task.

The selected machine learning methods were tested over the same set of samples (78 cases of melanoma + 80 cases of dysplastic nevus). Due to small statistics of individuals, for each classifier we applied the n-fold cross validation test to gain a reasonable estimation of its accuracy. The n-fold cross validation method randomly divides the set of all available patterns into n subsets. N different models are then trained using data from n−1 sets and validated using the remaining one (the holdout fold). The procedure is repeated for each of the n folds. Accuracy of both classifiers was measured by means of the Receiver Operating Characteristics (ROC) (Receiver operating characteristic (ROC) analysis 1998). The ROC curve shows the sensitivity, i.e. ratio of correctly identified melanoma lesions, in relation to the specificity, i.e. ratio of correctly identified dysplastic nevi. For the aforementioned methods we cite the area under the ROC curve (AUC) as the key indicator of the performance.

8.1. Neural network binary classifier

The MLP neural network classifier is easy in implementation and provides stable solutions. The topology of the network was subject to tests to determine an optimal configuration for maximizing its performance (both quality of the classifier and minimal training time).

The examined three-layer feed-forward back-propagated neural network had the following architecture. The input layer was composed of ten linear-output neurons with constant unity weights. The hidden layer made of ten neurons and an output layer formed by one neuron had logistic-like activation functions. In one training cycle some 150 000 iterations were performed until the network could classify the input with a defined precision.

To evaluate features according to their rank, which was crucial in effective teaching of the neural network, prior to the training of the neural network, the number of features was reduced by the Ridge model to twenty and ten (out of 231) most discriminating ones, to avoid the network overtraining effect.

The classifier was trained on the ten selected features. Ten most significant features selected with help of the Ridge linear model yielded sensitivity (TPF) of 89.2% and specificity (1-FPF) of 90% (Surówka et al. 2006). The cross validation set was shifted between the training runs. Bimodal distribution of the MLP-classified lesions is displayed in figure 14. As we can see, correctly identified lesions are those distributed around 1 on the ordinate and around 0 for the other ones. On average the classification performance yields sensitivity (true positive fraction) of 89.2% and specificity (1-false positive fraction) of 90%.

Figure 14.

Distribution of the recognized image classes. Input tags on the abscissa are: 1=melanoma, 0=dysplastic.

Lack of full separability of the two categories exhibits, as the major factor, the generalization error. Yet the inherent fuzzy nature of both feature sets also plays a role. The magnitude of both effects can be estimated by comparing appropriate performance factors above.

Due to statistics of our database we decided to perform experimentation with the oversampling technique. The image data were limited to the pigmented spot of the mole (background-free) and divided into 64x64 blocks of pixels. Each block was decomposed with the Daubechies-3 wavelets for all possible subbands of frequency (3 iterations). When the blocks were treated as individual images the accuracy of the feature set was estimated by AUC=87.8%. Segments grouped to their parent images yielded AUC=91.4%. This method appeared not to significantly improve the accuracy of the MLP-classifier and was a burden for our computational resources.

8.2 Support Vector Machine - C-SVM algorithm

The Support Vector Machine (SVM) technique (Vapnik et al. 1997, 1998) is well suited to search for an optimal binary classifier. In the Support Vector Machine approach input vectors are mapped into higher dimensional space to identify a hyperplane that best separates the positive and negative samples i.e. maximizes the margin between two classes. This principle is more intuitive than the MLP ‘black box’ approach. Since we cannot expect perfect separation of melanoma and dysplastic lesion images, the C-SVM algorithm was employed with the RBF kernel implemented in the osu-svm-3.0 toolbox for Matlab (Chang). The input data set was processed by scaling the values of each feature to [−1, 1]. The C-SVM approach with the radial base kernel (RBF) appeared to yield very good performance. Taking into consideration the whole feature set (231) the result for the scaled data was TPF=94.7% and (1-FPF)=95%. The most efficient discrimination between class melanoma and dysplastic nevus was found for the parameters C = 512 and = 0.000244, which were fixed by a grid search over the full parameter space. The feature set limited by the same Ridge model (Surówka et al. 2006) to the 10 most discriminating features yields a TPF=84.2% and(1-FPF)=85%, which is a very good result taking into account the reduction in the feature space. Tests with the polynomial kernel were not stable and yielded worse results.

In table 6 we show results of the above mentioned methods tested on the same data set.

	MLP	SVM RBS
#attributes	10	231	10
TPF; 1-FPF	89.2%; 90%	94.7%; 95%	84.2%; 85%

Table 6.

Classification performance of the MLP and SVM RBS methods tested on our data set. The quality of the methods is quoted in terms of TPF (true positive fraction) and FPF (false positive fraction). Tests were conducted with all 231 features and with 10 most important selected by the Ridge linear model (Surówka et al. 2006).

It must be stressed that unlike in (Patwardhan et al. 2005) we do not take advantage of any specialized instruments (Nevoscope) and our lesion images are taken in white light only.

9. Comparison of standard methods and feature sets with the new proposed approaches

Classification of the pigmented skin lesion images can vary from linear models through MLP, SVM, k-NN, to ensembles of models. Applicability of symbolic methods of classification to the problem of early detection of melanoma has not been thoroughly tested yet. However, one can show that inductive classification based on symbolic rules and knowledge inference (Meyrowitz et al. 1993) can be both simple and comprehensive for human understanding, and yield high classification performance to compare with ‘black box’ machine learning methods.

Symbolic methods of classification generate symbolic/declarative information about the classes. It is done usually by classification rules or sets/chains of notions/ideas inferred from the learning cases. On the contrary, classical numeric approaches as neural networks, fuzzy sets or genetic algorithms produce complex, less comprehensive (interpretable) knowledge (Michalski et al. 2006).

In this approach we study the classification methods developed by Michalski and applied in his software tool AQ21 (Michalski et al. 2006, Wojtusiak 2005). The learning phase of this classification i.e. construction of descriptions of the classes is the most important for the knowledge induction (Michalski 2004). In this step all positive and no negative examples form a description of a class (in the simplest way as an alternative of positive values and negation of counter-examples). Then such a definition is subject to maximally specific generalization to be able to classify new unknown examples (by elimination/addition of selectors, value range extension, value range closure, resolution, or negation extension). In the theory of Michalski class definitions should obey two conditions: completeness (each member of the class should comply with the class description), and consistency (if an object fulfills its class description, it cannot belong to other classes). Assigning an object to its proper class according to the generated rules can be represented in the form$D:\gg Y$ (D-description of class Y), where the description is a disjunction (an alternative) of the complexes:$D=C_1\vee C_2\vee \mathrm{...}{C}_n$. Such a form eliminates inseparable or incomplete descriptions by precedence of the class whose object membership probability is the biggest.

Reduction of the description, which is the goal of this calculus, leads to the maximum detailed set of complexes (called a star). To focus on the most promising candidate descriptions a heuristic LEF function (Lexicographic Evaluation Functional) is used.

Out of the lesion images the following features were extracted (altogether 45 features):

Set A:Geometry: area, perimeter, roundness, aspect ratio, fullness ratio, symmetry, border length, bounding rectangle.

Colors: average values and standard deviations, min, max for components {R,G,B}, {H,S,V}, {Y,Cr,Cb}, weighted color transforms, gray-blue, white, black, number of colors, area percentage.

Distributions: distance of the centre of gravity of a color component from the lesion centre of gravity and the bounding rectangle centre of gravity, zone: whitish, reddish, light-brown, dark-brown.

Set B:For wavelet features the (1D)x(1D) Daubechies wavelet packets (Mojsilovic et al. 2000) were applied sto the indexed images. An efficient pyramidal algorithm (Mallat 1989) (Chang et al. 1989) was adapted into a tree structure multi-level filter. In each iteration of the wavelet algorithm four subimages were produced: LL (upper left), LH (upper right), HL (lower left) and HH (lower right), where L and H denote the respective low- and high-pass filters (Porter et al. 1996). One iteration produced 11 coefficients calculated from energies of the four subimages (e1, e2, e3, e4), the maximum energy ratio (ei/emax, i=any three out of four except the subimage with the maximum energy) and the fractional energy ratio (e1/(e2+e3+e4) + its three permutations), where the term energy means the sum of the absolute values of the pixels normalized by the total number of pixels in the (sub)image (Patwardhan et al. 2003, 2005). Altogether 3 iterations were done yielding (1+4+16)*11=231 features. Both feature sets have been used separately.

AQ21 implements methods of the Attributional Calculus (Michalski 2004), so prior to the learning stage all continuous domains of the features have to be digitized by invoking the built-in ChiMerge algorithm (ChiMerge 5) (Kerber 1992). Since we had 231 features in the wavelet domain we decided to pre-select the most relevant ones. For feature selection we used the simplified ‘Promise’ method (Baim 1982) available in AQ21.

Pseudo-code of the simplified Promise algorithm looks as follows:

P = 0

For each value of the attribute v_i

S = {examples : e = v_i}

Find class C that has the largest number of examples in S

P = P+ #C ∩ S /#S

Return P

The ‘Attribute_selection_threshold’ parameter defining the minimum attribute discriminatory power was experimentally set to 0.9. This value limited our feature set B to less than 10 attributes out of 231.

The AQ21 program was run in two different learning modes: Theory Formation (TF) and Approximate Theory Formation (ATF) (Wojtusiak 2005). In the TF mode, learned rules are complete and consistent, whereas the ATF consists of rules optimized according to the Q measure that reads:$Q\left(R,w\right)=compl{\left(R\right)}^w\cdot cons{\left(R\right)}^{\left(l-w\right)}$, where the optimization parameter w may cause a loss of completeness (compl) and/or consistency (cons) but may increase the predictive power of the learned rule (R).

Due to limited statistics of individuals (78+80) the classification accuracy was tested by the n-fold cross validation: (testing(1M+1D)/learning(the rest of M and D).

A summary of other AQ21 settings used in our experiments are presented below.

Learning phase:

Consequent=[class=*],

MaxRule=1,

Cross_validation=80.

Testing:

Method=atest,

Threshold=0.5,

Eval_of_conjunction=min, Eval_of_disjunction=max,

Eval_of_selector=strict.

Overall results obtained for feature set A and B are shown in table 7. In case of multiple instances due to cross validation we quote below the average predictive accuracy as 1- (#errors+#no_class + #both_classes) /(#iterations).

	TF	ATF
Feature set A	88%	93%
Feature set B	94%	100%

Table 7.

Results of the AQ21 classification in feature set A and B

From table 6 we can see the quality of both feature sets. Our working hypothesis was that the geometric/colorimetric features, as more prone to local distortions/artifacts and also more sensitive to the segmentation process (Ogorzałek et al. 2009), should yield worse results than the wavelet-based signals which reproduce frequency modes of the skin lesion of melanoma and dysplastic-like origin. On this stage of work we cannot answer the question if the wavelet features carry better signals of melanoma than other direct measures, or they are just less biased by the image analysis process. However in view of the fact, that early stages of melanoma do not reveal apparent structures of malignancy it seems, that wavelet bases may probe those inherent signals.

Another conclusion is that methodology of Approximate Theory Formation (ATF) seems to fit the best to the inherent fuzzy nature of benign and malignant skin lesion images.

However comparison of the feature sets was not the primary goal of our study. Output generated from the tests includes learned rules consisting of selectors with a feature name, relation and the value/range. This is the full knowledge we can get from the symbolic classification methodology. We can observe the learned rules, and see the meaning and precedence of features (Meyrowitz et al. 1993). For instance from the feature set B we get an example rule:

[class=1]

# Rule 1

<-- [co188_M=0.421220..0.446530,0.520865..0.525190,0.596670..0.640195,0.642360..0.793450 : 78,0,100%,78,0,100%]

This rule gives a range of values where the selector “col88” (which corresponds to a fractional energy ratio e_2,3,4/(e_2,3,1+e_2,3,2,+e_2,3,3) has the total predictive power for the class “melanoma”. Of course this feature can be interpreted not only in terms of energy, but also structural interpretation may be deduced. The same analysis for the feature set A gives less predictive (94%) and more complex rules with some primary features like: “average level of red“ and „symmetry“ (Ogorzałek et al. 2009).

Comparisons with other machine learning methods tested on the same data set and analysis of the published results show that performance of the presented approach may be better than other learning methods and is promising in view of further research with better statistics of the melanoma data (Sboner et al. 2003)(Burroni 2004)(Surówka et al. 2007).

10. Performance analysis

78 (80) dermoscopy images of the melanoma (dysplastic) lesions, all confirmed by histological examinations, have been classified using a wavelet-based set of features. Discriminant power of those features has been determined by either Ridge regression models, or the ‘promise’ algorithm, and generalized in a three-layer back-propagated neural network/support vector machine, and by the Attributional Calculus.

Our results, presented in table 6 and table 7, confirm that neighborhood properties of pixels in dermoscopy images record the melanoma progression and together with the selected machine learning methods may be important diagnostic aids. Especially the inductive learning approach (AQ21, table 7) shows great classification potential.

The present accuracy of the classifiers is a promising factor for the further research in the field, especially considering the fact that all dermoscopy images were taken in white light without any specialized instruments. The constructed model of features for melanoma and dysplastic lesions can evolve, after some fine-tuning and improvements, to a computer-aided diagnostic system.

In future work, the most important factor is increase in statistics of the dermoscopic images. This is the most critical factor due to low rate (1-2%) of melanoma patients in the population (~1500) examined for the sake of this study.

Diversity of images should also cast a light on the role of image acquisition and quality.

Since the features of melanoma and dysplastic nevus form two joint sets, it would be reasonable to study soft membership functions to categorize the members of classification. Another interesting aspect that has not been investigated is considering the class melanoma as a few subsets of different morphology.

Last but not least, probing different bases of the wavelet decomposition should explain the role of the basis for optimal selection of the feature set.

Acknowledgments

This work has been supported by the Polish Ministry of Science and Higher Education research grant N N518 419038 "Application of Computational Intelligence Techniques in Analysis of Images for Computer-Assisted Melanoma Diagnosis".