Differential diagnosis of squamous cell carcinoma in situ using skin histopathological images

Computers in Biology and Medicine 70 (2016) 23–39

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Differential diagnosis of squamous cell carcinoma in situ using skin histopathological images Navid Noroozi n, Ali Zakerolhosseini Department of Computer Engineering and Science, Shahid Beheshti University, Tehran, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 29 September 2015 Accepted 29 December 2015

Differential diagnosis of squamous cell carcinoma in situ is of great importance for prognosis and decision making in the disease treatment procedure. Currently, differential diagnosis is done by pathologists based on examination of the histopathological slides under the microscope, which is time consuming and prone to inter and intra observer variability. In this paper, we have proposed an automated method for differential diagnosis of SCC in situ from actinic keratosis, which is known to be a precursor of squamous cell carcinoma. The process begins with epidermis segmentation and cornified layer removal. Then, epidermis axis is specified using the paths in its skeleton and the granular layer is removed via connected components analysis. Finally, diagnosis is done based on the classification result of intensity profiles extracted from lines perpendicular to the epidermis axis. The results of the study are in agreement with the gold standards provided by expert pathologists. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Skin cancer Cutaneous squamous cell carcinoma Actinic keratosis Differential diagnosis Classification

1. Introduction Cutaneous squamous cell carcinoma (SCC) is the second most common form of skin cancer, constituting about 20% of all reported cases of non-melanoma skin malignancies [1,2]. It arises from uncontrolled growth of squamous cells in the upper layer of the skin tissue and has the potential to metastasize to other organs of the body if left untreated. Squamous cell carcinoma in situ, also called intra-epidermal SCC or Bowen disease, is the first stage of SCC. “In situ” means that the malignant cells are still only in the epidermis and have not invaded deeper into the dermis [3]. Differential diagnosis is the procedure of distinguishing a particular disease from others that present similar clinical signs, e.g. differentiating between prostate cancer and benign prostatic hyperplasia. In the case of SCC in situ, there are number of skin lesions that should be considered in differential diagnosis, among which, actinic keratosis and basal cell carcinoma (BCC) are more remarkable [3]. Differential diagnosis of SCC in situ versus actinic keratosis is not a trivial task as both share similar clinical and histological features. Currently, differential diagnosis is done by pathologists based on examination of the histopathological slides under the microscope, which is laborious, time consuming and prone to inter and intra observer variability. Moreover, there are relatively few expert pathologists against the large number of tissue samples to n Correspondence to: Department of Electrical and Computer Engineering, Shahid Beheshti University, Velenjak, Tehran, Iran. E-mail addresses: [email protected], [email protected] (N. Noroozi).

http://dx.doi.org/10.1016/j.compbiomed.2015.12.024 0010-4825/& 2016 Elsevier Ltd. All rights reserved.

be investigated. So, automated systems for analyzing the histopathological slides are desirable. Actinic keratosis is commonly believed to be a precursor of squamous cell carcinoma in situ and is referred to as incipient SCC by some authors [4], because it can progress to invasive squamous cell carcinoma if left untreated. Hence, in this work, the problem of distinguishing SCC in situ from actinic keratosis is addressed. Skin is composed of three primary layers: epidermis, dermis and hypodermis (subcutaneous adipose layer). Epidermis is divided into three layers: Malpighian (basal and squamous) layer, granular layer and cornified layer. Fig. 1a shows a histopathological image of normal skin in which the two main layers, epidermis and dermis are indicated. Cells existing in basal and squamous layer are called keratinocyte. In pathology, SCC in situ is recognized by the presence of dysplastic keratinocytes which cover the full-thickness of epidermis and in advanced stages, invade the dermis [3]. Dysplastic keratinocytes are immature cells with a relatively high nucleus to cytoplasm (N:C) ratio, the ratio which indicates the maturity of a cell [5]. As a cell matures, this ratio generally decreases. Unlike healthy cells, cancerous cells reproduce very quickly and do not have a chance to mature. Other pathological cues to rule out SCC in situ are cells polymorphism, hyperkeratosis (thickened cornified layer), parakeratosis (presence of keratinocytes nucleus in cornified layer) and increase in epidermis layer thickness. Actinic keratosis is also characterized by the above mentioned cues specially dysplastic keratinocytes in epidermis layer, but it differs from SCC in situ in the way that abnormal cells present mainly in lower third of the epidermis [3,6].

24

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 1. Skin histopathological images of (a) normal case, (b) actinic keratosis, (c) SCC in situ, and (d) borderline case. White brackets in (b)–(d) indicate the areas covered by atypical keratinocytes. (e)–(g) Magnified views of vertical strips cropped from epidermis area in (a)–(c). Yellow arrows in (f) and (g) show cluttered cells with poorly defined membranes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 1b and c shows, respectively, skin histopathological images of actinic keratosis and SCC in situ. In contrast to the normal case of Fig. 1a in which immature cells just cover a relatively thin layer at the bottom of the epidermis (basal layer), in actinic keratosis and SCC in situ cases, they cover the lower one-third and full thickness of the epidermis, respectively. Note also the atypia which presents among keratinocytes in the last two cases, in the sense that they are more polymorphic and of different sizes. White bracket indicates the area with atypical keratinocytes. It should be mentioned that Fig. 1b and c are representative samples of actinic keratosis and SCC in situ which exhibit almost distinct appearances in texture of the epidermis. However, there are borderline cases of actinic keratosis in which dysplastic keratinocytes have covered, to some extent, the middle third part of the epidermis as well. Fig. 1d shows an example of this kind. This case is labeled as actinic keratosis by expert pathologists. Although parakeratosis is usually less pronounced in actinic keratosis, but there are cases in which it is relatively more severe, resembling the appearance of SCC in situ. Consequently, parakeratosis cannot be used alone as a cue for differential diagnosis. Also note that from the images in Fig. 1, it can be observed that the cells in epidermis not only are overlapping and non homogenous, but also they are too cluttered to enable simple segmentation of individual cells, especially in the SCC in situ case. The organization of this paper is as follows. A brief review of the works in the field of histopathological images analysis is presented in Section 2. The proposed technique is described in Section 3; results are demonstrated and discussed in Section 4, followed by conclusion in Section 5.

2. Related works Histopathological images processing is a popular research field the importance of which has rapidly grown over the years due to increasing incidence of malignant diseases worldwide. There are many works in the literature which have covered the processing of histopathological images from various tissues. Cervical tissue is analyzed in [7–10] to diagnose and grade cervical intraepithelial neoplasia (CIN). The automated method proposed in [7] consists of two stages. First, the input image is segmented using a block based method. Statistical features such as gray level co-occurrence descriptors and entropy are calculated for each block of the image and support vector machine is employed to classify blocks into background, stroma, squamous epithelium, columnar epithelium and red blood cells. Then, cells nuclei in the squamous epithelium are segmented by incremental thresholding. Finally, another set of features including nuclei size, number of nuclei, average edge length and mean area of the Delaney triangulation are obtained from image blocks located along the perpendicular orientation to the main axis of segmented epithelium. Resulting feature vectors are classified by a support vector machine to grade the tissue. Authors in [8] used Gabor filters with six frequencies and four orientations to segment histopathological images as background, stroma, normal and abnormal cells. Then, classified samples into CIN 1, CIN 2, CIN 3 and malignant according to the ratio of normal and abnormal cells. The method proposed in [9] employed marked watershed algorithm for segmentation of cells nuclei in cervix histopathological images. Descriptors obtained from Delaunay triangulation, such as area of the triangles, mean and entropy of the nodes degree are used to discriminate and grade CIN. Authors in [10] utilized four different segmentation techniques, each of which uses 18 combinations of parameters to obtain 72 nucleus segmentation results per input cervical tissue image and then, mean nuclear volume is used as a feature to perform a cancer/not cancer classification, considering the fact

25

that cancerous cells have larger nucleus than non-cancerous cells. The main challenge with these strategies is nuclei segmentation, as all of features used in diagnosis and grading stage depend on the accuracy of this task. Prostate samples are analyzed in [11–13] to accomplish Gleason grading and normal/abnormal classification, using fractal dimension, fractal code and color histogram features. Higher order spectra and local binary pattern, along with fuzzy classifier are employed in [14] for oral cancer detection and grading. The problem of breast cancer diagnosis using La*b* color space and Gaussian mixture model is addressed in [15] and an overview of methods for analyzing breast histopathological images is presented in [16]. There are also some works which have covered colon and lung cancer detection [17–19]. There are relatively few works in the literature which cover skin histopathological images processing. The problem of melanocyte detection in epidermis area has been addressed for the first time in [20]. The authors employed morphological operation and adaptive thresholding to segment cells nuclei and then recognized melanocytes by analyzing the gradient magnitude and direction of pixels located on radial lines originating from centroid of cells nuclei. Another approach is presented in [21] for segmentation of melanocytes which is based on a double ellipse descriptor. After performing mean shift algorithm to decompose the image into more homogenous regions, nucleus candidates are selected using a recursive split and merge based algorithm. In the next stage, an ellipse is fitted to these candidate regions to filter out false positives based on an ellipticity parameter. At last, another ellipse with larger axes is built around the first ellipse and melanocytes nuclei are recognized by analyzing the histogram of intensity values of surrounded pixels. Skin epidermis segmentation problem is addressed in [22–25] using Otsu’s algorithm, shape features, statistical analysis and thickness measurement. Based on the work in [20] for melanocyte detection and the approach in [22] for epidermis segmentation, authors in [26] presented a method for melanoma diagnosis using histopathological images. In [27], adaptive thresholding, ellipse descriptor analysis and marked watershed algorithm are used to segment cells nuclei in skin images. A voting procedure based on the image gradient direction is applied to select the seeds needed for marked watershed algorithm. The problem of melanocytic tumor depth measurement, which is necessary for melanoma grading, was first addressed in [24] and further investigated and improved in [28]. There are many works in the field of medical image processing regarding differential diagnosis of SCC using clinical and dermoscopy images (for example, see [29–32]), but to the best of our knowledge, differential diagnosis using histopathological images has not yet been addressed. Since histopathological image analysis is the gold standard for diagnosing and grading SCC and other skin tissue malignancies, in this paper we have proposed a method for automating this procedure.

3. Materials and methods 3.1. Materials For conducting this study, 30 histopathological images were obtained from skin tissue samples of unknown patients with the age between 40 and 60, sixteen images of SCC in situ and fourteen images of actinic keratosis. After sectioning the tissue samples at a thickness of about 5 mm and staining using Hematoxylin and Eosin (H&E), high resolution images were digitally captured by Aperio ScanScope slide scanner from histopathology slides at 60
magnifications with pixel size of 0.28 mm. From these high resolution images, regions of interest (regions which contain epidermis area)

26

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

are cropped by pathologist and stored in 24-bit BMP format. All images were evaluated by expert pathologists to rule out SCC in situ or actinic keratosis and labeled accordingly. The decisions made by pathologists were treated as ground truths in all experiments. From the total of 30 images, 6 images (3 images from each tissue type) were used for training procedures. 3.2. Proposed method The proposed approach is based on extracting the intensity profile of straight lines perpendicular to the epidermis axis and capturing the pattern of the profiles by analyzing their frequency distribution, as follows. As an illustration, the method will be applied to the sample in Fig. 1b in a step by step manner. 3.2.1. Extraction of granular layer boundary In the first stage, epidermis region is segmented based on the approach in [24] and then, granular layer boundary is extracted by the method we presented in our previous work [28]. After cornified layer segmentation and performing morphological opening and closing operations by a disk shaped structuring element, boundary pixels of the epidermis (with removed cornified layer) which have no neighbor from dermis region are segmented and the largest connected component is marked as the granular layer boundary. Granular layer boundary is extracted two times, which results in two distinct shapes for the boundary. The first one is obtained using the aforementioned method and the other is obtained in the same manner, except that in morphological opening and closing step, the disk shaped structuring element is replaced by a square

shaped one to preserve the “blockiness” form of the boundary. The former is used in epidermis axis extraction while the latter is suitable for granular layer segmentation stage, as will be discussed in the following sections. Fig. 2a shows epidermis segmentation result for the sample in Fig. 1b after cornified layer removal. Boundary of the granular layer obtained using the method in [28] (with disk shaped structuring element) is shown in blue. 3.2.2. Specifying the epidermis axis Authors in [25] defined epidermis axis as the longest path between endpoints of its skeleton. As will be shown shortly, this is not accurate, especially in cases in which epidermis thickness is large in some parts or there are hair follicles appended to the epidermis. So, here we define epidermis axis as a sequence of image pixels which are located approximately at the central part of the epidermis with an orientation similar to that of the granular layer boundary. The basics of how to extract the epidermis axis are illustrated in Fig. 2b–d. At first, skeleton of the segmented epidermis (with cornified layer removed) is extracted using morphological operations. Note that before skeleton extraction, a simple preprocessing is performed on the image: Boundary of the segmented epidermis, except the part corresponding to granular layer boundary, is traced and small square regions of size 5*5, centered on pixels on the boundary, with a predefined distance dsq from each other are removed (set to black) to create porosities on the boundary. This is done in order to increase the number of branches in the skeleton (an object with smooth boundary has fewer number of branches in

Fig. 2. (a) Result of epidermis segmentation for the image in Fig. 1b, after cornified layer removal. Shown in blue is the granular layer boundary (thickened for illustration purpose). (b) Skeleton of the segmented epidermis (shown in green). (c) Graph of the skeleton (shown red) along with the nearest points on the granular layer boundary, p0i and p0j , to points pi and pj on the skeleton. Orientation of the line segment p0i p0j (shown yellow) is considered as the granular layer boundary orientation near the edge (pi,pj). (d) Epidermis axis, shown dark blue. The path between points A and B is the longest path between end points of the epidermis skeleton, which is considered as epidermis axis by the method in [25]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

its skeleton). In this way, the chance of finding a more suitable axis will increase.
 Let P ¼ p1 ; ::::; pN be the set of all endpoints and branch points of the skeleton, then a graph G ¼ ðV; EÞ is built in which V is the set of nodes in the graph, such that V ¼ P (end points and branch points of the skeleton are considered as the nodes of the graph). E, the set of graph edges, is defined as n o ð1Þ E ¼ ðpi ; pj Þj pi ; pj A P 4 contigðpi ; pj Þ ¼ True where contigðpi ; pj Þ is an indicator such that  False ( kn j kn A Br 4 kn A pathðpi ; pj Þ contigðpi ; pj Þ ¼ True otherwise:

ð2Þ

in which Br is the set of all branch points and pathðpi ; pj Þ is the set of all pixels located in the path between pi and pj in the skeleton. In other words, contigðpi ; pj Þ is true if the points pi and pj in the skeleton are contiguous, in the sense that there is not any branch point along the path connecting them. Fig. 2b and c shows the skeleton obtained from the segmented epidermis in Fig. 2a and the corresponding graph, respectively. For each edge ðpi ; pj Þ in the graph, we define an orientation which is simply the orientation of the straight line connecting points pi and pj in the skeleton. Now, since we are looking for sequence of pixels which goes through the epidermis length and in parallel with the granular layer boundary, we assign a cost to each edge ðpi ; pj Þ of the graph G calculated as   L½ðpi ; pj Þ L½ðpi ; pj Þ ΔO½ðpi ; pj Þ þ w2 U U C½ðpi ; pj Þ ¼  w1 U ML ML 90    L½ðpi ; pj Þ ΔO½ðpi ; pj Þ ¼ U  w1 þw2 U ð3Þ ML 90 where w1 and w2 are weighting constants (adjusting their values is discussed in Section 4.3), L½ðpi ; pj Þ is the length of the path in the skeleton connecting points pi and pj , ML is the maximum value of L½ðpi ; pj Þ in the skeleton (used for normalizing path length and prevent large values dominate the result) and ΔO½ðpi ; pj Þ=90 is the normalized difference between the orientation of the edge ðpi ; pj Þ and that of the granular layer boundary limited to the region near the points pi and pj , calculated using the following procedure. Let p0i and p0j denote the nearest points on the granular layer boundary to nodes pi and pj of edge ðpi ; pj Þ in the graph, as is

27

shown in Fig. 2c. ΔO½ðpi ; pj Þ is defined as the difference between the orientation of line segment p0i p0j and the orientation of the edge ðpi ; pj Þ with the range of ½0; 90. For the special case of p0i ¼ p0j , k points before and after p0i on the boundary are considered to establish line segments (in our implementations, we set k ¼ 3, i.e. we considered three points before and after p0i , with a distance of 32 pixels from each other, so that three different values for the orientation of boundary is obtained. Finally, the average of these three values is considered as the boundary orientation near the edge P i P j ). Finally, among the paths between each pair of leaf nodes in the graph, the path with lowest cost is selected and all the nodes in this path are marked. The path in the skeleton which goes through these marked nodes is chosen as epidermis axis. Note that in the second term of Eq. (3) (in the first line), multiplying byL½ðpi ; pj Þ=ML is to prevent an edge with small Lð:Þ and large ΔOð:Þ dominating the total cost of a path. The edge ðpr ; pc Þin Fig. 2c is an example of this case. Fig. 2d shows the location of epidermis axis for the skin histopathological image of Fig. 1b obtained using the aforementioned procedure. In this case, the path between points A and B is the longest path between end points of the epidermis skeleton, which is considered as epidermis axis by the method in [25]. Another sample along with the corresponding epidermis axis is shown in Fig. 3. Note that, part of the epidermis has protruded to the dermis and surrounded a skin appendage (in this case, a hair follicle). Obviously, the path between points A and B is the longest path between end points. It is evident from the images in Figs. 2d and 3b that the longest path in the skeleton is not always a suitable choice for epidermis axis. Our method selects the best path by penalizing the paths which are too short or deviate from the orientation of granular layer boundary. 3.2.3. Granular layer removal As mentioned at the beginning of this section, differential diagnosis is based on capturing the pattern of the intensity profiles extracted from straight lines perpendicular to the epidermis axis. Obviously, some of the perpendicular lines coincide with areas containing granular layer. Since we are going to analyze the frequency distribution of each line's profile, image blocks containing granular layer areas should be recognized and removed, as they

Fig. 3. (a) Original skin histopathological image with an appendage (b) Result of epidermis segmentation after cornified layer removal along with the skeleton and axis of the epidermis obtained by the proposed method, shown in green and dark blue, respectively. Light blue indicates boundary of the granular layer (thickened for illustration purpose). The path between points A and B is the longest path between end points of the epidermis skeleton, which is considered as epidermis axis by the method in [25]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

28

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 4. (a) Result of epidermis segmentation after cornified layer removal, processed with a square shaped structuring element. (b)–(e) Magnified view of image blocks 1–4, marked with yellow boxes in (a). (f)–(i) Binary version of images in (b)–(e), with values of (MA,NCC) equal to (112,15), (506,3), (314,6) and (176,2), respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

may disturb the frequency distribution of the extracted profiles. To this aim, we propose to use the properties of connected components in binarized version of image blocks which will be shown to have a good discriminating power in classification of granular layer versus other parts of a skin histopathological image. Starting from boundary blocks which have a neighbor from cornified layer, each block bl is converted to gray scale and binarized using active contours model [33]. In this stage, pixels with gray level values in the range ½min; min þ a
ðmax minÞ are considered as seeds, where min and max are the minimum and maximum gray level values in bl and a is a predefined parameter (0.1 in our implementations). Fig. 4a shows the result of epidermis segmentation after cornified layer removal, as in Fig. 2a, but processed with a square shaped structuring element instead of a disk shaped one to preserve the “blockiness” form of the boundary (otherwise it would not be a trivial task to establish image blocks in boundary areas). Four image blocks of size 64*64 are magnified, one from granular layer and the others from squamous layer. Binary versions of these blocks are shown, too. Considering blocks 1 and 2, it is obvious that there are many connected components with small size in one corresponding to granular layer (because of “granulated” pattern of granular layer), and a few number with relatively large size in the other. Hence, two features, namely as MA (rounded mean of connected components area) and NCC (number of connected components) are used for discrimination. Note that employing just NCC is not sufficient, because due to staining imperfections, there are cells with some fragments in their wall, which results in slightly increasing the number of connected components in binary version, as is the case for block 3 (note that in general, segmentation of cell’s wall is not a trivial task), but in such cases, the mean of the area is still relatively large. Also, MA alone is not sufficient, for example, block 4 contains a small part of a cell. However, this block doesn’t belong to granular layer. So number of connected components should be taken into consideration when MA is small. Note that for the almost smooth areas of squamous layer, the foregoing procedure results in a few numbers of connected components with large area, so that they will be classified as “notgranular”. Fig. 5 shows the distribution of about 75 blocks from granular layer and 75 blocks from squamous layer in MA-NCC

Fig. 5. Distribution of about 75 blocks from granular layer and 75 blocks from squamous layer in MA-NCC space, along with the decision boundary obtained using support vector machine.

space, along with the decision boundary obtained using support vector machine, which is built as follows. The decision function in a support vector machine is given by [34] gðXÞ ¼ W T U ΦðXÞ þb

ð4Þ

where ΦðXÞ is a linear or non-linear mapping from input feature space into hidden feature space, W and b are optimal weight vector and bias, determined from the training set. To classify granular and not-granular blocks, we employed a linear innerproduct kernel (first order polynomial kernel), defined as KðX; X i Þ ¼ X T U X i

ð5Þ

where X i is a sample’s feature vector from training data. By using a linear kernel and fine tuning the parameter C of the classifier (which controls the tradeoff between complexity of the classifier and number of non-separable samples), we minimized the risk of over fitting the classifier to the training data. If C is too large, we have a high penalty for non-separable points, the classifier may store many support vectors and over fit. If it is too small, too

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

29

Fig. 6. Result of granular layer removal performed on the image in Fig. 1b. Red arrows show blocks of the image that cornified layer segmentation algorithm failed to detect (they are remaining parts of cornified layer). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Five perpendicular lines to the epidermis axis (shown in yellow) used to extract intensity profiles for subsequent sample-frequency distribution analysis, superimposed on gray scale version of the image in Fig. 1b, processed with anisotropic diffusion filter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

simple solutions will be found which causes the classifier to under fit [34]. To adjust the parameter C, a train-test procedure was performed, in which, training set consisted of 150 image blocks from granular layer and Malpighian layer (equal proportions from each class), with the same number of blocks for testing. Blocks for train and test were obtained from 6 images in train group and 24 images in test group, respectively. As mentioned previously, granular layer removal is started from boundary blocks. Once a block is detected as granular, the neighbor blocks are also investigated to be classified, until no further block is detected as granular layer. Result of granular layer removal performed on the image in Fig. 1b is shown in Fig. 6. 3.2.4. Intensity profiles extraction and analysis After granular layer removal, the input image is converted to gray scale and processed using anisotropic diffusion filter to suppress the effect of noise while preserving the edges [35]. Then, perpendicular lines to the epidermis axis with distance d from each other are established and the values of image pixels along these lines are sampled to create Nip intensity profiles. Fig. 7 shows five perpendicular lines to the epidermis axis superimposed on the image in Fig. 1b and the corresponding intensity profiles are plotted in Fig. 8 (note that, the beginning of the intensity profiles are at the bottom of the Malpighian layer). Similarly, Fig. 9 depicts plots, obtained in the same manner, from the SCC in situ case of Fig. 1c. Obviously, intensity profile which is extracted from actinic keratosis sample can be thought as a finite-duration signal which has many large amplitude oscillations at the beginnings of its duration with decrease in their number as the signal proceeds. On

the other hand, in a profile extracted from SCC in situ sample, number of large amplitude oscillations remains relatively constant over whole of the signal duration. In fact, the former contains high frequency components at the beginning and low frequency components at the end, while the high frequency content of the latter doesn’t change significantly over the signal duration. So, it is a natural thought that “time–frequency distribution” concept, which is referred to as “sample1  frequency distribution” in the context of our work, can be appropriate for analyzing the extracted profiles. Originally devised for non-stationary signals analysis, time– frequency distribution (abbreviated as TFD) makes it possible to study a signal in both the time and frequency domains simultaneously. Specially, TFD for a signal xðtÞ, is a joint function of time and frequency, stated as Wðf ; tÞ, which describes the contribution of sinusoid with frequency f in construction of the signal at time t. There are several approaches to compute TFD of a signal, among which, the Cohen’s class formulations [36] and short-time Fourier transform [37] are the most widely used techniques. Here, we employ short-time Fourier transform which were found empirically to give more satisfactory results: Z þ1 xðτÞ U hðτ  tÞ U e  i2π f τ dτ ð6Þ Wðf ; tÞ ¼ 1

1 The meaning of the word “sample” used throughout of this section regarding “sample-frequency distribution” should not be confused with its meaning as specimens in a dataset. By “sample”, here we mean distinct values of the independent variable in a signal.

30

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 8. Signals obtained from intensity profile of perpendicular lines 1–5 shown in Fig. 7.

where hðtÞ is a window function which is nonzero for only a short period of time. TFD is an ideal tool to deal with a situation where the frequency composition of a signal changes over

time and has been widely used in the field of electroencephalogram signal processing and epilepsy diagnosis [38,39]. In the context of our work, “time” variable is replaced

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

31

Fig. 9. Signals obtained from intensity profiles of the SCC in situ sample in Fig. 1c.

by “spatial location” or “x” variable, so we use the term “sample-frequency distribution” (SFD) instead of “time–frequency distribution”.

As an illustration of how SFD analysis could be used to differentiate between intensity profiles extracted from SCC in situ and actinic keratosis images, the SFD for the intensity profiles shown in

32

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Figs. 8 and 9 are calculated and then, the values of instantaneous frequencies at each spatial location x are obtained as fP ¼π

IFðxÞ ¼

f ¼0

given by [40] PðC i j XÞ ¼

f U Wðf ; xÞ

fP ¼π f ¼0

ð7Þ Wðf ; xÞ

According to the above equation, instantaneous frequency at location x is a weighted average of all constituent frequencies at location x. Values of instantaneous frequencies for the intensity profiles of Figs. 8 and 9 are plotted in Fig. 10. To better visualize the overall shape of the curves, a second order polynomial is fitted to each one (shown in red). Note that a value of w ¼ 100 for width of the window hðxÞ was used and SFD was calculated for x A ½w=2; L  w=2, where L is the length of the intensity profile (it does not mean that first and last w=2 samples of the signal are ignored, as they do have contribution in frequency content of their nearby samples in SFD). As can be seen, values of instantaneous frequency for the actinic keratosis sample can be thought as a decreasing function (considering its overall shape); on the other hand, there are no significant changes in the values obtained from SCC in situ case. Consequently, the first feature used to discriminate between these two types of intensity profile is defined as 0P 1 x UIFðxÞ wA 1 x @ P U ð8Þ  NCoG ¼ 2 N IFðxÞ x

where N is the number of spatial locations at which instantaneous frequencies are calculated (length of the IFðxÞ’s support). The first term in the brackets in Eq. (8) is simply recognized as the center of gravity of IFðxÞ, from which w=2 is subtracted and multiplied by 1=N to normalize it to the range ð0; 1Þ. The value of normalized center of gravity for instantaneous frequencies shown in left column of Fig. 10 are, from top to bottom, 0.4278, 0.3998, 0.3812, 0.4214 and 0.4375, while for those shown in right column, this feature has the values 0.5221, 0.5277, 0.5147, 0.5069 and 0.5049, respectively. The second feature used in this stage is the average of instantaneous frequency values in the second half of IFðxÞ, as the intensity profiles from actinic keratosis cases are expected to have smaller values of instantaneous frequencies in this interval: AIF ¼ AverageðIFðxÞÞ; x A ½L=2; L  w=2

ð9Þ

The value of average instantaneous frequency for the plots shown in left column of Fig. 10 are, from top to bottom, 0.1145, 0.1991, 0.1271, 0.2003 and 0.1540, while for those shown in right column, this feature has values 0.2875, 0.2544, 0.2887, 0.3458 and 0.3323. Fig. 11 shows the estimated probability density functions for distribution of about 900 intensity profiles obtained from 3 actinic keratosis images and 3 SCC in situ images, in (NCoG-AIF) feature space. Note that features are considered as continuous variable to estimate the PDFs, so values greater than 1 for PDF are expected. 3.2.5. Classification In this stage, each of the intensity profiles extracted from segmented epidermis is classified as either “actinic keratosis” (class “0”) or “SCC in situ” (class “1”), using a Bayesian classifier, based on the features described in previous section. According to Baye’s rule, an intensity profile is classified to a specific class that has the highest posteriori probability among all the others. If the prior probability PðC i Þ and class conditional probability density function PðX j C i Þ for feature vector X is known, then posteriori probability is

PðX j C i Þ U PðC i Þ NC P PðX j C j Þ UPðC j Þ

ð10Þ

j¼1

where NC is the number of classes. In this work, we assumed that prior probabilities for the two classes, “actinic keratosis” and “SCC in situ” are equal, i.e. PðC 1 Þ ¼ PðC 2 Þ. The class conditional probability density functions, PðX j C 1 Þ and PðX j C 2 Þ are modeled as a linear mixture of bivariate Gaussian probability density functions, with maximum of four Gaussian kernels to avoid over fitting the classifier to training data. About 900 intensity profiles from 6 images (3 images from each class) were used to train the Bayesian classifier. Fig. 12a shows the result of this classification for 187 intensity profiles obtained from the image in Fig. 1b, in the form of a stem plot. Horizontal axis indicates the profile number and the vertical axis corresponds to the possible classes: “0” or “1”. Note that there are intensity profiles that are classified as “1”, because there are occasionally concentrations of cells at the upper parts of the epidermis, even in an actinic keratosis case. But a sample is recognized as SCC in situ if there are many contiguous intensity profiles which are classified as “1”. Hence, the following procedure is applied: First, the class label of each single intensity profile with label “0” which is located between two members from class “1” is set to “1” (as they are more likely to be misclassified) and then, runs of “1”s with length larger than a predefined threshold TL are indicative of regions corresponding to SCC tumor. We have found empirically that the value of TL ¼25 is appropriate, which gives the best classification results. A sample is classified as SCC in situ if contains at least one tumor region. Note that special cases in which, all of the intensity profiles are alternatively labeled “0” and “1” (before setting “0”s to “1”s) will be rejected, as they cannot be classified with high confidence. Result of performing this procedure for the image in Fig. 1b is shown in Fig. 12b, with the length of longest run of “1”s equal to 6, so no tumor is detected. On the other hand, processing the intensity profiles extracted from the SCC in situ image of Fig. 1c, yields 4 runs of length 32, 89, 35 and 33.

4. Results and discussion Performance of the proposed method was evaluated in two parts: classification of samples and tumor area segmentation. 4.1. Classification As this is the first work in which differential diagnosis of SCC in situ from actinic keratosis is addressed, performance of the proposed method is compared to the gold standard provided by expert pathologists. As stated previously, about 900 intensity profiles obtained from 6 samples (3 samples from each class) were used to train a Bayesian classifier and the remaining samples were employed in test phase. The resulting confusion matrix is presented in Table 1. Note that the borderline case shown in Fig. 1d was classified correctly as actinic keratosis by our system, with the length of longest run of “1”s equal to 14. If we adopt the convention that “SCC in situ” corresponds to “positive” class, the values for sensitivity, specificity, precision and accuracy can be computed from confusion matrix, as presented in Table 2. To further illustrate the proposed classification method, results of applying the method on the SCC in situ case of Fig. 1c are depicted in Fig. 13.

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

33

Fig. 10. Instantaneous frequencies of the signals depicted in Fig. 8 (left column) and Fig. 9 (right column), shown in blue, along with the corresponding second order polynomials, shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

34

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

4.2. Tumor region segmentation

accuracy ¼

As stated previously, all runs of “1”s in the intensity profiles classification result longer than the threshold TL are indicative of SCC tumor. Consequently, the region in the epidermis enclosed between the perpendicular lines corresponding to the first and last intensity profiles of such runs is considered as a tumor region. To evaluate this stage, the following metrics are used: sensitivity ¼

precision ¼

j GT \ SEGj j GT j

j GT \ SEGj j SEGj

specif icity ¼

j GT \ SEG j j GT j

ð11Þ

ð12Þ

ð13Þ

j GT \ SEG j þ j GT \ SEGj

ð14Þ

j GT j þ j GT j

where GT is the ground truth (set of pixels in the regions diagnosed as tumor by pathologist) and SEG is the set of pixels segmented as tumor region by the method. The results obtained using the above metrics are shown in Table 3. Note that the presented results are averaged among all SCC in situ samples, except those two recognized as actinic keratosis (false negatives), because it is obvious that for these two samples, sensitivity and specificity are equal to 0% and 100%, respectively, with undefined precision (no tumor region is detected). Fig. 14a depicts a skin histopathological image in which SCC tumor is located in a limited part of the epidermis (shown with red arrow). Tumor region segmented manually by the pathologist and the result of segmentation by the proposed method are shown in Fig. 14b and c, respectively. Note that the relatively low sensitivity of tumor region segmentation method is due to the gaps between runs of “1”s in the intensity profile classification result for some samples. On the other hand, high value of specificity indicates small amount of over segmentation by the algorithm. Table 1 Confusion matrix of the proposed method for differential diagnosis of SCC in situ from actinic keratosis. Actual class

Predicted

SCC in situ Actinic keratosis

SCC in situ

Actinic keratosis

11 2

2 9

Table 2 Quantitative results for classification stage of the proposed method. Fig. 11. Estimated probability density functions for distribution of about 900 intensity profiles in (NCoG-AIF) space, obtained using Gaussian kernels. The red mesh corresponds to SCC in situ samples and the blue one corresponds to actinic keratosis samples. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Sensitivity

Precision

Specificity

Accuracy

84.6%

84.6%

81.8%

83.3%

Fig. 12. (a) Result of classification of intensity profiles extracted from the image in Fig. 1b, into class “0” (actinic keratosis) or “1” (SCC in situ). (b) Result obtained by setting to “1” the class label of intensity profiles with label “0” which are between two members from class “1”.

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

35

Fig. 13. (a) Boundary of granular layer for the SCC in situ image in Fig. 1c along with the epidermis skeleton and epidermis axis. (b) Result of intensity profiles classification. Table 3 Quantitative results of tumor region segmentation stage. Sensitivity

Precision

Specificity

Accuracy

81.2%

88.1%

97.6%

92.1%

4.3. Notes on parameters adjustment The value of three parameters that are used in epidermis axis specification stage, i.e. w1 and w2 in Eq. (3), as well as dsq (distance of porosities on the boundary of segmented Malpighian layer, in Section 3.2.2) are determined by means of a grid search on the parameter space, such that w1 A f0; 0:1; 0:2; 0:3; :::; 1:0g, dsq A f20; 30; :::; 250g and w2 ¼ 1  w1 . For this purpose, an expert pathologist was asked to manually segment the epidermis, exclude cornified layer and specify an axis for epidermis. The objective function to be minimized in grid search is the mean absolute distance, defined as MADðaxisGT ; axisSK Þ ¼

 N  1 X SK U minðdistðaxisGT i ; axisj ÞÞ N i¼1 j

ð15Þ

where axisGT is the axis specified by pathologist, axisSK is the axis axisSK obtained from skeleton automatically, axisGT i and j are, respectively, ith pixel on the ground truth axis and jth pixel on the axis obtained automatically, N is total number of pixels in the

ground truth axis and distðp; qÞ is the Euclidean distance between pixels p and q: Using foregoing procedure, the values of parameters w1 , w2 and dsq are determined to be 0:4, 0:6 and 20, respectively. This is intuitively reasonable that, as dsq decreases, the number of branches in the skeleton and thus, the chance for finding a more accurate axis increases. But on the other hand, with a large number of branches in the skeleton, computational complexity of the epidermis axis specification stage increases significantly. So, to lessen the computational burden, we increased the value of dsq past 20 in a step by step manner (with increments of 10) until the value of MAD grew by 10% of its minimum value which was obtained previously. Finally, a value of dsq ¼ 140 found to be appropriate, resulting in a good tradeoff between accuracy and computational complexity. There are four other parameters used in the algorithm which are determined empirically: a (in Section 3.2.3 regarding seed selection), d (in Section 3.2.4 for intensity profile extraction), w (in Section 3.2.4 regarding short time Fourier window size) and T L (in Section 3.2.5 for classification stage). To evaluate the effect of changes in the value of this set of parameters, the algorithm is executed several times by varying the value of each parameter separately and accuracy of the classification stage is calculated. The results are shown in Table 4. Change in the value of parameter a alters the performance of classification by affecting that of the granular layer detection stage. As a increases, more (false) seeds with larger areas are introduced in the mask image, which are used by active contour algorithm to

36

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 14. Illustration of tumor region segmentation: (a) original skin histopathological image with a SCC tumor shown by arrow; (b) tumor region marked manually by the pathologist and (c) result obtained by the proposed method. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Table 4 Effect of different selections for the proposed algorithm's parameters on samples classification result. Parameter name a

Parameter values 0.1 0.2 0.3 0.4 0.5

Sensitivity (%)

Precision (%)

Specificity (%)

Accuracy (%)

84.6 84.6 84.6 84.6 84.6

84.6 84.6 78.5 78.5 73.3

81.8 81.8 72.7 72.7 63.6

83.3 83.3 79.1 79.1 75.0

d

1 10 19 28

84.6 84.6 84.6 69.2

84.6 84.6 84.6 81.8

81.8 81.8 81.8 81.8

83.3 83.3 83.3 75.0

w

40 60 80 100 120 140

46.1 69.2 76.9 84.6 84.6 84.6

60.0 90.0 83.3 84.6 84.6 78.5

63.6 90.9 81.8 81.8 81.8 72.7

54.1 79.1 79.1 83.3 83.3 79.1

TL

10 15 25 45 55 65

84.6 84.6 84.6 76.9 76.9 53.8

73.3 84.6 84.6 100 100 100

63.6 81.8 81.8 100 100 100

75.0 83.3 83.3 87.5 87.5 75.0

perform segmentation on image blocks. This causes some parts of the granulated structure of granular layer merge together and, consequently, leads to presence of relatively less number of connected components with larger areas in the segmentation result. So, some granular layer blocks are misclassified as “not-granular”.

These “missed” granular blocks disturb the frequency content of the intensity profiles, extracted from actinic keratosis samples, by increasing the contribution of high frequency sinusoids near the end point of the intensity profile and thus, increasing the probability of misclassifying the profile as SCC in situ. The rows in

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Table 4 corresponding to a ¼0.4 and a ¼0.5 confirms this reasoning, in which specificity degrades to 72.7% and 63.6%, respectively. Note that the sensitivity of the result remains intact. Intuitively, the best choice for parameter d is 1, since it provides the ability to capture the pattern of all parts of the Malpighian layer, throughout the main axis. As d increases past 1, width of vertical strips of Malpighian layer that are missed, increases as well, which may lead to failure in detection of some parts of SCC tumor. Although a choice of d ¼19 still causes no change in performance of the proposed algorithm on our dataset, however, for general applicability of the method, we chose d ¼10 pixels which results in a good tradeoff between accuracy and computational complexity. Perhaps the effect of changes in parameter w on the classification result is more remarkable. When the window size is too small, the short time Fourier transform fails to capture the frequency content of signal segments encompassed by the window accurately. Hence, some of the intensity profiles from actinic keratosis class are misclassified as SCC in situ and vice versa. Relatively low values of sensitivity and specificity in Table 4, in the row corresponding to w¼ 40, confirm this reasoning. On the other hand, large values for w, leads to failure in capturing the frequency content in spatial locations near the start point and end point of the signal, since sample frequency distribution is calculated for x A ½w=2; L  w=2 (L is the length of the intensity profile). This concept is illustrated in Fig. 15. Fig. 15a shows an intensity profile from a borderline case of actinic keratosis. Red arrow shows the range of spatial locations for which SFD is calculated when w¼ 140. Fig. 15b and c shows, respectively, plots of instantaneous frequencies for w¼100 and w ¼140. The former is classified as actinic keratosis while the latter is misclassified as SCC in situ. When w¼140, in spatial locations near x ¼ L  w=2, still some parts of the signal with large oscillations are encompassed by the window, so high frequency content of the signal in these locations increases slightly, which in turn, increases the probability of misclassification of the intensity profile, while this is not the case for lower values of w. Of course, this occurs only in border line cases

37

of actinic keratosis, which is confirmed by the results in row corresponding to w¼140 in Table 4. As can be seen, the value of specificity has dropped to 72.7%, since one borderline case of actinic keratosis is misclassified as SCC in situ. It should be emphasized that spatial locations after x ¼L  w/2, do have contribution in instantaneous frequency values, but their contribution is less pronounced with higher values of w. Effect of changes in parameter TL on classification performance can be seen in the last part of Table 4, which is further investigated by plotting the corresponding ROC curve in Fig. 16. The ROC curve is obtained by decreasing the value of parameter TL, from 128 (maximum value for the length of runs in the samples) to 0, in a step by step manner, and calculating the operating point of the diagnosis system, in terms of true positive rate and false positive rate. Relatively high value of area under the curve (AUC), which is 0.8462, indicates satisfactory performance of the proposed method.

Fig. 16. Receiver operating characteristic curve for the proposed diagnosis system.

Fig. 15. Effect of change in parameter w on the classification of intensity profiles. (a) An intensity profile from a borderline case of actinic keratosis. The red arrow shows the range of spatial locations for which SFD is calculated with w¼ 140 and the green one indicates spatial locations encompassed by a window of size w¼ 140, when calculating SFD for the last location in the range, i.e. x ¼330. (b) Plot of instantaneous frequencies with w¼ 100, for which (NCoG,MA) ¼ (0.4403,0.2013). (c) Plot of instantaneous frequencies with w¼140, for which (NCoG,MA) ¼(0.4852,0.2454). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

38

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

Fig. 17. (a) Part of a misclassified SCC in situ sample with blurring artifact, (b) intensity profile extracted from the red line superimposed on (a), and (c) corresponding instantaneous frequency plot. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.4. Execution time The proposed approach was implemented using MATLB R2014b and run on a PC (CPU Corei5-3.00 GHz, 4.00 GB RAM) with execution time of about 6.15 min for an image of size 6300*4500 pixels (using a not optimized code), which is relatively high. Cornified and granular layer segmentation stages were performed on the image with 60
magnification, while other parts of the algorithm, including epidermis axis construction, intensity profile extraction and classification were done using the down sampled version of the image with a factor of 2, to lessen the computational burden. The most computationally extensive part of the algorithm is cornified layer segmentation stage, which took about 4 min. 4.5. Discussion The most important aspect of our work is that we have employed minimum number of features in each stage of the algorithm, which not only eliminates the need to feature reduction approaches, but also enables the classifiers to be trained with a relatively small set of training samples. This is a valuable characteristic for medical systems, as the number of available samples for training the system is usually small. Moreover, differential

diagnosis is done without individual cell segmentation, which is a challenging task due to the cluttered arrangement of cells and fragments in their membrane. Note also that diagnosis based on texture classification would not be a trivial task because of variability of the cells shape and cytology, which results in different possible textures and necessity of larger number of features to discriminate between them. One source of error in our algorithm is the blurring artifact which presents in some of H&E histopathological images. Blurred regions in dermis layer have no negative effect on the algorithm’s result, but when blurring occurs in epidermis, it may disturb the frequency content of the intensity profiles and consequently, lead to misclassification of these profiles. Fig. 17a shows a part of a SCC in situ sample in which the upper region of the epidermis has been degraded by blurring artifact. Intensity profile extracted from the values of pixels along the red line superimposed on the image (after converting to gray scale and anisotropic diffusion filtering) and its instantaneous frequency plot are depicted in Fig. 17b and c, respectively. This profile is classified as actinic keratosis. According to the literature, histopathological images deblurring has not been well studied (for example, see [41]) and can be a field for research in future to improve the performance of medical systems which process such images.

N. Noroozi, A. Zakerolhosseini / Computers in Biology and Medicine 70 (2016) 23–39

5. Conclusion and future works In this paper, a method for differential diagnosis of squamous cell carcinoma in situ from actinic keratosis is proposed. After epidermis segmentation and cornified layer removal, epidermis axis is specified using the paths in its skeleton and granular layer is removed via connected components analysis. Finally, diagnosis is done based on the classification result of intensity profiles extracted from lines perpendicular to the epidermis axis. Evaluation results have shown that our method is globally efficient, but some improvement is needed to increase sensitivity in tumor area detection stage, as well as reducing computational cost of the algorithm. The presented system can be used as a part of a multiexpert system which we are going to propose in future, for automated diagnosis of non-pigmented skin lesions. However, it should be mentioned that computerized methods for assessment of histopathological images only provide the pathologists a second consultative opinion and their prediction can be used as support of diagnostic decision. They never eliminate the role of the pathologists and their administration in judgment and decision making.

Conflict of interest statement None declared.

Acknowledgment The authors would like to thank Dr. Maryam Almassi, M.D., A.P. C.P, from Vanak Pathobiology Laboratory (Tehran, Iran) for insightful advices and helping us in medical aspects of this research.

References [1] S. Kalia, M.L. Haiducu, The burden of skin disease in United States and Canada, Dermatol. Clin. J. 30 (1) (2012) 5–18. [2] B.K. Armstrong, A. Kricker, Skin cancer, Dermatol. Clin. J. 13 (1995) 83–94. [3] D.E. Elder, R. Elenitsas, B.L. Johnson, Lever's Histopathology of the Skin, second edition, 2007. [4] C.J. Cockerell, Histopathology of incipient intra epidermal squamous cell carcinoma (“actinic keratosis”), J. Am. Acad. Dermatol. 42 (1) (2000) S11–S17. [5] M.L. Turgeon, Clinical Hematology: Theory and Procedures, fourth edition, Lippincott Williams & Wilkins, 2005. [6] J. Roewert Huber, E. Stockfleth, H. Kerl, Pathology and pathobiology of actinic (solar) keratosis – an update, Br. J. Dermatol. 157 (Supplement s2) (2007) S18–S20. [7] Y. Wang, D. Crookes, O.S. Eldin, S. Wang, p Hamilton, J. Diamond, Assisted diagnosis of cervical intraepithelial neoplasia (CIN), IEEE J. Sel. Top. Signal Process. 3 (1) (2009) 112–121. [8] G. Rahmadwati Naghdy, M. Ros, C. Todd, E. Norahmawati, Cervical cancer classification using Gabor filters, in: Proceedings of the First IEEE International Conference on Healthcare Informatics, Imaging and Systems Biology (HISB), July 2011, pp. 48–52. [9] G.H.B. Miranda, J. Barrera, E.G. Soares, J.C. Felipe, Structural analysis of histological images to aid diagnosis of cervical cancer, in: Proceedings of the 25th Conference on Graphics, Patterns and Images (SIBGRAPI), August 2012, pp. 316–323. [10] H. Ahmady Phoulady, B. Chaudhury, D. Goldgof, L. Hall, P. Mouton, A. Hakam, E. Siegel, Experiments with large ensembles for segmentation and classification of cervical cancer biopsy images, in: Proceedings of the International IEEE Conference on System, Man and Cybernetic, San Diego, CA, 2014. [11] A. Tabesh, M. Teverovskiy, H.Y. Pang, V.P. Kumar, D. Verbel, A. Kotsianti, O. Saidi, Multi feature prostate cancer diagnosis and Gleason grading of histological images, IEEE Trans. Med. Imaging 26 (10) (2007). [12] P.W. Huang, C.H. Lee, Automatic classification for pathological prostate images based on fractal analysis, IEEE Trans. Med. Imaging 28 (7) (2009) 1037–1050. [13] S.K. Tai, C.Y. Li, Y.C. Wu, Y.J. Jan, S.C. Lin, Classification of prostatic biopsy, in: Proceedings of the International Conference on Digital Content, Multimedia and Technological Applications, 2010, pp. 354–358.

39

[14] M.M. Krishnan, V. Venkatraghavan, U.R. Acharya, M. Pal, R.R. Paul, L.C. Min, A. K. Ray, J. Chatterjee, C. Chakraborty, Automated oral cancer identification using histopathological images: a hybrid feature extraction paradigm, Micron 43 (2012) 352–364. [15] M.M. Dundar, S. Badve, G. Bilgin, V. Raykar, R. Jain, O. Sertel, M.N. Gurcan, Computerized classification of intraductal breast lesions using histopathological image, IEEE Trans. Biomed. Eng. 58 (7) (2011) 1977–1984. [16] M. Veta, J.P.W. Pluim, P.J. Van Diest, M.A. Viergever, Breast cancer histopathology image analysis: a review, IEEE Trans. Biomed. Eng. 61 (5) (2014) 1400–1411. [17] B. Akbar, V.P. Gopi, V.S. Babu, Colon cancer detection based on structural and statistical pattern recognition, in: Proceedings of the 2nd International Conference on Electronics and Communication Systems (ICECS), 2015, pp. 1735– 1739. [18] S. Rathore, M. Hussain, M.A. Iftikhar, A. Jalil, Ensemble classification of colon biopsy images based on information rich hybrid features, Comput. Biol. Med. J. 47 (2014) 76–92. [19] M.V. Dass, M.A. Rasheed, M.M. Ali, Classification of lung cancer subtypes by data mining technique, in: Proceedings of the International Conference on Control, Instrumentation, Energy and Communication (CIEC), 2014, pp. 558–562. [20] C. Lu, M. Mahmood, N. Jha, M. Mandal, Detection of melanocytes in skin histopathological images using radial line scanning, Pattern Recognit. 46 (2013) 509–518. [21] C. Lu, M. Mahmood, N. Jha, M. Mandal, Automated segmentation of the melanocytes in skin histopathological images, IEEE J. Biomed. Health Inform. 17 (2) (2013) 284–296. [22] C. Lu, M. Mandal, Automated segmentation and analysis of the epidermis area in skin histopathological images, in: Proceedings of the Annual International Conference of Engineering in Medicine and Biology Society (EMBC), 2012. [23] C. Lu, M. Mandal, Efficient epidermis segmentation for whole slide skin histopathological images, in: Proceedings of the 36th Annual International Conference of the Engineering in Medicine and Biology Society, 2014. [24] M. Mokhtari, M. Rezaeian, S. Gharibzadeh, V. Malekian, Computer aided measurement of melanoma depth of invasion in microscopic images, Micron 61 (2014) 40–48. [25] Xu. H. Hongming, M. Mandal, Epidermis segmentation in skin histopathological images based on thickness measurement and k-means algorithm, EURASIP J. Image Video Process. 18 (2015). [26] C. Lu, M. Mandal, Automated analysis and diagnosis of skin melanoma on whole slide histopathological images, J. Pattern Recognit. 48 (2015) 2738–2750. [27] H. Xu, C. Lu, M. Mandal, An efficient technique for nuclei segmentation based on ellipse descriptor analysis and improved seed detection algorithm, IEEE J. Biomed. Health Inform. 18 (5) (2014) 1729–1741. [28] N. Noroozi, A. Zakerolhosseini, Computerized measurement of melanocytic tumor depth in skin histopathological images, Micron 77 (2015) 44–56. [29] K. Korotkov, R. Garcia, Computerized analysis of pigmented skin lesions: a review, Artif. Intell. Med. 56 (2) (2012) 69–90. [30] Q. Abbas, I. Fondón, M. Rashid, Unsupervised skin lesions border detection via two-dimensional image analysis, Comput. Methods Prog. Biomed. 104 (3) (2011) e1–e15. [31] N.M. Shakya, R.W. LeAnder, K.A. Hinton, S.M. Stricklin, R.K. Rader, J. Hagerty, W.V. Stoecker, Discrimination of squamous cell carcinoma in situ from seborrheic keratosis by color analysis techniques requires information from scale, scale-crust and surrounding areas in dermoscopy images, Comput. Biol. Med. 42 (12) (2012) 1165–1169. [32] K. Shimizu, H. Iyatomi, M.E. Celebi, K.A. Norton, M. Tanaka, Four-class classification of skin lesions with task decomposition strategy, IEEE Trans. Biomed. Eng. 62 (1) (2014) 274–283. [33] T.F. Chan, L.A. Vese, Active contours without edges, IEEE Trans. Image Process. 10 (2) (2001) 266–277. [34] E. Alpaydin, Introduction to Machine Learning, second edition, MIT Press, Cambridge, Massachusetts, London, England, 2010. [35] P. Perona, J. Malik, Scale-space and edge detection using anisotropic diffusion, IEEE Trans. Pattern Anal. Mach. Intell. 12 (7) (1990) 629–639. [36] L. Cohen, Time–frequency distributions – a review, Proc. IEEE 77 (7) (1989) 941–981. [37] Jont B. Allen, Short time spectral analysis, synthesis and modification by discrete Fourier transform, IEEE Trans. Acoust. Speech Signal Process. (ASSP) 25 (3) (1977) 235–238. [38] M. Musselman, D. Djurdjanovic, Time–frequency distributions in the classification of epilepsy from EEG signals, J. Expert. Syst. Appl. 39 (13) (2012) 11413–11422. [39] A. Tzallas, M. Tsipouras, D. Fotiadis, Epileptic seizure detection in EEGs using time–frequency analysis, IEEE Trans. Inf. Technol. Biomed. 13 (5) (2009) 703–710. [40] R.O. Duda, P.E. Hart, D.G. Stork, Pattern Classification, second edition, Wiley, 2000. [41] X. Moles Lopez, E. D'Andrea, P. Barbot, A.-S. Bridoux, S. Rorive, I. Salmon, An automated blur detection method for histological whole slide imaging, PLoS One 8 (12) (2013) e82710.