2.1.

Extraction of Cell Nuclear Area

We used cervical smear samples collected at the Department of Gynecologic Oncology, Saitama Medical University International Medical Center. These samples were applied to slides, fixed with 95% alcohol, and subjected to Papanicolaou staining. Squamous cells in these samples were observed and imaged at $1000\times$ ($10\times$ eyepiece and $100\times$ objective lenses) magnification with an optical microscope (AXIO imager A1; Carl Zeiss Ltd., Oberkochen, Germany) attached to a cooled charge-coupled device camera (256 shades of gray) and three transmission filters of red, green, and blue. In each shooting, exposure time and white balance were fixed. In this study, we targeted squamous cells only.

Cells were estimated from these images by a pathologist and two cytotechnologists according to the Bethesda system.²¹ The classifications are as follows: negative for intraepithelial lesion or malignancy (NILM), atypical squamous cells of undetermined significance (ASC-US), low-grade squamous intraepithelial lesion (LSIL), high-grade squamous intraepithelial lesion (HSIL), atypical squamous cells but cannot exclude HSIL (ASC-H), and squamous cell carcinoma (SCC). ASC-US and ASC-H, respectively, represent intermediate classifications between NILM and LSIL and between LSIL and HSIL. In this paper, we avoided these intermediate classifications and used only typical cells (NILM, LSIL, HSIL, and SCC). In addition, we divided NILM cases into three types: normal cell (NOR), metaplastic cell (MET), and regenerative cell (REG). MET and REG include reactive nuclear atypia, which is occasionally difficult to discriminate from neoplastic nuclear atypia [LSIL, HSIL, carcinoma in situ (CIS), and SCC]. Therefore, it is important to classify these cell types using image-based characterization system. Although several detail studies included the image-based cell classification system, few were investigated using the detailed NILM classification.^16,20,27 We also divide SCC cases as CIS and SCC because cytologists usually distinguish these two categories. We, therefore, evaluated seven types of cells: NOR, MET, REG, LSIL, HSIL, CIS, and SCC. Figure 1(a) shows representative images including CIS cells.

Fig. 1

Examples of cervical cytological images. (a) Original image and (b)–(f) mask extraction image of each cell nuclei: NOR, MET, REG, LSIL, HSIL, CIS, and SCC, respectively.

Subsequently, we manually extracted cell nuclear regions from images to generate masking images and transformed these from RGB color to gray-scale according to the $Y$ value in the YCbCr color system. Figures 1(b)–1(h), respectively, show examples of the gray-scale masking images of NOR, MET, REG, LSIL, HSIL, CIS, and SCC. Although automatic methods for extracting cell and cell nuclear regions have been proposed,⁸ these are not completely accurate. We, therefore, manually extracted the studied cell nuclear regions.

2.2.

Conventional Feature Values

Previous studies of the feature quantification of cervical cytology images have used feature values related to the cell nuclear size and shape.^8,13,20,28 We also use eight feature values (nuclear area $f_1$, nuclear perimeter $f_2$, nuclear longest diameter $f_3$, nuclear shortest diameter $f_4$, nuclear convex hull area $f_5$, nuclear convex hull perimeter $f_6$, nuclear circularity $f_7$, and nuclear extension $f_8$). In this paper, convex hull refers to a convex line (i.e., shape of a rubber band) surrounding the nuclear outline. $f_7$ and $f_8$ are represented by the following equations:

Murata et al.¹⁶ used other nuclear shape values to evaluate images of thyroid tumor cytology specimens. These values can be expressed using the following equations:

Feature values have also been proposed for chromatin distribution in the cell nucleus, including the average value, the number of maximum value, and the number of minimum value of the image pixel intensities in nuclear regions.^8,13,28 In addition, the RD value²⁰ represents the difference in average intensity between the center and periphery of the cell nucleus and has been suggested. Regarding the gray-scale masking images described in Sec. 2.1, we also use these four chromatin distribution feature values, represented as $f_{11}$, $f_{12}$, $f_{13}$, and $f_{14}$, respectively. Murata et al. also used skewness, kurtosis, the coefficient of variation, and the upper 20 percentile ratio of the intensity histogram. The coefficient of variation is defined as the ratio between the average value and standard deviation. We also use these values and denote as $f_{15}$, $f_{16}$, $f_{17}$, and $f_{18}$, respectively.

Murata et al.^10,16–18 used 15 texture feature values, in which 10 are Haralick feature values¹⁴ calculated using cooccurrence matrices and 5 are run-length feature values¹⁵ calculated using run-length matrices. Both matrices are calculated from gray-scale intensities within cell nuclei. The 10 Haralick feature values are contrast ($f_{19}$, contrast of intensity), energy ($f_{20}$, uniformity of intensity and texture), correlation ($f_{21}$, correlation of intensity and texture), variance ($f_{22}$, variance of intensity), entropy ($f_{23}$, diversity of intensity and texture), sum variance ($f_{24}$, contrast of intensity and texture), sum entropy ($f_{25}$, diversity of intensity), difference variance ($f_{26}$, variance of texture), difference entropy ($f_{27}$, diversity of texture), and inverse difference moment ($f_{28}$, homogeneity). The five run-length feature values are gray level nonuniformity ($f_{29}$, ununiformity of intensity), run percentage ($f_{30}$, ununiformity of intensity and texture), short run emphasis ($f_{31}$, magnitude of high frequency), long run emphasis ($f_{32}$, magnitude of low frequency), and run-length nonuniformity ($f_{33}$, nonuniformity of texture).

The cooccurrence matrix represents the appearance frequency $\mathbf{P}$ $[=P(i,j)]$ of pixel intensities on a gray-scale image, where $i$ is the intensity of a pixel of interest $\mathbf{A}$ and $j$ is the intensity of a pixel $\mathbf{B}$ near $\mathbf{A}$. We used gray-scale images of 256 gradations, and the matrix size became $256\times 256$. Multiple cooccurrence matrices can be generated using the differences in distance values ($r$) and argument values ($\theta$) between $\mathbf{A}$ and $\mathbf{B}$. We used four types each of $r$ ($r=1$, 2, 4, and 8 pixels) and four types each of $\theta$ ($\theta =0\mathrm{deg}$, 45 deg, 90 deg, and 135 deg), generated 16 cooccurrence matrices, and used the averages of Haralick feature values calculated by their 16 matrix as $f_{19}$ to $f_{28}$.

The run-length matrix represents the appearance frequency $\mathbf{R}$ $[=R(i,l)]$ of run $l$ in the pixel of interest $i$. Run indicates the number of consecutive identical intensity values in the scanning direction $\theta$. The intensity gradient is frequently subjected to quantization before creating a run-length matrix. We, therefore, used four types each of $\theta$ ($\theta =0\mathrm{deg}$, 45 deg, 90 deg, and 135 deg) and four type each of quantization values (gradations 256, 16, 4, and 2), generated 16 run-length matrix, and used the averages of 5 run-length features calculated by their 16 run-length matrix as $f_{29}$ to $f_{33}$.

Furthermore, Kiyuna et al.¹⁹ previously quantified the complexity of chromatin distribution on a nuclear image from mammary gland cells as a CC value and a fractal feature. These features are represented as $f_{34}$ and $f_{35}$ in the following equations:

Another four feature values had been proposed by Haralick: sum average ($f_{36}$, average of intensity), information measures of correlation 1 ($f_{37}$, uniformity of texture), information measures of correlation 2 ($f_{38}$, diversity of texture), and maximal correlation coefficient ($f_{39}$, uniformity of intensity and texture). These features can be expressed using the following equations:

We designed feature values $f_{1–39}$ as conventional feature (Cf) values.

2.3.

Proposed Feature Values

2.3.1.

Convex hull contour complexity values (Pf.1)

Kiyuna et al.¹⁹ quantified the complexity of chromatin distribution using feature $f_{34}$. However, this feature is not counted if the perimeter $L(i)$ is less than $f_2$; in other words, $f_{34}$ does not consider the chromatin complexities in small regions. Therefore, we previously proposed the following feature value $F$²⁹

Eq. (11)

$$F=\sum _i\left\{1\right[\frac{L(i)}{L_C(i)}\ge 1.2\left]\right\},$$

where $L_C(i)$ and $\frac{L(i)}{L_C(i)}$ represent the convex hull perimeter and convex hull ratio, respectively, of an image binarized using threshold $i$, and the expression in $\mathrm{\Sigma }$ represents the indicator function. $F$ is a variable that counts the number of binarization threshold values $i$, in which the convex hull ratio is 1.2 or more. The value of $F$ increases with chromatin distribution complexity such that $F$ is counted even in small chromatin regions with sufficient complexity. However, $F$ had a large correlation with $f_{10}$, which is the convex hull ratio of the outer shape of the nucleus.

We, therefore, use the following $L_r(i)$, which is a chromatin distribution complexity divided by $f_{10}$, and propose following new feature values $f_{40}$ to $f_{43}$

Eq. (12)

$$L_r(i)=\frac{L(i)}{L_C(i)·f_{10}},$$

Eq. (13)

$$f_{40}=\sum _i\{L_r(i)-1[L_r(i)\ge 1]\},$$

Eq. (14)

$$f_{41}=\sum _i\{1[L_r(i)\ge 1.1]\},$$

Eq. (15)

$$f_{42}=\sum _i\{1[L_r(i)\ge 1.2]\},$$

Eq. (16)

$$f_{43}=\sum _i\{1[L_r(i)\ge 1.3]\}.$$

$f_{40}$ is shown as the fill area in Fig. 2 for which the horizontal axis is the intensity threshold for binarization and the vertical axis is $L_r(i)$. $f_{41}$, $f_{42}$, and $f_{43}$ are shown as intensity widths of the graph when $L_r(i)\ge 1.1$, 1.2, and 1.3 in Fig. 2. For a more detailed representation of the shape of the graph shown in Fig. 2, multiple feature values are used. We name values of $f_{41}$, $f_{42}$, and $f_{43}$ as convex hull (CH) CC, convex hull intensity-width 1.1 (CW1.1), CW1.2, and CW1.3, respectively. We designed these four features $f_{40–43}$ as proposal feature values 1 (Pf.1).

Fig. 2

Convex hull contour complexity values.

2.3.2.

Chromatin distribution spreading value (Pf.2)

We also propose a method for quantifying the chromatin distribution spreading (CDS) value. First, we create an intensity histogram using the gray-scale intensity set from the cell nucleus on the input image and obtain a threshold $T_h$ via a linear discriminant analysis³⁰ of the histogram. $T_h$ is the threshold used to distinguish dark-stained (i.e., assumed chromatin) and light-stained regions (i.e., nonchromatin). Here, a coordinate on the input image is designed as $\mathbf{Z}$ ($\mathbf{Z}$ $=[x,y]$), and the 256-level gray-scale intensity of $\mathbf{Z}$ is designated as $I(\mathbf{Z})$. Next, a chromatin image is generated by replacing the $I(\mathbf{Z})$ of all pixels in the image with $I^{\prime }(\mathbf{Z})$, calculated using the following equation:

Eq. (17)

$$I^{\prime }(\mathbf{Z})=\{\begin{array}{ll}255-I(\mathbf{Z})-T_h,& [I(\mathbf{Z})+T_h<255]\\ 0,& (\mathrm{otherwise})\end{array}.$$

$I^{\prime }(\mathbf{Z})$ is a value obtained by inverting the image negative and positive [$255-I(\mathbf{Z})$] and subtracting the bias value $T_h$. As a result, high-density stained pixels such as nucleoli and chromatin appear as high values. Figures 3(a) and 3(b) show representative chromatin images based on those in Figs. 1(c) and 1(h), respectively, and show at fivefold intensity to enhance visualization.

Fig. 3

Chromatin image. (a) and (b) Image calculated from Fig. 1(c) and 1(h), respectively.

Next, the center of a gravity of the chromatin region $\mathbf{G}$ is obtained using the following equation:

Eq. (18)

$$\mathbf{G}=(\frac{\sum _x\sum _{y}x{I}^{\prime }(\mathbf{Z})}{\sum _x\sum _y{I}^{\prime }(\mathbf{Z})},\frac{\sum _x\sum _{y}y{I}^{\prime }(\mathbf{Z})}{\sum _x\sum _y{I}^{\prime }(\mathbf{Z})}).$$

Finally, we used $\mathbf{G}$ to obtain the CDS value, denoted as $f_{44}$ in the following equation:

Eq. (19)

$$f_{44}=\sqrt{\frac{1}{f_1}\sum _x\sum _y\{{\left|\frac{\mathbf{G}-\mathbf{Z}}{f_3}\right|}^2{I}^{\prime }(\mathbf{Z})\}}.$$

We designed this feature $f_{44}$ as Pf.2.

4.1.

Comparative Experimental Results and Discussion Between Cell Types

To verify the characteristics of our proposed feature values, we calculated feature values $f_{1–46}$ from 633 masking images of cervical smear samples described in Sec. 3. Thereafter, features $f_{1–46}$ were linearly normalized to yield $f_{1–46}^{\prime }$ such that the maximum and minimum values of each feature became 1 and 0. Figure 5 shows an experimental result calculated using the 75% tiles, median values, and 25% tiles of features $f_{1–46}^{\prime }$ for each cell types. However, we note that some Cf values were omitted.

Fig. 5

Experimental results from a feature value analysis of seven cell types.

Figure 5 shows many differences in feature values among NOR, NET, and REG, which we classified as NILM. In particular, REG had a large area ($f_1^{\prime }$ was high), and MET exhibited high texture homogeneity ($f_{28}^{\prime }$ and $f_{39}^{\prime }$ were high, whereas $f_{23}^{\prime }$ and $f_{27}^{\prime }$ were low). All proposed value $f_{40–46}^{\prime }$ were small for NOR and large for SCC. In particular, $f_{44}^{\prime }$ and $f_{46}^{\prime }$ were also large for CIS. Cancer cells possess a chromatin structure that differs from the normal structure.⁵ These results suggest that our proposed method represents the features of this chromatin structure.

In addition, we calculated the absolute $|{\rho }_{i,j}|$ values of the correlation coefficients between $f_i^{\prime }$ and $f_j^{\prime }$ ($i\in \{1,2,\dots ,46\}$, $j\in \{1,2,\dots ,46\}$) using SPSS software (IBM Corporation, Armonk, New York). ${\rho }_{i,j}$ is expressed by the following equation:

Table 2 shows a list of feature numbers ($k$) with high correlation coefficients ($|{\rho }_{i,j}|\ge 0.95$ or $0.95>|{\rho }_{i,j}|\ge 0.90$) with other features. There were strong correlations between size and shape features $f_{1–6}^{\prime }$ and the run-length features $f_{30,33}^{\prime }$. There were also strong correlations among intensity-related features $f_{11,18,24,36}^{\prime }$, all of which are Cfs. In contrast, proposed features $f_{40–46}^{\prime }$ did not show strong correlations with any other features; therefore, our proposed features are highly original. These various feature values are useful for improving the accuracy of machine-learning-based cellular classification, which we will discuss further in Sec. 5.

Table 2

List of feature numbers i and j with high correlation coefficient.

Furthermore, we used SPSS software to perform a one-dimensional ANOVA of the cellular classification corresponding to each feature value. Accordingly, all values ($f_{40–46}^{\prime }$) differed significantly (significance level: 1.0%) among the cellular classifications, and therefore, any proposed or Cf values could potentially improve the accuracy of cellular classification accuracy.

Next, we used a $t$-test to evaluate whether each feature value differed significantly with respect to reactive (MET and REG) and neoplastic nuclear atypia (LSIL, HSIL, CIS, and SCC). The results are shown in Fig. 5(upper): here, the * and ** symbols indicate that the corresponding feature values had significant differences at respective significance levels of 5.0% and 1.0%. This test revealed significant differences in many feature values related to chromatin distribution, including the proposed values $f_{44–46}^{\prime }$. These could, therefore, be considered useful for distinguishing between reactive and neoplastic nuclear atypia.

4.2.

Verification of Experimental Results and Discussion of Representative Images

We next calculated some of the normalized feature values $f_{1,34,35,40–46}^{\prime }$ corresponding to the representative images in Figs. 1(b)–1(f). Figure 6 shows the results, with feature numbers indicated on the horizontal axis.

Fig. 6

Feature values of the images in Fig. 1.

In Fig. 6, many values in area $f_1^{\prime }$, the conventional complexity values $f_{34–35}^{\prime }$, and the proposed complexity values $f_{40–44}^{\prime }$ exhibited similar tendencies; however, in the representative image of SCC, $f_{34–35}^{\prime }$ were moderate, whereas $f_{40–44}^{\prime }$ were high. In addition, the proposed value $f_{45}^{\prime }$ was also high. Although these findings are subjective, the chromatin distribution in Fig. 1(h) appears to be complex and widely spread. We consider that our proposed method reflects this trend.

5.1.

Validation Method

Next, we verified the cellular classification accuracy using machine learning and a CV method.²² These verifications were compared among eight different models, a conventional model (Cf) and seven models combining Cf with models including our proposed values: Cf + Pf.1, Cf + Pf.2, Cf + Pf.3, Cf + Pf.1 + Pf.2, Cf + Pf.1 + Pf.3, Cf + Pf.2 + Pf.3, and Pf.A (= Cf + Pf.1 + Pf.2 + Pf.3).

For machine learning, we used the SVM;^23,24 however, we note that this method is intended for two-class classification. We, therefore, implemented a one-versus-one method²⁵ to expand the classification from two-class to multiclass using a round-robin method of classes. In addition, we selected variables for SVM using a stepwise (floating) method.²⁶ In Sec. 4.1, some of the Cfs showed high correlation coefficients. If we use all of these features directly to create a model of the SVM, the accuracy of identification may decrease due to over-learning.³² The stepwise method we use can mitigate the reduced classification accuracy caused by over-learning, because the possibility of simultaneously selecting features exhibiting high correlation coefficients in the method is low.

Two methods could be used to combine multiclass classification and variable selection. The first involves selecting the same type of features for each comparison, and the second involves selecting different types of features for each comparison. In this paper, we used the second method, which is capable of more detailed feature selection.

A machine learning protocol based on these methods is depicted in Fig. 7(a) as stepwise SVM (SSVM). Before performing the procedure described in Fig. 7(a), we calculated the normalized feature values of $f_1^{\prime }$ to $f_m^{\prime }$ for all 633 cervical smear samples described in Sec. 4 (NOR for 164, MET for 86, REG for 74, LSIL for 36, HSIL for 155, CIS for 84, and SCC for 34). They were a number of imbalanced samples, which can cause incorrect answer rates. We, therefore, virtually matched the sample number of each class using a oversampling method “adaptive synthetic sampling approach for imbalanced learning”³³ to increase the number of each class up to 200 (for a total of 1400 samples) and assumed the value sets to be the feature vectors ${\mathbf{f}}_{1-m}^{\prime }$, where $m$ represents the number of features types used for calculation. For example, $m$ becomes 39 when calculating the machine learning model Cf.

Fig. 7

Flowchart of the (a) SSVM and (b) accuracy evaluation.

In (1) of Fig. 7(a), the order of the normalized feature values ${\mathbf{f}}_{1-m}^{\prime }$ changed randomly with the initialization of some variables. We assumed the changed values to be feature vectors ${\mathbf{F}}_{1-m}$. Here, $\mathbf{C}$ is a set representing the round-robin selection of seven classes, $\mathbf{W}$ is a set representing the kernel functions and cost parameters among the SVM parameters, ${\mathbf{E}}_j=\{{\mathbf{E}}_{j,1},\dots ,{\mathbf{E}}_{j,i},\dots \}$ is a set of feature vectors ${\mathbf{F}}_{1-m}$ selected for model $C_j$, $V_j$ is the SVM parameter selected for model ${\mathbf{C}}_j$, $\mathrm{CV}({\mathbf{E}}_j,C_j,V_j)$ is a function used to calculate the accuracy rate from the CV of the SVM, $A_M$ is the maximum accuracy rate calculated by the CV, and $u$ is an updated flag of the maximum accuracy rate. In addition, LIN, RBF, and the numeric values of the elements of $\mathbf{W}$ in Fig. 7(a), respectively, represent a linear function of the kernel, a radial basis function of the kernel, and the cost parameters included among the SVM parameters. In this paper, we used 10-fold as the number of CV divisions.

In Fig. 7(a), (2), (3), and (4), respectively, represent procedures involving forward feature selection, backward feature selection, and SVM parameter optimization. These values were calculated based on the CV of the SVM. Here, $A_C$ is an accuracy rate calculated by the CV, $\mathbf{E}$ is a power set of $E_j$, and $\mathbf{V}$ is a power set of $V_j$. Feature selection was implemented by calculating these procedures until the maximum accuracy rate was no longer updated, and SSVM was implemented by calculating procedures using the round-robin selection of seven classes.

Figure 7(b) shows an accuracy evaluation procedure based on the SSVM, where ${\mathbf{U}}_k$ is a set of $\mathbf{E}$ and $\mathbf{V}$, SSVM is a function used to calculate the procedure in Fig. 7(a), MCV is a function used to calculate the accuracy rate from the multiclass CV according to the one-versus-one method, $B_{k,h}$ is an accuracy rate obtained from the MCV, $M_k$ is a mean accuracy rate obtained by repeatedly ($50\times$) calculating the MCV, $\mathbf{M}$ is a set of $M_k$ obtained by repeating these procedures, and $d$ is an index number of the maximum value of $\mathbf{M}$. SSVM is likely to fall into a local solution, and ${\mathbf{U}}_k$ is not necessarily the optimum value when obtained from the calculation of a single SSVM. In other words, the selected features and accuracy may be affected by the order of the initial data set. We, therefore, randomly exchanged data sets and extracted the optimum value by repeating the SSVM from $k=1$ to 40 to eliminate the fall into a local solution as much as possible.

Finally, the optimum accuracy rate set ${\mathbf{B}}_d=\{B_{d,1},\dots ,B_{d,k},\dots ,B_{d,50}\}$ and parameter set ${\mathbf{U}}_d$ were outputted, and the results of eight classification models (Cf, Cf + Pf.1, Cf + Pf.2, Cf + Pf.3, Cf + Pf.1 + Pf.2, Cf + Pf.1 + Pf.3, Cf + Pf.2 + Pf.3, and Pf.A) were compared.

5.2.

Validation Results and Discussion

We calculated the averages (Ave.) and standard deviations (SD.) of the accuracy rate set ${\mathbf{B}}_d$ of the eight classification models, using the validation method shown in Sec. 5.1. Table 3 presents the results of a comparison of these values, as well as the Dunnett’s test (D-test) results for each model. Here, D-test 1 represents the D-test results of comparisons between each models and Cf, D-test 2 represents the D-test results of comparisons between each model and Pf.A, and ** indicates a significant difference (significance level = 5.0%). D-test is a multiple comparison, many-to-one procedure (i.e., compares each of many treatment groups with one control group) and is used to verify differences between the average values from each group.^34,35 We used SPSS software to perform this procedure.

Table 3

Experimental accuracy rates for each model.

The average accuracy rates of all proposed models except Cf + Pf.2 were higher than the conventional model (Cf) and exhibited statistically significant differences from Cf by the D-test. Therefore, our proposed models Pf.1 and Pf.3 (features $f_{40–43,45–46}$) are useful features for cervical cell classification by machine learning. Although there was no significant difference between Cf + Pf.2 and Cf, there was a significant difference between Pf.A and Cf + Pf.1 + Pf.3. Therefore, our Pf.2 (feature $f_{44}$) is also a useful feature for cervical cell classification. These results show the usefulness of incorporating our features into the diagnostic support system of the cytology. In addition, these results indicate that our features are different from the Cfs; therefore, our features have the possibility to be useful features in cell diagnosis by the cytologist.

As shown in 5.1, to extend the SVM to a multiclass classification of seven classes, we performed the SVM 21 times in the round-robin selection format; in other words, we obtained 21 selected feature sets $\mathbf{E}$ ($\mathbf{E}\subset {\mathbf{U}}_d$) to create a single machine learning model. We, therefore, extracted the 21 selected feature sets $\mathbf{E}$ of Pf.A, which had the highest accuracy rate and calculated the frequencies as shown in Fig. 8 (cumulative bar chart). Red, magenta, blue, cyan, green, brown, and black colors indicate selected features from comparisons related to NOR, MET, REG, LSIL, HSIL, CIS, and SCC, respectively. Based on Fig. 8, the selection of all proposal features indicates that all contributed to improve the classification accuracy.

Fig. 8

Frequency of selected features in the Pf.A model.