A Whole Slide Image Grading Benchmark and Tissue Classification for Cervical Cancer Precursor Lesions with Inter-Observer Variability

TL;DR

This paper introduces a comprehensive benchmark for grading cervical cancer precursor lesions using whole-slide images, addressing inter-observer variability and proposing a novel tissue classification method based on morphological features.

Contribution

It presents a new whole-slide image grading benchmark, a tissue classification method utilizing a novel morphological feature, and an analysis of inter-observer variability among pathologists.

Findings

The tissue classification method achieved promising accuracy.

Inter-observer variability significantly impacts diagnosis consistency.

The benchmark facilitates standardized evaluation of grading algorithms.

Abstract

The cervical cancer developing from the precancerous lesions caused by the Human Papilloma Virus (HPV) has been one of the preventable cancers with the help of periodic screening. There are two types of grading conventions widely accepted among pathologists. On the other hand, inter-observer variability is an important issue for final diagnosis. In this paper, a whole-slide image grading benchmark for cervical cancer precursor lesions is introduced. The papillae of the cervical epithelium and overlapping cell problems are handled and a tissue classification method with a novel morphological feature exploiting the relative orientation between the BM and the major axis of all nuclei is developed and its performance is evaluated. Besides, the inter-observer variability is also revealed by a thorough comparison among pathologists' decisions, as well as, the final diagnosis.

Tables

Table 1. Table 1: Number of epithelium pieces in each class depending on the CIN-based grading

		Normal	CIN1	CIN2	CIN3	Total
CIN-Based Grading	SEP	471	240	107	139	957
	LEP	150	79	34	57	320

Table 2. Table 2: Number of epithelium pieces in each class depending on the SIL-based grading

		Normal	LSIL	HSIL	Total
SIL-based Grading	SEP	471	240	246	957
	LEP	150	79	91	320

Table 3. Table 3: List of morphological features extracted in the proposed tissue classification method

Features		Description
Average nucleus area (ANA)		the average nucleus area of each region
Average cytoplasm area (ACA)		the region which represents the subtraction of total nucleus area from total area of each region
Nucleus-cytoplasm ratio (NCR)		the division result of total nucleus area to the total the cytoplasm area in each region
Nucleus perimeter (NP)		the length of the line which surrounds the nucleus in each region
Border irregularity (BI)		obtained by dividing length of the uniform ellipses that fit the nucleus to the circumference of the related nucleus
Hyperkromosis index (HI)		represents the standard deviation value of the parabasal cells with respect to the cells of the same lesion
Polarity lossindex (PI)		The average angle between the basal membrane and the major axis of all nucleus

Table 4. Table 4: Agreement/disagreement of the pathologists experts in diagnosis of SEP with respect to CIN-based grading.

	Pathologist 2
Pathologist 1		Normal	CIN1	CIN2	CIN3	Total
	Normal	354	37	0	0	391
	CIN1	88	156	8	0	252
	CIN2	13	52	71	15	151
	CIN3	4	9	22	128	163
	Total	459	254	101	143	957

Table 5. Table 5: Agreement/disagreement of the pathologist 1, pathologist 2 and the proposed method with final diagnosis in diagnosing of SEP with respect to CIN-based grading.

CIN Grading		P1				P2				Proposed
		N	C1	C2	C3	N	C1	C2	C3	N	C1	C2	C3
Final Diagnosis	N	377	74	18	2	424	44	3	0	386	61	17	7
	C1	13	176	37	14	26	193	16	5	116	91	19	14
	C2	0	2	87	18	8	17	79	3	13	30	47	17
	C3	1	0	9	129	1	0	3	135	4	15	18	102
	Tot.	391	252	151	163	459	254	101	143	519	197	101	140

Table 6. Table 6: Agreement/disagreement of the Pathologists with final diagnosis in the diagnosis of SEP with respect to SIL-based grading.

SIL Grading		P1			P2
		N	LSIL	HSIL	N	LSIL	HSIL
Final Diagnosis	N	377	74	20	424	44	3
	LSIL	13	176	51	26	193	21
	HSIL	1	2	243	9	17	220
	Total	391	252	314	459	254	244

Table 7. Table 7: Agreement/disagreement of the proposed method and Delaunay Triangulation(DT) with final diagnosis in the diagnosis of SEP with respect to CIN-based grading.

CIN Grading		Proposed				DT
		N	C1	C2	C3	N	C1	C2	C3
Final Diagnosis	N	386	61	17	7	471	0	0	0
	C1	116	91	19	14	240	0	0	0
	C2	13	30	47	17	107	0	0	0
	C3	4	15	18	102	139	0	0	0
	Total	519	197	101	140	0	0	0	0

Table 8. Table 8: Agreement/disagreement of the proposed and DT method with final diagnosis in diagnosis of SEP with respect to SIL-based grading.

SIL Grading		Proposed
		N	LSIL	HSIL
Final Diagnosis	N	381	57	33
	LSIL	106	91	43
	HSIL	16	27	203
	Total	503	175	279

Equations

$$N=\sqrt{\frac{w\times h}{k}}$$

$$d_{rgb}=\sqrt{(r_j-r_i{)}^2+(g_j-g_i{)}^2+(b_j-b_i{)}^2}$$

$$d_{xy}=\sqrt{(x_j-x_i{)}^2+(y_j-y_i{)}^2}$$

$$d_s=d_{rgb}+m/(N\times d_{xy})$$

$$p(\mathbf{x})=\sum _{i=1}^K{\varphi }_i\mathcal{N}(\mathbf{x}|{\mu }_i,{\Sigma }_i)$$

$$\sum _{i=1}^K{\varphi }_i=1$$

Full text

Introduction

Cervical cancer is one of the most commonly seen cancer type in the world and the 4th most common cause of death, which develops from precursor lesions [1]. Studies indicate that, almost all cervical cancer cases develop by the effect of Human Papilloma Virus (HPV), which reaches epithelial basal layer cells with the help of micro-injuries in the cervical epithelium. Carcinogenic effect of the virus occurs when HPV’s genome integrates with the cell genome [2, 3, 4]. This effect, which requires a certain period of time, appears as morphological changes in the cervical epithelium. These precancerous lesions characterized by dysplastic changes are called squamous intraepithelial lesions (SIL).

Impact of HPV on the cervical epithelium varies throughout the life cycle of the virus, which in turn results in different morphological changes. Early diagnosis can be made possible by the analysis of these morphological structures [1, 5]. After being infected by HPV, basal cells proliferate and the epithelium loses its maturation. As well as the loss of maturation which results in polarity loss in the epithelium, cells show nuclear enlargement, nuclear irregularity, and hyperchromasia. Depending on the proliferation process, the number of mitoses also increases. The effect of viral proteins on the cyto-skeleton reveals halo cells with characteristic perinuclear halo named “koilocytes” (it means "hollow" in Greek). These dysplastic changes are graded according to whether they are seen among the lower, middle and upper part of the epithelium. They are reported according to Cervical Intraepithelial Neoplasia (CIN) 1-3 in the CIN-based grading and LSIL (low-grade squamous intraepithelial lesion) and HSIL (high-grade squamous intraepithelial lesion) in the SIL-based grading [6, 7, 8]. Currently, the use of SIL-based grading is recommended, yet the CIN-based grading is also used.

Pathologic diagnosis of cervical biopsies varies depending on biopsy ingestion or artifacts due to laboratory steps, and pathological interpretation. Due to the spread of women health screening programs, the diagnosis of cervical biopsies is frequently encountered and this diagnosis variability has become a more important problem. Cervical biopsy interpretation has inter- or intra-observer variability, which means that a biopsy may have different diagnoses by different pathologists or by the same pathologist at different times, and it is accepted to an extent in the literature [9, 10, 11]. Studies have been made to overcome this problem with classification systems suitable for the nature of HPV or with the help of immunohistochemical techniques [12, 2, 7].

Increasing role of the information technology (IT) on the area of medicine has a positive impact on the pathology. Digital pathology that includes diagnosis, education, consultation, archiving, and also morphometric evaluation tools [13, 14]. Studies about morphometric analysis are available for different tissues and systems [15], as well as for cervical lesions [16, 17, 18, 19] in the literature. De et al. [16] studied image analysis methods on 62 digital images of cervical epithelial with Normal, CIN1, CIN2 and CIN3 labeled lesions. The cervical regions are manually marked by the pathologist on selected epithelial images, and these regions are divided into vertical segments by calculating the medial axis. The obtained epithelial segments are examined in terms of structural, geometric, and profile-based properties. Contrast, energy levels, pixel correlation values, and neighborhood features of the pixels within the vertical segment are studied as the structural features. Geometric features include the distances between nuclei centers and Delaunay triangulation. In the profile-based feature extraction, correlation values and the brightness values of all pixels of each row of vertical segment is calculated. Linear Discriminant Analysis (LDA) and Support Vector Machines (SVM) are utilized to classify feature vectors of vertical segments. First, each of the vertical segments is classified individually, then these decisions are fused to obtain a whole epithelium classification result. The effect of individual decisions of vertical segments on the whole epithelium classification result is also examined. One-to-one correspondence between the system result and pathological diagnosis is named as "Exact Class Label" (1st approach), only one class difference between system result and pathological diagnosis is named as "Windowed Class Label" (2nd approach) and bigger differences between system result and pathological diagnosis is named as Normal versus CIN (3rd approach). Different classification performances are calculated using different approaches and features. Using all the structural, geometric and profile-based features, recognition rate of 62.3% on vertical segments and 39.3% on whole epithelium is reached.

Guo et al. [17] tried to develop enhanced image analysis methods on the cervix image data set that was formed in De et al. [16]. They have increased the success of classification by adding structural features of the nucleus and cytoplasm in addition to the features extracted from similar vertical segments as in De et al. [16]. These features consist of nucleus, cytoplasm, and acellular areas and ratios, color scale (red, green, blue) brightness values, numbers of triangles obtained by Delaunay triangulation at upper, middle, and lower epithelium regions. The features are classified by the classification methods of the previous study. The name "Windowed Class Label" used in the second approach in the previous study is changed to "Of-By-One Class Label". In this study, as well as using the same data set as De et al., they made a difference of examination by 2 different pathologists. The diagnostic success rates of the extracted features are determined by the Attribute Information Gain Ratio (AIGR) algorithm. They evaluated the success of features according to two different classification approaches. As a result of adding structural features of cervical regions, they have increased their classification success up to 82-88.5%.

Wang et al. studied morphometric analysis methods on 31 digital images of cervical biopsies [19]. Their study consists of two steps as the automated segmentation of squamous epithelium and the CIN classification. In the first step, the epithelium is segmented using the difference of the visual properties of five different regions consisting of squamous epithelium, columnar epithelium, stroma, background and erythrocytes. The medial axis is drawn parallel to the basal and upper membrane borders after the epithelial region is segmented. Square windows with $250\times 250$ pixel dimensions are created on normal lines of medial axis. The features of the nucleus average area and number, the average area of the triangles obtained by Delaunay triangulation, and the average edge length are analyzed within each window. Obtained feature vectors are fed to different classification methods, and have reached accuracy rates ranging from 60% to 95%.

Keenan et al. proposed a study to analyze 230 digital cervix images consisting of normal, koilocytosis, CIN1, CIN2, and CIN3 labeled lesions [20]. The features of the nucleus area, the nucleus cytoplasm ratio, the ratio of nucleus area to cytoplasm, and the edges / areas of the Delaunay triangles are analyzed. The Kappa value for the observer difference between the two pathologists involved in the study is 0.415. The classification performance of the system in distinguishing normal and CIN lesions is 98.7%. The overall success rate is 62.3% where the worst performance is achieved on CIN2 labeled patterns.

Nagdhy et al. proposed a study to classify a total of 475 cervical biopsies with normal, CIN1, CIN2, CIN3, and invasive carcinoma using three different methods [18]. The nucleus area, core cytoplasm ratio, core boundary irregularity, and areas of Delaunay triangles are analyzed. They reached up to 97% with respect to specificity and 100% with respect to sensitivity using different methods including Gabor-based texture descriptor, GLCM (Gray-level co-occurrence matrix) texture descriptor, and pre-trained convolutional neural network.

In this study, morphometric analysis methods for the cervical SIL diagnosis is investigated on a new digital cervical image data set. The numerical values of the morphological features used by the pathologists in the diagnosis are extracted. The statistical significance of their contribution to diagnosis is examined. Within the scope of the study, a Computer Aided Diagnostic Auxiliary System (CADAS) is developed and its performance is evaluated.

Contributions of this study are as follows:

• •

A new whole slide image grading benchmark for grading of cervical dysplasias is created and introduced to histopathological image analysis community.

• •

Images obtained from the data set are labeled by two pathologists to mention the inter-observer variability in cervical dysplasia grading.

• •

Pathologists diagnosed each image patch stained with hematoxylin $\epsilon$ eosin (H$\epsilon$E) in the data set independently. In the likely case of inconsistent diagnoses, the image patches that are stained with p16 and Ki67 immunhistochemical dyes are analyzed to decide a final diagnosis.

• •

A morphometric analysis method for cervical SIL diagnosis is proposed.

• •

The presence of papillaries in the dataset that leads to tangential sections is one of the important parameters that pathologists give account for when diagnosing.

1 Materials and Methods

This study is conducted by a group of scientists and medical researchers. Cervical tissue slide samples with diagnosis results are collected in the pathology laboratory of Istanbul Medipol University (IMU) Hospital, Istanbul, Turkey. Fig. 1 shows the processing steps followed in this study for the proposed CADAS to grade cervical cancer precursor lesions.

1.1 Data Collection and Image Acquisition

Within the scope of the study, 127 high resolution slides from 54 patients are scanned at the pathology laboratory of IMU Hospital. Fig. 2 represents the whole slide images obtained from the data set. The images stained with H$\epsilon$ E, Ki67, and p16 immunhistochemical dyes. All high-resolution images are then divided into 957 small epithelium pieces by the pathologist. Each slide in the data set is diagnosed after splitting into smaller epithelial pieces. Totally, 957 epithelial pieces are obtained from the whole slides. 471 of the 957 images are diagnosed as normal, 240 of them are diagnosed as CIN1, 107 of them are diagnosed as CIN2 LSIL, 57 of them are diagnosed as CIN3.

The images of the hematoxylin and eosin (H$\epsilon$E), p16 and Ki67 preparations are acquired by an off-the-shelf whole slide scanner (See Fig.3). The scanner has a capability of up to $20\times$ optical and $40\times$ digital zoom. The whole slide images obtained with the high-resolution scanner are transferred to the digital platform to be processed by several image processing techniques and also to be interpreted by the expert pathologists. The images are saved in TIFF format without any loss. The size of the images obtained by the scanner are varied from $10,000\times 10,000$ to $80,000\times 80,000$ and there are more than one diagnosis in a single lesion.

1.2 Ethics Statements

Authors confirm that all samples taken from patients were prepared in accordance with the legislation prepared by the Ministry of Health of Turkey and in accordance with international agreements & European Union standards. All experimental protocols were approved by the Istanbul Medipol University’s licensing committee. Informed consent was obtained from all subjects whose tissue samples were used in experiments. In tissue sample collection for data set, there were no subjects under 18.

1.3 Annotation and Image Labeling

A graphical user program is developed within the scope of this study for pathologists to mark/label the basal membranes (BM) and papillae of the cervical epithelium. After marking the membranes and the papillae, the program extract the hot spot region from the background. Fig. 4 represents an original input image taken from the data set and a clean image after marking the coordinates of the epithelium. Further image processing and analysis algorithms use this clean image.

At first, two pathologists made the diagnosis independently for each small epithelial piece (SEP) image patch. A final diagnosis is then made by observing the same lesions stained with p16 and Ki67 immunohistochemical dyes in case of disagreement. According to the final diagnosis, 471 of SEP (%49.2) are labeled as normal, 240 of them (%25.1) are CIN1, 107 of them (%11.2) are CIN2 and 139 of them (%14.5) are CIN3. However, 150 of large epithelial piece (LEP) (%46.9) are labeled as normal, 79 of LEP (%24.7) are CIN1, 34 of LEP (%10.6) are CIN2 and 57 of LEP (%17.8) are CIN3 (see Table 1). Diagnostic distributions of the SIL-based grading are shown in Table 2. Similarly, distribution of final diagnosis in SIL-based grading are as follows: 471 of SEP (%49.2) are normal, 240 of SEP (%25.1) are LSIL and 246 of SEP (%25.7) are HSIL. Similarly, 150 of LEP (%46.9) are normal, 79 of LEP (%24.7) are LSIL, 91 of LEP (%28.4) are HSIL.

1.3.1 Inter-Observer Variability

Interpretations of morphologic changes representing dysplasia may differ between pyhsicians or for the same physician in different time intervals. This variety can be interpreted as Inter/Intra-observer agreement / disagreement. Artifacts associated with the biopsy procedure and tangential sections in the microscopic examination are also effective on this variety. Inter- and intra-observer agreement rates are in the range of 0.20 and 0.47 in the literature [10, 11]. Regarding to CIN-based grading, SIL-based grading provides higher inter-observer and intra-observer agreement rates. The highest diagnosis diversity is reported between the groups of CIN2, while the lowest is CIN3. The disagreement rates are smaller between normal and CIN1 groups. McCluggage et  al.  reported weak inter-observer agreement in the CIN-based grading with Kappa value of 0.2. Although the compatibility rates are reported low, Kappa value is found as 0.3 in the SIL-based grading. Failure to achieve the expected high agreement rates is interpreted by the pathologists involved in the study not being familiar with SIL-based grading. The same test is repeated with the observers who have been experienced to use the SIL-based grading for six months more and new Kappa values are calculated as 0.33 (intra-observer) and 0.47 (inter-observer). Galgano et al. tried to maximize the agreement rates between the observers with P16 and Ki67 immunohistochemical methods (54). The Kappa value is found to be 0.68 by immunohistochemistry examination while standard hematoxylin and eosin (H$\epsilon$E) detection has the kappa value of 0.47. In the study, it is stated that the low agreement rates associated with diagnostic differences can be increased by using SIL-based grading rather CIN-based grading, or utilizing some immunohistochemical methods aiding diagnosis.

1.3.2 Final Diagnosis

Immunohistochemical examinations are used as an assistive method to obtain the diagnosis in case the morphological features are not clearly interpreted. P16, Ki67 and ProExC are the most widely used immunohistochemical studies for cervical precursor lesions [21, 22, 23, 24, 25]. Staining pattern with p16 is important in immunohistochemical evaluation, and block-like and strong staining demonstrates HrHPV association, with at least 1/3 of the epithelium (7,56,58). Ki67 is an indicator of proliferation. Positivity may also be seen in other proliferating cells such as inflammatory cells as in keratinocytes. For this reason, it must be interpreted carefully in the presence of inflammation. ProExC is similar to Ki67 in terms of being a proliferation indication and its staining type. P16 and Ki67 are frequently used in routine practice. HrHPV-associated lesions show strong “nuclear” or “nuclear and cytoplasmic”, block-like staining with P16. The squamous metaplasia, atrophy, reactive regenerative changes that appear in the SIL discriminator pattern show a negative staining pattern. While Ki67 normally stains parabasal cells, positivity is also observed in higher epithelial sections in relation to the grade of dysplasia in SIL. Reactive cells and inflammatory cells may also exhibit immunoreactivity, it should be very careful when interpreting these tissues.

1.4 Morphometric Feature Extraction and Tissue Classification

In this study, a morphological analysis based feature extraction method is used for the grading of cervical cancer precursor lesions. The processing steps followed in the study is represented in Fig. 1. The first row of the diagram describes the pathological pre-process which is handled in the pathology laboratory. Whole slide scanning and the filing process mentioned in the second row is done by the medical researcher and the computer scientist in collaboration. This section describes the CADAS framework which is mentioned in the third row.

1.4.1 Creating the Small Epithelial Pieces (SEP)

In Fig. 5, the red line corresponds to the coordinate information of the BM, while the green line corresponds to the coordinate information in the upper membrane. Determining the basal and upper membrane coordinates allows to know in which region of the epithelium the cells are located. After the region of interest (squamos epthitelium) has been obtained, an interface developed within the scope of the study is used to divide the whole epithelium to SEP which can be assumed equal in length (see Fig. 6).

The image patches which are analyzed in this study are represented in Fig. 6. The coordinates of basal and upper membranes of the epithelium are marked by the pathology experts with the use of a graphical interface. Coordinates data information of the papillae which represented with yellow line are also stored in separated files.

1.4.2 Obtaining Cells by Simple Linear Iterative Clustering(SLIC) Superpixels Segmentation Algorithm

After small epithelial pieces are obtained, the high resolution histopathological images are ready for further analysis. First, a median filter of dimension $9\times 9$ is applied to the image to remove the artifacts without effecting the boundaries. Then, cellular structures have been obtained by simple linear iterative clustering (SLIC) superpixels segmentation algorithm, which is one of the methods that have not been widely used yet in histopathological images in recent years. This method performs the segmentation process based on the color similarities and neighbour relations of the pixels in the image [26]. The grid size $N$ is expressed as

[TABLE]

where $k$ is the number of superpixels for a given input image, $w$ and $h$ represent the width and height of the given image patch, respectively. The euclidean distance of the related pixel to the superpixel center is

[TABLE]

where, $i$ represents the value to be clustered and $j$ represents the center pixel. Here, $r,g$ and $b$ represent the brightness values of red, green, and blue color of the respective pixels. RGB color space is used in this study instead of using Lab color space as mentioned in Achanta et al[26]. The Eq. 3 also represents the distance of the coordinates of each pixel to the related cluster center:

[TABLE]

where, $x_{j}$ and $y_{j}$ are the horizontal and vertical coordinate information of each center pixel, and $x_{i}$ and $y_{i}$ values are the coordinate information of each pixel to be clustered:

[TABLE]

the value of $d_{s}$ is the sum of the (x, y) plane distance normalized by the grid interval N and the RGB distance. Here, normalization is done so that the calculation of the coordinate information does not directly affect the brightness interval. The value of $m$ is defined to set the compactness of superpixels.

According to the SLIC method applied in this study, the cellular structures become more compact and can be separated from the background when each obtained SEP image is expressed with 3000 superpixels. A cellular structure in cervical precursor lesion is approximately $20\times 20$. A crucial point to note here is that the superpixels’ sizes should not exceed the size of cellular structure. As can be calculated from the Eq.1, choosing at least 2000 superpixels will guarantee most of the superpixels to not exceed the size of a cellular structure. Less number of superpixels cause overlapping cells. It can be quite difficult to distinguish cellular structures, especially those close to the BM. Since the superpixels which represent the cellular structures are darker than the superpixels which represent the background, the elimination of the unwanted pixel groups (fat-like tissue) above a certain threshold level clears the background. At this stage, small artifacts similar to the cellular structures and some inflammation can remain with the cells as a foreground information. These structures can be eliminated with a morphological size operation that can be applied to the binary image after segmentation stage. Final segmentation result of a sample SEP image patch obtained from the data set is shown in Fig. 8.

1.4.3 Handling the Overlapping Cells Problem

Following the morphological operations, overlapping cells are separated. In the literature, the problem of overlapping nuclei in a plurality of nucleus segmentation studies has been encountered. Because the presence of overlapping cell structures significantly reduces the success of CADAS. Solving the problem of overlapping nucleus at this point is very crucial as a significant contribution in this area. In our study, it is observed that after the segmentation process, there are a large number of overlapping nucleus structures, especially around the BM.

Overlapped cells are intensively present on the SEP image patches. The problem of overlapping of cells should be handled in order to obtain the morphological characteristics of the cell nuclei. Binary images segmented by using SLIC algorithm usually consist of small cellular-like noisy parts. These unwanted small pixel groups are eliminated by an automatic method, which clears pixel groups smaller than 50 pixel. Therefore, the circumference of the cells is also quite rough after the segmentation process. As shown in Fig. 9(a), the closing process has been applied to make the binary image more compact. In order to obtain cell centers, the distance transform is applied to get local maxima shown in Fig 9(c). Local maxima are estimated by using Algorithm-1.

The ellipse form is able to model the cell shapes mathematically well. Thereby, Gaussian Mixture Model (GMM) is one of the best candidate for ellipse fitting over cell heaps. GMM is one of the most common algorithms for statistical data modeling [27]. The main purpose of the algorithm is to express the distributions as the sum of the weighted Gaussian mixtures (see Eq. (5)). $\mathcal{N}(\mathbf{x}|\mathbf{\mu}_{i},\Sigma_{i})$ and $\phi$ intend normal distribution which has mean $\mu_{i}$ and covariance matrix $\Sigma_{i}$ and its weight parameter.

[TABLE]

Fig. 10 shows an example of overlapping nucleus and how these overlaps are resolved in a small patch of the image obtained from the data set. The basic structure of the algorithm that determines the cell overlapping is based on the determination of local maxima from the distance of the cell centers to the boundaries. Once the cell centers are determined, the distance transformation yields the value of the multiplexing for each pixel. The processing steps applied for center estimation, distribution generation and multiplexing are given in Algorithm-1. $\mathbf{C}$ and $\mathbb{H}$ refer to row-wise cell centers matrix, and multiplexed coordinate set, respectively.The $\alpha$ parameter obtained from the distance transform indicates how many times the corresponding coordinate will be repeated in the set $\mathbb{H}$. Thereby, the distance value of the pixels away from the border are higher, so the amount of these in $\mathbb{H}$ will be more. Pixel distribution becomes more suitable for GMM. For example, if each pixel is far away from a boundary, the coordinate information of the related pixel is multiplexed. By applying GMM on the obtained multiplexed coordinate distribution, suitable ellipses are obtained for each cell (see Fig. 11).

1.4.4 Obtaining the Morphological Features of Each SEP

After the cell segmentation and elimination of the cell overlap problem, several morphological features of each cell are extracted. Table 3 represents the morphological features extracted for grading each SEP image patch in this study. Average nucelus area (ANA), average cytoplasm area (ACA), nucleus - cytoplasm area ratio (NCR), nucleus perimeter (NP), border irregularity (BI), hyperkromasis index (HI), and polarity loss index (PLI) are the features represented from the first row to the end, respectively. ANA defines the average nucleus area while ACA is the average cytoplasm area. NCR describes the division results of nucleus ratio to the cytoplasm ratio. NP is the average length of the border pixels of nucleus. BI is the divison of the surrounding length of each pixel to the ellipse that fit each nuclei. The value which represent the hyperkromasis of cell is calculated by taking standard deviation of the pixel intensity values of the cellular structure. PLI is determined by calculating the magnitude of each cellular structure to the BM. All the features are extracted for each region shown in Fig. 12

Since the emphasized morphological features change depending on their distance to the basal and the upper membranes, each image segment is divided into three main regions as represented in Fig. 12. Then, morphological features related to each region are extracted and stored for further analysis in grading the SEP image patch.

The data set of feature vectors is imbalanced. Different classifiers have been proposed in the literature for imbalanced data sets [28, 29, 30]. The Weighted k-Nearest Neigbour (w-kNN) algorithm is one of these. In this study, w-kNN algorithm is preferred because of its fast operation and practical use. It is also another important criterion for selecting that it gives successful results for imbalanced data sets [31, 32]. The w-knn algorithm looks at the k closest neighbors class as the k-NN algorithm. In addition, for each neighbor, the weight w defined in $\textit{w}=\frac{1}{d(\textbf{x}_{i},\textbf{x}_{q})^{2}}$ is assigned to classify according to the weight of the classes. $d(.,.)$ is Euclidean distance function. If neighbour sample $\textbf{x}_{i}$ is far from query sample $\textbf{x}_{q}$, the effect on the classification is weak, and vice versa.

2 Results

The similarities and differences in the diagnoses of SEP image sections in the data set given by the two pathologists with respect to CIN-based grading. Diagonal values refer to the number of images which have the same diagnosis of two pathologist. The agreement ratio of the pathologists is 73% in the classification of the SEP images with respect to CIN-based grading are shown in Table 4.

The agreement/disagreement of pathologist 1 and pathologist 2 to the final diagnosis of each SEP image patch are presented in Table 5. Final diagnosis is determined according to the disagreement of pathologists for a SEP image patch. The SEP that are not labeled as the same by the pathologists are then observed from the same tissue stained with p16 and Ki67 immunhistochemical dyes.

Table 5 represents the agreement between the final diagnosis and two pathologists with respect to CIN-based grading system. It can be observed that pathologist 2 has more compatible diagnosis result than pathologist 1 with final diagnosis. However, pathologist 1 has more consistent diagnosis in CIN2 SEP image patches. An important information to be drawn from the table is that the number of windowed classes (labelling CIN1 instead of Normal tissue; Normal or CIN2 instead of CIN1; CIN1 or CIN3 instead of CIN2 and CIN2 instead of CIN3) is high over.

Another system that pathologists pay attention to while diagnosing tissues is the SIL-based grading system. In this system, the CIN2 grade is assumed to be two level, CIN3-like and CIN1-like. CIN2 lesion which resembles CIN3 and CIN3 are expressed as HSIL; CIN2 which resembles CIN1 and CIN1 are expressed as LSIL. The treatment of precursor lesions of cervical cancer varies according to LSIL and HSIL. Table 6 represents the agreement between the final diagnosis and two pathologists with respect to SIL-based grading system. Pathologist 1 has more accurate results than pathologist 2 in normal and LSIL while pathologist 2 has more accurate result in diagnosing HSIL. If the diagnosis agreement of the pathologists according to the Table 5 and Table 6 are compared, it can be observed that pathologists makes more consistent diagnosis in SIL-based grading system.

Table 7 represents the agreement between the final diagnosis with the proposed method and the DT method which is the one of the best known algorithms used in diagnosis of cervical cancer grading depending on CIN-based grading system.

Table 7 represents the agreement between the final diagnosis and the proposed method with respect to CIN-based grading system. Normal and CIN3 SEP patches prediction are classified accurately. However, predicting CIN1 and CIN2 SEP patches is lower while comparing with CIN1 and CIN3. The classification accuracy of the proposed method is approximately $65.4\%$. If the results obtained from the proposed method is compared with the pathologists, it can be said that the CAD system developed in this study should be improved in order to be used as a secondary decision system to help the pathologist in predicting cervical cancer precursor lesion grade.

Table 8 represents the agreement between the final diagnosis and the proposed method with respect to SIL-based grading system. The classification performance of the proposed study is again less accurate than the pathologist. However, the results obtained from the SIL-based grading system of the proposed method is improved to $70.5\%$ comparing to the CIN-based grading system.

3 Discussion

In this study, it is aimed to translate the evaluations of pathologists which have subjectivity on cervical dysplasia to into numerical values and to develop a ‘Computer Assisted Diagnostic Auxiliary Systems (CADAS)’ as a result. Our study on cervical dysplasia has the largest data set according to the similar studies available in the literature. Furthermore, the fact that the diagnoses are given by two pathologists, and the reassessment and determination of the definitive diagnosis during inconsistent cases increased the reliability of the CADAS training set. The developed CADAS promises to be used as an assistant system in the future because of numerical values that are found to be in parallel with the diagnostic parameters used by the pathologists (such as ratio of nucleus to cytoplasm, nucleus boundary irregularity, polarity loss and hyperchromaticity) and statistically significant. The studies in the literature are mostly designed by engineers and the contribution of pathologists is very limited. For this reason, there are some shortcomings when viewed from the perspective of pathology and clinical approach. In the development of a CADAS to be used in pathology, the presence of pathologists at every step is a necessary requirement.

4 Conclusion

In this paper, we present a new benchmark data set of cervical cancer precursor lesions, which we make available to the scientific community for grading the cervical intraepithelial neoplasia. Each image in the data set is labeled by two pathologist experts to reveal the inter-observer variability. In case of different diagnoses, p16 and Ki67 immunohistochemical dyes are used to decide a final diagnosis (ground truth). There are also papilla areas that seriously affect the performance of automated methods which makes this study unique as far as we know. A morphological analysis based feature extraction method is also proposed in the study for the grading of cervical cancer precursor lesions. The result of the study is also compared with each pathologist expert and the ground truth. The results show that CAD systems could be used as a secondary decision system for experts with some improvement. It is aimed to improve the classification performance of our CAD system by developing up-to-date image processing and machine learning algorithms especially types of deep learning.

5 Data Availability

The materials and datasets generated during and/or analysed during the current study are available from the http://simplab.yildiz. edu.tr/?q=sources on reasonable request.

Conflict of Interest

All authors declare that they have no conflict of interest.

Competing Interests

The authors declare no competing interests.

Acknowledgement

Authors also would like to thank Argenit Company and Istanbul Medipol University Hospital for providing and annotating the whole slide histopathological images of cervical cancer precursor lesions image data set. The authors state no conflict of interest and have nothing to disclose.

Author contributions statement

A.A. , A.U. and N.C. designed and performed the experiments, G.B., L.D.A and B.U.T were involved in planning and supervised the work, A.C., I.T. and B.M. contributed to the design and implementation of the research. All authors reviewed the manuscript.