Ensemble Learning-Based Solutions: An Approach for Evaluating Multiple Features in the Context of H&E Histological Images

Abstract

In this paper, we propose an approach based on ensemble learning to classify histology tissues stained with hematoxylin and eosin. The proposal was applied to representative images of colorectal cancer, oral epithelial dysplasia, non-Hodgkin’s lymphoma, and liver tissues (the classification of gender and age from liver tissue samples). The ensemble learning considered multiple combinations of techniques that are commonly used to develop computer-aided diagnosis methods in medical imaging. The feature extraction was defined with different descriptors, exploring the deep learning and handcrafted methods. The deep-learned features were obtained using five different convolutional neural network architectures. The handcrafted features were representatives of fractal techniques (multidimensional and multiscale approaches), Haralick descriptors, and local binary patterns. A two-stage feature selection process (ranking with metaheuristics) was defined to obtain the main combinations of descriptors and, consequently, techniques. Each combination was tested through a rigorous ensemble process, exploring heterogeneous classifiers, such as Random Forest, Support Vector Machine, K-Nearest Neighbors, Logistic Regression, and Naive Bayes. The ensemble learning presented here provided accuracy rates from 90.72% to 100.00% and offered relevant information about the combinations of techniques in multiple histological images and the main features present in the top-performing solutions, using smaller sets of descriptors (limited to a maximum of 53), which involved each ensemble process and solutions that have not yet been explored. The developed methodology, i.e., making the knowledge of each ensemble learning comprehensible to specialists, complements the main contributions of this study to supporting the development of computer-aided diagnosis systems for histological images.

1. Introduction

Histopathological analysis involves procedures that aim to investigate tissue samples that are commonly stained with specific dyes, such as hematoxylin and eosin (H&E) [1,2]. In these processes, specialists identify unusual alterations in the structures of cells and highlight potential abnormal health conditions, such as the diagnosis of cancer. It is noteworthy that this type of disease is a significant cause of early deaths worldwide and has high social and economic costs [3]. For instance, cancer is the second leading cause of death in the United States [4]. Thus, the early detection of diseases often enables less invasive treatments and increases the possibility of finding a cure and/or patient survival. The required steps in the preparation process of H&E images can influence the presentation of histological aspects, further increasing the difficulty of accurately diagnosing diseases under investigation. In addition to this problem, the analysis process takes time and may be susceptible to subjective interpretations by specialists [5,6,7]. These interpretation problems are mainly caused by human issues, such as subjectivity and fatigue. On the other hand, computer-aided diagnosis (CAD) methods play a fundamental role in this task since they can support specialists with second opinions [8,9,10], especially regarding H&E images [5,8,11,12,13,14,15,16,17].

In this regard, two categories of descriptors are typically investigated for the development of CAD systems. The first category consists of handcrafted features (HFs) defined by distinct extraction methods, usually aiming to overcome specific problems [18,19,20,21,22]. Among the HFs, it is possible to highlight techniques based on fractal geometry that use multiscale and multidimensional methods (Higuchi fractal dimension, probabilistic fractal dimension, box fusion fractal dimension, lacunarity and percolation) [23,24,25,26,27,28,29], Haralick [30] and local binary patterns (LBPs) [31]. For instance, Haralick and LBPs have been applied in several imaging contexts [32,33,34], exploring the identification of lung cancer subtypes [35], the presence of cancerous characteristics in breast tissue samples [18,36] and the classification of colorectal cancer [6]. In addition, techniques that involve fractals at multiple scales and/or dimensions have also been applied to quantify the pathological architectures of tumors [23,25,26], demonstrating relevant results in the pattern recognition of prostate cancer [37], lymphomas [38], intraepithelial neoplasia [39], breast tumors [40], colorectal cancer [13] and psoriatic lesions [41]. Moreover, fractal methods are important for texture analysis because they provide information about the complexity of textures and patterns that are similar at various levels of magnification [42].

The second group of descriptors consists of deep-learned (DL) features obtained using convolution neural networks (CNNs) [43]. This group has been useful for defining different CAD approaches [15,17,22,44,45,46,47,48] that consider data representations with multiple levels of abstraction [8,47,49]. The most explored models have provided the best accuracy rates on the ImageNet dataset [50], such as AlexNet [51], GoogleNet [52], VGGNet [53], ResNet [54], Inception [52] and other applications [20,32,55].

Despite the advances provided individually by DL features and HFs, investigations have been carried out to develop models based on combinations of these features [18,21,32,55,56,57,58], generating strategies known as ensemble learning (EL) [20,59]. Moreover, these studies have indicated no single solution for distinct datasets. EL models can also consider distinct classification algorithms in order to obtain more accurate single decisions by applying ensembles of classifiers [20,60]. This type of association has provided important results in the study of cervical cell imaging [20]. Another highlight is the model presented by [59], who conducted a comparative study between a logistic regression classifier trained only with DL features, an ensemble of HFs and an ensemble of all features. The authors concluded that the ensemble involving all features delivered the best distinction rates. In the context of using fractal techniques with CNN models, the method presented by [29] considered an ensemble involving two CNN architectures: one pre-trained with histological images and the other pre-trained with artificial images, which were generated using features from fractal techniques. The authors concluded that their proposal outperformed classification algorithms and CNN models when applied separately.

It is essential to observe that EL models can be developed with the most relevant features, exploring a single selection stage to reduce the search space and increase the accuracy of the system [13,15,18]. The use of feature selection through two stages can also be implemented [61,62,63]. In these cases, the strategies were applied to the ensemble of DL, and the results obtained were better than those obtained via a single stage. Thus, this approach can provide a reduced number of main compositions for developing CAD systems, making knowledge even more comprehensive for specialists. Moreover, the use of this combination is not restricted to image analysis with promising solutions in the frequency domain [64].

The versatility of EL strategies was observed with different combinations between features and classifiers to investigate some types of histological images [20,21,29,32,55,56,59,65,66], but not as described here, which were addressed to define patterns of techniques through multiple H&E datasets. Some examples of EL strategies that can still be investigated are the HFs, DL features, and the HFs and DL features, all of which are in a classifier ensemble context. The best configuration can be compared via classifications only with the use of CNN models, which are useful to indicate the pertinence of using ensemble learning in pattern recognition of various H&E images (colorectal cancer, oral dysplasia, non-Hodgkin’s lymphoma, and liver tissue). Moreover, the previously indicated EL models can also be explored via feature selection in two stages, which is a valuable approach to present more optimized solutions, in addition to significantly reducing a high-dimensional search space, such as those explored here. Thus, known problems such as overfitting or underfitting are minimized [67]. An EL model capable of providing the main solutions for various H&E images, with robust computational approaches for pattern recognition, can significantly improve CAD systems and make knowledge more comprehensive for specialists.

This work presents an EL approach to classify histological images from different contexts. The proposal explored multiple handcrafted features through multidimensional and multiscale fractal techniques (Higuchi fractal dimension, probabilistic fractal dimension, box fusion fractal dimension, lacunarity, and percolation), Haralick and LBPs descriptors, and deep-learned descriptors, which were obtained from several convolutional neural network architectures. Moreover, a two-stage feature selection (ranking with metaheuristic algorithms) with a heterogeneous ensemble of classifiers completed the proposed method to indicate the best solutions. The first stage of selection was defined through the ReliefF algorithm. In the second phase, the approach to discover the most effective features within each reduced subset was employed, exploring particle swarm optimization, genetic algorithm, and binary gray wolf optimization. Each result was verified through a robust ensemble process with a Support Vector Machine, Naive Bayes, Random Forest, Logistic Regression, and K-Nearest Neighbors. This proposal provided the following contributions:

• An EL approach not yet explored in H&E image classification, able to identify the primary combinations of features via two-stage feature selection (ranking with meta-heuristics) with a heterogeneous ensemble of classifiers;
• Best ensembles of descriptors to distinguish multiple histological datasets that have been stained with H&E;
• An analysis of the proposal’s usefulness concerning relevant models available in the specialized literature with indications of the best performances concerning colorectal cancer, oral epithelial dysplasia, and gender classification from liver tissue. This was achieved by utilizing a limited number of features, ranging from 11 to 29 attributes;
• A more robust baseline approach, with solutions without overfitting, which is useful in evaluating and composing new approaches for pattern recognition in histological images;
• A breakdown of the main descriptors present in the best ensembles, making the knowledge comprehensive for specialists and helpful in improving CAD systems.

Section 2 presents the proposed methodology, providing information about the techniques used to compose the ensemble learning approach. Section 3 shows the results and engages in a discussion following the application of this approach. Finally, Section 4 indicates the main findings and suggestions for future exploration.

2. Methodology

In this proposal, the first stage considered techniques for feature extraction. The explored models to compose the handcrafted descriptors included multidimensional and multiscale fractal approaches as well as Haralick and LBPs. The DL features were collected via different CNN architectures with the transfer learning strategy. The explored CNN models were ResNet-50, VGGNet-19, DenseNet-121, Inception v3, and EfficientNet-B2. In the next stage, the descriptors were organized to define the ensembles of features, which were evaluated through a feature selection process based on two stages: ranking with metaheuristic approaches. Lastly, an ensemble of classifiers was determined to indicate the main EL solutions within the scope of various H&E datasets. Figure 1 summarizes the main stages of our proposal. Details are presented in the next sections.

2.1. Handcrafted Fractal Features

Among the various fractal techniques found in the specialized literature, models from a multidimensional and multiscale perspective, such as fractal dimension (probabilistic [24,41], box-merging [68], and Higuchi [12,28,39]), lacunarity [24,41], and percolation [38] were considered in this proposal, since they allow a complementary quantification of color images. Some image types, such as histological medical images, are spectrally and spatially complex and often show certain similarities at different spatial scales. From the chosen approaches, fractal geometry allows the study and description of irregular or fragmented forms of elements in nature as well as complex objects that Euclidean geometry cannot analyze [26]. Moreover, the combinations of fractal features explored here are highly capable of quantifying histological information, such as those existing in the H&E datasets investigated in this study [13,29].

Fractal dimension is often applied to evaluate the irregularity and complexity of a region under analysis, enabling the quantification of the fractional filling of a structure in some scale interval. The lacunarity attribute quantifies the deviation of the translational invariance of a geometric object, indicating how similar the parts of different regions of the objects are to each other. Thus, images with low-value lacunarity are more homogeneous regarding size distribution and the spatial arrangement of gaps. They are also translation invariant, since the sizes of the holes are equal [24,41]. Finally, percolation theory is useful in characterizing many disordered systems, as the percolation process is purely random. Therefore, the locations of the system under analysis to be occupied or that will remain empty are randomly chosen with probability p, allowing the formation of clusters [69]. When a cluster presents a path connecting the two ends of the system, it is considered that there has been percolation. The obtained topology from this process has structures highly related to the fractals. These concepts are also applied in image analysis, considering percolation through connectivity among neighboring pixels [38,70].

2.1.1. Probabilistic Fractal Dimension

The fractal dimension based on a probabilistic approach ($D{F}_p$) was determined following the method outlined in [24,41], which involves the gliding box process. In this technique, for a given color image provided as input, utilizing the RGB color model, each pixel was represented by a 5D vector $(x,y,r,g,b)$ with spatial coordinates $(x,y)$ corresponding to the color components $(r,g,b)$. Using a hypercube with side $L=3$, the image was scanned in steps of one pixel from the top left to the bottom right. After this, the size of the hypercube was incremented by two units until it reached $L_{max}$. An illustration of this process is presented in Figure 2.

Each pixel i of a hypercube, $F_i=f(x_i,y_i,r_i,g_i,b_i)$, was compared to the corresponding central pixel $F_c=f(x_c,y_c,r_c,g_c,b_c)$ using a distance measure. Here, the chessboard (${\Delta }_{ches}$), Euclidean (${\Delta }_{eucl}$) and Manhattan (${\Delta }_{manh}$) distances were considered, according to Equations (1)–(3), respectively. Each pixel$F_i$ with a distance $\Delta$ less than or equal to L was labeled as 1 to indicate that it belongs to the hypercube under analysis. Otherwise, the assigned label was 0. By counting the pixels labeled as 1, it was possible to define a matrix $P(m,L)$, which characterized the probability P that m points belong to a hypercube with side L. The result was a structure according to the illustration available in

Table 1. Illustration of a probability matrix P.

	3	5	…	${\mathit{L}}_{\mathit{max}}$
1	$P(1,3)$	$P(1,5)$	…	$P(1,L_{max})$
2	$P(2,3)$	$P(2,5)$	…	$P(2,L_{max})$
⋮	⋮	⋮	⋱	⋮
$L^2$	$P(L^2,3)$	$P(L^2,5)$	⋮	$P(L^2,L_{max})$

.

$${\Delta }_{ches}=max\left(\right|F_i\left(o_i\right)-F_c\left(o_c\right)\left|\right),o\in r,g,b,$$

(1)

$${\Delta }_{eucl}=\sqrt{\sum _o{(F_i\left(o_i\right)-F_c\left(o_c\right))}^2},o\in r,g,b,$$

(2)

$${\Delta }_{manh}=\sum _o|F_i\left(o_i\right)-F_c\left(o_c\right)|,o\in r,g,b.$$

(3)

Thus, the matrix $P(m,L)$ was normalized according to Equation (4), ensuring that the sum of the elements of a column is equal to 1.

$$\sum _{m=1}^{L^2}P(m,L)=1,\forall L.$$

(4)

From the matrix $P(m,L)$, it was possible to obtain the local fractal dimension $N\left(L\right)$ from each size L through Equation (5). This quantification allowed obtaining the angular coefficient of the linear regression, defined by $log$$L\times log$$N\left(L\right)$, as the probabilistic fractal dimension $D{F}_p$ of the image under analysis.

$$N\left(L\right)=\sum _{m=1}^{L^2}\frac{P(m,L)}{m}.$$

(5)

2.1.2. Box-Merging Fractal Dimension

The fractal dimension with box fusion approach ($D{F}_n$) was applied as described by [68] to quantify a color image within the RGB color space. Each image axis was divided into s partitions to establish a partition table. For instance, the partitions on the x-axis occurred according to Equation (6):

$$t_x=\lfloor \frac{x}{{ϵ}_x}\rfloor =\lfloor \frac{xs}{L_x}\rfloor ,$$

(6)

where $t_x$ indicates the x axis partition, x is the coordinate of any pixel in the box, L is the observation scale, and $ϵ$ refers to the $L/s$ ratio.

The resulting table considered the coordinates of all partitions with at least one element. Identical lines were grouped, resulting in a total of z distinct lines. Subsequently, the regression $log$$z\times log$s was defined, and the corresponding angular coefficient indicated the value of $D{F}_n$ for an image.

2.1.3. Multidimensional and Multiscale Higuchi Fractal Dimension

The fractal dimension, proposed by Higuchi ($D{F}_H$) in 1988 [71], enables the analysis of time series in the 1D domain. The method proved relevant results for biological signal analysis [27,72,73]. The approach explored here considered the approach presented in [28], which expanded the analysis process by including a multidimensional and multiscale strategy.

The applied model considered sets of 1D series from an image I, with dimensions h x w. The multiscale step was defined using different observation scales. The procedure involved sliding vectors to analyze each pixel series. The initial position of the sliding vector was determined from the first position of the pixel series under analysis, and it was incremented by one position to the right until its last position coincided with the last position in the series. The vector lengths were represented as l, such that $3\le l\le min(h,w)$, given the need for each sliding vector to have a central pixel.

In this process, as previously indicated, each pixel of the series under analysis was represented as a 5D vector $(x,y,r,g,b)$ [28,41], allowing an analysis between the pixels involved in each transition of the sliding vector by the series. This procedure characterized the multidimensional approach, which is a useful strategy to quantify the patterns in color images. In each iteration in the multidimensional approach, the central pixel of the sliding vector was compared to every other pixel inside it based on the Manhattan distance (Equation (3)). If $\Delta$ was less than or equal to the current sliding vector length, the central pixel composed an auxiliary subseries referring to the observation scale.

The Higuchi fractal dimension ($D{F}_H)$ was defined through each finite auxiliary subseries extracted from the input series of an image I. In this approach, a finite series of discrete points can be defined as $X=x\left(1\right),x\left(2\right),x\left(3\right),\dots ,x\left(N_{DH}\right)$, such that x indicates an element of the image I and $N_{DH}$ indicates the maximum number of points available in the series under analysis. Thus, the method generated d new series $X_{dj}$, considering that the starting point occurred from j, and d is a granularity factor:

$$X_{dj}:x\left(j\right),x(j+d),x(j+2d),\dots ,x\left(j+\lfloor \frac{N_{DH}-j}{d}\rfloor ·d\right),$$

(7)

where $j=1,2,\dots ,d$.

For each $X_{dj}$, the curve length, $Le{n}_j\left(d\right)$, was calculated such that:

$$Le{n}_j\left(d\right)=\frac{1}{d}\left\{\left(\sum _{i=1}^{\lfloor \frac{N_{DH}-j}{d}\rfloor }\left|x(j+id)-x(j+(i-1)·d)\right|\right)\frac{N_{DH}-1}{\lfloor \frac{N_{DH}-j}{d}\rfloor ·d}\right\}.$$

(8)

Thus, $Len$ was defined as the length of the curve for the interval d, which was calculated as the average value across the d sets of $Le{n}_j\left(d\right)$:

$$Len\left(d\right)=\frac{1}{d}\sum _{j=1}^{d}Le{n}_j\left(d\right).$$

(9)

The $D{F}_H$ value for a subseries was given by the angular coefficient obtained from the linear regression of a $lnLen\left(d\right)\times lnd$ plot, which was obtained through linear least squares fitting. The t auxiliary pixel subseries from each initial input set was employed to calculate the descriptor $D{F}_H$ [28]. Consequently, there are t values of $D{F}_H$. The fractal dimension value of the input series, $D{F}_{HSerie}$, was defined as the average of these t values:

$$D{F}_{HSerie}=\frac{1}{t}\sum _{i=1}^{t}D{F}_{H_i},$$

(10)

where t indicates the number of values that l can take.

The fractal dimension value for the color image was given by averaging the values of all $D{F}_{HSerie}$ [28].

2.1.4. Lacunarity

The lacunarity ($Lac$) was based on the approach from [41,74], considering the same probability matrix indicated previously in Section 2.1.1 to represent the multidimensional and multiscale strategy. This descriptor was obtained through the first-order (Equation (11)) and second-order (Equation (12)) moments based on the distribution measure given by Equation (13).

$$\lambda \left(L\right)=\sum _{m=1}^{L^2}mP(m,L),$$

(11)

$${\lambda }^2\left(L\right)=\sum _{m=1}^{L^2}{m}^{2}P(m,L),$$

(12)

$$\Lambda \left(L\right)=\frac{{\lambda }^2\left(L\right)-{\left(\lambda \left(L\right)\right)}^2}{{\left(\lambda \left(L\right)\right)}^2}.$$

(13)

2.1.5. Multidimensional and Multiscale Percolation

The percolation attribute ($Perc$) was computed using the approach outlined in [38]. Percolation theory was employed to analyze connected pixel paths stretching from one edge to another within an image. The applied method explored a multiscale approach using the gliding box technique. Thus, hypercubes were initially defined with $L=3$ [13,38] were incremented by two units after the complete scan of the image. The relation of the quantity T of hypercubes that traversed an image with height H and width W, as a function of L, was given by the following:

$$T\left(L\right)=(H-L+1)\times (W-L+1),L\le min(H,W).$$

(14)

For each hypercube with size L, a multidimensional approach was applied considering the most relevant color channel, according to the RGB model, aiming to perform a comparison against the central pixel $P_c$, as presented in [24,41] and Section 2.1.1. The comparison was also defined from the three distances mentioned in Equations (1)–(3). Therefore, when the distance $\Delta$ assumed a value less than or equal to L, the pixel P was labeled with $-1$, indicating that the pixel represents a pore.

The percolation clusters were obtained based on the labeling of the Hoshen-Kopelman algorithm, as described in [38]. Once this labeling was given, a cluster was defined by considering neighboring pixels with the same label. When the first cluster has been identified, the algorithm advances to the next unverified pore. From this process, three functions were extracted: the average cluster C; the percolating box ratio Q; and the average coverage ratio of the largest cluster M. The calculation of the average number of clusters per box $C\left(L\right)$ was obtained by the number of clusters in a single box ($c_i$), given a scale L, divided by the total number of boxes T, according to the equation:

$$C\left(L\right)=\frac{{\sum }_{i=1}^{T\left(L\right)}{c}_i}{T\left(L\right)}.$$

(15)

The value of the percolating box ratio Q was obtained by counting the number of percolating boxes for each value of L. A box $q_i$ was counted to increment $Q\left(L\right)$ if the ratio between the number of pixels labeled as pores (${\Omega }_i$) and the total number of pixels ($L^2$) is greater than the percolation threshold p:

$$q_i=\left\{\begin{array}{c}\hfill 1,\frac{{\Omega }_i}{L^2}\ge 0.59275,\hfill \\ \hfill 0,\frac{{\Omega }_i}{L^2}<0.59275.\hfill \end{array}\right.$$

(16)

The ratio of percolating boxes as a function of L ($Q\left(L\right)$) was obtained by dividing the total number of percolating boxes $q_i$ by the total number of boxes T on a scale L:

$$Q\left(L\right)=\frac{{\sum }_{i=1}^{T\left(L\right)}{q}_i}{T\left(L\right)}.$$

(17)

Finally, the average coverage ratio of the largest cluster (M) was calculated by identifying the coverage ratio of the largest cluster in each box evaluated at the L scale, according to Equation (18), where ${\gamma }_i$ indicates the largest cluster in a i box.

$$M\left(L\right)=\frac{{\sum }_{i=1}^{T\left(L\right)}\raisebox{1ex}{${\gamma }_i$}\!\left/ \!\raisebox{-1ex}{$L^2$}\right.}{T\left(L\right)}.$$

(18)

2.1.6. Metrics Obtained from the Descriptor Curves

Fractal descriptors based on the probabilistic approach, lacunarity, and percolation were calculated with L scale variations, according to the gliding box method. In these cases, the value $L_{max}=41$ was considered in this investigation [29]. The multidimensional and multiscale Higuchi fractal dimension was determined with l scale variations of the sliding vector, also with a maximum value of 41, and a granularity factor $d=8$ [28]. The values from each approach were used to define feature curves for each image under analysis. In the Higuchi fractal dimension and lacunarity approaches, the curves were composed of their respective local values concerning each size of the sliding vector/hypercube. Regarding the percolation approach, the obtained curves were C, Q, and M, referring to the percolating regions.

In addition, each curve was represented by scalar values to compose the descriptor vectors [38,75]. We used the following extracted metrics: area under the curve (A), skewness (S), area ratio ($\Gamma$), maximum point ($MP$), and maximum point scale ($MPS$). The area under the curve (A) was calculated by the trapezoidal numerical integration method given by Equation (19).

$$A={\int }_a^{b}f\left(x\right)dx\approx \frac{b-a}{2N}\sum _{n=a}^{b-1}(f\left(x_n\right)+f\left(x_{n+1}\right)),$$

(19)

where a and b denote the minimum and maximum L values, respectively, and N denotes the number of samples.

The skewness (S) metric defines the asymmetry of a sample concerning its mean value. A negative skewness corresponds to a sample where the more frequent values occur on the lower part of the interval range, and the opposite applies if the skewness is positive. A completely symmetric sample yields a skewness of 0. Given a sample with N values, skewness was defined by Equation (20), where $\overline{x}$ denotes the mean value of the sample, and $x_i$ represents the i-th value of x.

$$S={\displaystyle \frac{\frac{1}{N}{\sum }_{i=a}^b{(x_i-\overline{x})}^3}{\sqrt{{\left[\frac{1}{N-1}{\sum }_{i=a}^b{(x_i-\overline{x})}^2\right]}^3}}}.$$

(20)

The area ratio ($\Gamma$) was calculated as the ratio of the area of the right half of the curve to the area of the left half of the curve. This was calculated by considering the area under the curve ($Are{a}_{a,b}$) between two points a and b on the x axis, as represented by Equation (21).

$$Area={\int }_a^{b}f\left(x\right)dx.$$

(21)

From this definition, the area ratio $\Gamma$ was calculated according to Equation (22).

$$\Gamma ={\displaystyle \frac{A_{(\raisebox{1ex}{$b$}\!\left/ \!\raisebox{-1ex}{$2$}\right.+1,b)}}{A_{(a,\raisebox{1ex}{$b$}\!\left/ \!\raisebox{-1ex}{$2$}\right.)}}}$$

(22)

The two last measures were obtained from the maximum point ($MP$) and its respective observation scale L ($MPS$).

2.2. Haralick Features

Haralick descriptors were derived from the co-occurrence matrix [30]. This matrix allowed us to verify the transition between two pixels in gray levels, u and v, considering a distance measure $\delta$ and an arbitrary angle $\theta$. It is important to highlight that in the histology area, there have been studies exploring the interpixel distance variation GLCM textural parameter approach, such as the work proposed by [76]. This study has shown that Haralick descriptors present pattern differences between the cell nuclei for one interpixel distance. Moreover, the one interpixel distance allows the definition of fine-grained texture and regularity and demands reduced computational effort, especially when dealing with large images. Thus, in the proposed framework, the parameter values of distance $\delta =1$ and direction $\theta =0^{\circ },\theta ={45}^{\circ },\theta ={90}^{\circ }$ and $\theta ={135}^{\circ }$ [13,62] were applied. Considering the co-occurrence matrices and the widely known descriptors from [30], the features investigated here were: (1) Angular Second Moment; (2) Contrast; (3) Correlation; (4) Sum of Squares; (5) Inverse Difference Moment; (6) Sum Average; (7) Sum Variance; (8) Sum Entropy; (9) Entropy; (10) Difference Variance; (11) Difference Entropy; (12) and (13) are Information Measures of Correlation; and (14) Maximal Correlation Coefficient. For the final 14 descriptors, the averages of the descriptors obtained in the four matrices were considered.

2.3. Local Binary Patterns

LBPs is used to generate a set of features that represent how binary patterns are distributed in a circular neighborhood around the center pixel. The result is an LBP code for each pixel of the image [31,77]. The neighborhood is defined by a radius R and the number of neighbors $PN$ [78]. If the intensity of the neighboring pixel is greater than or equal to the reference intensity, the binary value 1 is assigned. Otherwise, the value 0 is assigned. Concatenating all the values in the region under analysis, starting from the top left corner and moving clockwise, results in a binary number. The corresponding decimal value takes values between 0 and 255 and represents the analyzed region. Therefore, the LBPs strategy was applied as described in [78,79] with a radius value of $R=1$ and considering $PN=8$ neighbors. From these parameters, it was possible to obtain a total of 59 descriptors for each of the images under analysis, according to the method “extractLBPFeatures(Image)”, which is available in the Matlab tool [80].

2.4. Deep-Learned Features

One of the main problems CNN models face is the number of samples required in the training stage. In this work, this issue was resolved through transfer learning, which is a strategy that uses knowledge obtained in one or more source tasks to improve the learning process in a target task [81,82]. CNN architectures belong to the class of inductive learning algorithms, mapping input features between classes to obtain a generalization of the data. Thus, inductive learning can be transferred from an architecture trained on the source task to the target class by adjusting the model space and correcting the bias values. Usually, this procedure is performed by replacing the last layer of the model (classifier) [32], providing shorter training time and suppressing the possibility of overfitting due to a small number of samples. The different CNN models explored here followed this strategy.

Therefore, the deep-learned descriptors were computed using transfer learning [82], as it allowed analyses involving image sets with a reduced number of samples. Here, the deep-learned features were obtained from CNN architectures pre-trained on the ImageNet dataset [29,32]. Five distinct CNN architectures were considered: ResNet-50 [54], VGG-19 [53], Inception v3 [83], DenseNet-121 [84], and EfficientNet-B2 [85]. Additional information is summarized in Table 2. These networks have shown relevant results in medical image classification problems in varied contexts [32,55,56] and also in histological image classification [20,29,86].

For the ResNet-50 architecture, the deep-learned descriptors were obtained from the last convolutional layer (avgpool) at the point just before the fully connected layer, as explored in other works [32,59,86]. Regarding the VGG-19 architecture, the features resulted from the second fully connected layer (FC2) [87,88]. In the Inception v3 and DenseNet-121 architectures, the descriptors were extracted at the point just before the last fully connected layer [89,90], indicating the last avgpool and global avgpool, respectively. Finally, the deep-learned features explored in the EfficientNet-B2 architecture were extracted from the last fully connected layer (FC) [91].

Table 2. Summary of each CNN model investigated in this study.

CNN Model	Layers	Parameters [92]	Input Dimensions
ResNet-50	50	$25\times {10}^6$	$224\times 224\times 3$
VGG-19	19	$143\times {10}^6$	$224\times 224\times 3$
Inception v3	48	$23\times {10}^6$	$299\times 299\times 3$
DenseNet-121	121	$8\times {10}^6$	$224\times 224\times 3$
EfficientNet-B2	342	$9.2\times {10}^6$	$260\times 260\times 3$

Table 2. Summary of each CNN model investigated in this study.

CNN Model	Layers	Parameters [92]	Input Dimensions
ResNet-50	50	$25\times {10}^6$	$224\times 224\times 3$
VGG-19	19	$143\times {10}^6$	$224\times 224\times 3$
Inception v3	48	$23\times {10}^6$	$299\times 299\times 3$
DenseNet-121	121	$8\times {10}^6$	$224\times 224\times 3$
EfficientNet-B2	342	$9.2\times {10}^6$	$260\times 260\times 3$

2.5. Ensemble of Descriptors

In this study, we explored combinations of handcrafted and deep-learned features through an ensemble of descriptors to broaden the strategies of [21,32,55,56,60,66] with new patterns of techniques in multiple H&E images.

The following frameworks were used: (i) ensemble of handcrafted descriptors with an ensemble of classifiers; (ii) ensemble of deep-learned descriptors with an ensemble of classifiers; (iii) ensemble of handcrafted and deep-learned descriptors with an ensemble of classifiers. Moreover, CNN architectures were directly employed to classify H&E images to validate the effectiveness of the suggested approach within the context of colorectal cancer, oral epithelial dysplasia, non-Hodgkin’s lymphoma, and liver tissues.

The total number of descriptors that comprised the ensembles depended on each category explored here. The handcrafted descriptors tally 462 values, including 389 fractal descriptors distributed as follows: percolation (225), lacunarity (75), probabilistic fractal dimension (63), Higuchi fractal dimension (25), and box-merging approach for fractal dimension (1). Additionally, there are 59 LBPs descriptors and 14 Haralick descriptors. For the deep-learned features, they amount to 10,624 values. This total depends on the layers used in the analyses. Here, as described in Section 2.4, the chosen layers and their corresponding number of descriptors are summarized in

Table 3. Number of extracted descriptors and their respective layer in each architecture.

CNN Architecture	Number of Descriptors	Extracted Layer
ResNet-50	2048	avgpool
VGGNet-19	4096	FC2
Inception v3	2048	last avgpool
DenseNet-121	1024	global avgpool
EfficientNet-B2	1408	FC

.

2.6. Two-Stage Feature Selection and Ensemble of Classifiers

The ensembles of descriptors were evaluated through the application of a two-stage feature selection process. The aim is to reduce the dimensionality of the feature space and identify the most suitable descriptors within each subset [62,93,94]. Initially, our approach involves ranking the feature sets through the ReliefF algorithm [95] and then applying a threshold to decrease the number of potential matches. The threshold considered the 100 best-ranked descriptors in each ensemble [13]. Therefore, the second phase intended to identify the optimal matches within each narrowed subset, which is accomplished through wrapper selection and the exploration of various metaheuristics, namely: genetic algorithm (GA) inspired by genetic evolution, particle swarm optimization (PSO) inspired by particle swarm behavior, and binary gray wolf optimization (bGWO) inspired by the hunting strategy of gray wolves [96]. We used the K-nearest neighbors algorithm to assess the effectiveness of each of these feature selection methods. This proposal allowed the investigation of possible patterns resulting from the combination of techniques in multiple H&E datasets [61,62,63]. Due to the stochastic characteristic of the wrapper selection algorithms, the final performance of a combination was determined from the arithmetic average of 10 different runs of the selection step. Also, the most selected features were defined by their frequencies in the ten runs [97].

2.6.1. Ensemble Classification

The final step of the method focused on evaluating each ensemble of descriptors derived from the preceding step. This assessment was conducted using a strong ensemble of classifiers, exploring five algorithms from distinct categories (based on the function, probability, decision tree, and instance-based K-nearest neighbors). In this case, the ensemble of classifiers was built upon the following algorithms: Support Vector Machine (SVM), K-Nearest Neighbors (KNN), Random Forest, Logistic Regression, and Naive Bayes. The approach allowed combining results from individual classifiers to avoid overfitting [29,59,60]. The classifiers were individually trained, and the ensemble was applied taking into account the sum rule [29,56,60]. Figure 3 illustrates the evaluation process with the ensemble approach. The blocks indicated by dashed lines indicate that only one of the internal blocks was used in the composition of the model.

The assessment process for each classification algorithm involved the use of the cross-validation k-fold method, with $k=3$ [62]. In this approach, the instances were divided into k independent groups, with each group further divided into a training set and a test set. In each iteration of the training process for k, one fold was used for evaluation, while the remaining $k-1$ folds served as training data. Instances within each group were classified, resulting in k performance measures, one for each fold. The ultimate performance was determined by calculating the arithmetic average of these partial performances. The metrics used to evaluate each solution were the area under the ROC curve (AUC) and accuracy [98].

2.7. Application Context

The proposed approach was employed to assess five histological image datasets, namely colorectal cancer (CR), oral epithelial dysplasia (OED), non-Hodgkin’s lymphoma (NHL), and liver tissues, covering liver gender (LG) and liver aging (LA).

The CR dataset consists of histological images extracted from 16 sections stained with H&E, specifically focusing on T3 or T4 stage colorectal cancer. An expert in pathology labeled the images, categorizing them into benign or malignant groups. The dataset includes 165 images, with 74 images representing benign cases and 91 images representing malignant cases [99].

The OED dataset was created using 30 H&E histological sections from the tongues of mice that had been previously exposed to a carcinogen [100]. This study was approved by the Ethics Committee, protocol 038/39 of the Federal University of Uberlândia. A total of 148 ROIs divided into two classes, healthy (74) and severe dysplasia (74), were considered in this work.

The NHL dataset consists of 374 H&E images of malignant non-Hodgkin’s lymphoma subdivided into three classes: Chronic Lymphocytic Leukemia (CLL) with 113 images; Follicular Lymphoma (FL) with 139 images; and Mantle Cell Lymphoma (MCL) with 122 images. The dataset is a collection of samples prepared by different histologists at various hospitals [101].

The last two datasets comprise liver tissues from 48 male and female mice obtained through the National Institute on Aging, Atlas of Gene Expression in Mouse Aging Project (AGEMAP). Images were acquired manually by a Carl Zeiss Axiovert 200 microscope at 40× magnification. The liver gender (LG) dataset consists of 265 images, which were categorized into male (150) and female (115) classes of mice under an ad libitum diet. The liver aging (LA) dataset, on the other hand, contains 528 images of female mice at four ages: 1 (100), 6 (115), 16 (162) and 24 (151) months, also on an ad libitum diet [102].

An overview of the datasets is given in

Table 4. General information about the datasets used in this investigation.

Dataset	Image	Number of Images	Number of Classes	Resolution
CR	Colorectal Tumor	165	2	567 × 430 to 775 × 522
OED	Oral Epithelial Dysplasia	148	2	450 × 250
LA	Liver Tissue	528	4	417 × 312
LG	Liver Tissue	265	2	417 × 312
NHL	Non-Hodgkin’s Lymphoma	374	3	1388 × 1040

, indicating the names and their respective types, number of images, and number of classes. Some examples of images with their respective classes are presented in Figure 4.

2.8. Software Packages and Execution Environment

In this work, the handcrafted descriptors were implemented and extracted via MATLAB^® R2019a [103]. The pre-trained convolutional neural network architectures were implemented via the Pytorch v1.10 library [104] in the Google Collaboratory execution environment [105]. Considering the CNN models, the package explored here uses a strategy in the first layer that allows the input of images of any size. This strategy avoided the loss of significant information from the images. Also, the feature selection and classification algorithms were employed using the Weka 3.8 [106] software. The default values suggested in each package were used except for specific ones mentioned in the text. The experiments were performed on an Intel Core i5 notebook, 8 GB RAM and 64-bit architecture operating system.

3. Results and Discussion

The ensemble learning approach was evaluated on five datasets of H&E histological images, as described in Section 2, with comparisons involving the different classes of each set. It is important to note that 45 types of tests were performed to explore different compositions of ensembles, including three associations of wrapper methods, in order to provide the main compositions among the 100 best-ranked descriptors with ReliefF. Each composition was evaluated via a heterogeneous ensemble of classifiers (Section 2.6.1). In

Table 5. Average accuracy values obtained through the proposed approach with two-stage feature selection (ReliefF + Wrapper).

Ensemble	Wrapper	CR	OED	LA	LG	NHL
HFs	PSO	91.94%	100%	98.30%	99.36%	78.02%
	GA	91.39%	99.66%	98.58%	99.36%	77.59%
	bGWO	92.30%	99.66%	98.73%	99.36%	77.94%
DL	PSO	99.45%	96.89%	93.50%	98.75%	90.64%
	GA	99.15%	97.64%	93.01%	98.75%	90.56%
	bGWO	99.76%	96.96%	93.39%	98.87%	90.13%
HFs+DL	PSO	99.58%	98.45%	96.86%	99.40%	90.51%
	GA	99.45%	98.45%	96.72%	99.36%	90.72%
	bGWO	99.45%	99.05%	96.44%	99.47%	90.43%

The best accuracy rates for each dataset are highlighted in bold.

, the average performances for each ensemble, considering the HFs and DL attributes, are shown. The best rates are highlighted in bold.

From Table 5, it is observed that the HFs and HFs+DL ensembles were responsible for the best results in four H&E datasets (OED, LA, LG, and NHL) out of the five investigated here. The accuracy values ranged from 90.72% to 100%. Thus, it is possible to indicate that the handcrafted descriptors explored here (via HFs ensemble) are relevant for the classification process whether used separately or in combination with DL. The HFs ensemble provided the highest distinction rate on the OED dataset regardless of the wrapper selector taken as reference, indicating the optimal match on this dataset, which represents a further contribution of the proposed approach. On the other hand, HFs presented the lowest accuracy values (approximately 78%) in the NHL dataset, with three classes (CLL×FL×MCL), demonstrating a possible limit of HFs. When this category was combined with DL (HFs+DL ensemble), the result was the most expressive for the NHL dataset with an accuracy of 90.72%. Even so, this result represents an accuracy of at least 8% lower than the values achieved on other datasets. This is an important indication of the difficulties in distinguishing CLL×FL×MCL groups, especially considering only the HFs ensemble. Finally, when the DL ensemble is considered, the best solution was achieved in a single dataset (CR) but with an expressive rate (99.76% accuracy), illustrating its importance for the development of strategies to support the diagnosis of colorectal cancer. In addition, on this dataset, it is worth mentioning the HFs+DL ensemble as another potential solution, which achieved an accuracy of 99.58%, which was very close to that provided by the DL combination. This configuration represents an acceptable and common solution for different types of histological samples.

To summarize the results discussed here, the best combinations of descriptors and selection algorithms are presented in

Table 6. Main combinations in each H&E-stained histological dataset based on the criterion of higher accuracy with the lowest number of descriptors.

Dataset	Ensemble	Feature Selection	⌈Average of Descriptors⌉	Average Accuracy	Average AUC	Top 1 Accuracy	Top 10 Accuracy
CR	DL	ReliefF + bGWO	21 ± 3	99.76% ± 0.30	1	100%	99.39%
OED	HFs	ReliefF + PSO	11 ± 3	100%	1	100%	100%
LA	HFs	ReliefF + bGWO	40 ± 7	98.73% ± 0.36	0.999 ± 0.4 × ${10}^{-3}$	99.24%	98.11%
LG	HFs+DL	ReliefF + bGWO	29 ± 6	99.47% ± 0.30	1	100%	98.87%
NHL	HFs+DL	ReliefF + GA	53 ± 5	90.72% ± 1.04	0.980 ± 1.5 × ${10}^{-3}$	92.25%	89.57%

, including the total number of descriptors, AUC, and accuracy averages. The accuracy values of the top 1 and top 10 solutions are also indicated, making it possible to observe the existing variation for the first and tenth solutions in each dataset, since the averages were calculated from the 10 best-ranked compositions in each dataset. It is important to emphasize that this ranking indicated the highest accuracy with the lowest number of descriptors.

Based on the previously stated criterion, it can be reiterated that the HFs and HFs+DL ensembles were responsible for the best results in four of the five H&E-stained histological datasets investigated here. In these cases, the ReliefF + bGWO selection processes stand out with three occurrences. This indicates another pattern for the CR, LA, and LG datasets, with expressive average AUC rates of 0.999 (LA) and 1 (CR and LG). In addition, the lowest top 10 accuracy was significant, with a rate above 98% (LA), on the dataset that involves four groups for the classification process. Moreover, the main solutions for these three datasets indicated a reduced number of descriptors, an average of 21 features for CR, 29 for LG, and 40 for LA.

Regarding the two-stage ReliefF + PSO, this strategy was the main solution for the OED dataset, providing maximum accuracy with the lowest number of descriptors among all solutions with only 11 features. When the two-stage ReliefF + GA is observed, this approach constitutes the best solution on a single dataset (NHL). In this case, the solutions explored 53 features on average, identifying the highest value among all solutions. The top 10 accuracy was 89.57%, and the top 1 was slightly better, 92.25%, reinforcing the difficulties present in this set. The NHL dataset comprises three classes (CLL×FL×MCL) and, possibly, with less heterogeneous histological patterns, implying more difficulties in the constitutions of the solutions. Even so, in this case, the average AUC was 0.98, which is an important value under the exposed conditions. For instance, considering that the original feature space had a range of 462 to 11,086 values, the outcome achieved in this study is another contribution that is capable of providing expressive average performances with few descriptors but highly relevant to the classification process.

3.1. Feature Occurrences in the Main Solutions: An Overview

To identify the descriptors present in the top solutions, as summarized in Table 6, a survey of the occurrence of each category of features in the first 10 solutions of each H&E-stained histological dataset was performed. To better understand the origin of each descriptor in each solution, the occurrences of the deep-learned descriptors are in Figure 5, CR, LG, and NHL datasets, and the handcrafted ones are in Figure 6, indicating solutions for LA, LG, OED, and NHL. It is important to highlight that the best solution for the LA dataset involved only HFs. Occurrences in NHL and LG also involved HFs due to the HFs+DL ensemble, justifying the representation of these datasets in Figure 6. Also, in these two datasets, the occurrence percentages were calculated relative to the total number of HFs in the HFs+DL ensemble, disregarding the percentages of deep-learned features.

Considering the distributions illustrated in Figure 5, it is possible to verify some behaviors. The lowest occurrences occurred concerning the Inception v3 network descriptors, a maximum of 3.37% for the CR dataset, and there were no instances among the best 10 solutions for the LG set. On the other hand, the descriptors via DenseNet-121 and EfficientNet-B2 networks have the highest occurrences, especially for the NHL dataset, in which 63.45% of the features originated from the DenseNet-121 model. Descriptors via the EfficientNet-B2 architecture stood out in the solutions for the CR dataset with 38.46% of the occurrences surpassing more homogeneous occurrences (from 17.79% to 21.63%) for the deep-learned features from the ResNet-50, VGG-19, and DenseNet-121. Another homogeneous distribution can be seen in the LG dataset, involving the same descriptor origins as the CR set. In this case, occurrences ranged from 15.17% to 23.45%. When DL versus HFs totalization is considered, it can be seen that DL attributes predominated in the solutions for the LG and NHL datasets with occurrences of 65.17% and 97.16%, respectively. Despite these differences, it is not possible to indicate that these were the most important in the classifications.

Concerning the occurrences of handcrafted descriptors (Figure 6), the lowest occurrence was of the box-merging fractal dimension ($D{F}_n$ from [68]), since it was not selected for the top 10 solutions in three out of four histological datasets stained with H&E. This descriptor was present in the solutions for OED but with the lowest occurrence, only 1.84%. The probabilistic fractal dimension descriptor ($D{F}_p$ from [24,41]) was the second lowest occurrence but constituted the solutions for three of the four histological datasets. Another interesting result involves the enhanced version of Higuchi fractal dimension descriptor ($D{F}_H$ [28]) with occurrences that surpassed those of the $D{F}_n$ and $D{F}_p$ approaches, which are widely explored in the literature [24,41,68], contributing to advancements in this particular research field. Finally, it is possible to observe the descriptors with the highest occurrences for each H&E dataset: lacunarity (45.80%) for LA; percolation (46.36%) for OED; LBPs (39.60%) for LG; and with a highlight, Haralick as the only ones that constituted the solutions for NHL.

3.2. Performance Overview against Different Approaches

The best performances were observed respecting those obtained via traditional CNN architectures applied directly to the H&E-stained histological images. ResNet-50 [54], VGG-19 [53], Inception v3 [83], DenseNet-121 [84] and Efficient-Net-B2 [85] were the models tested in this process, using the following: a fine-tuning process; cross-validation k-folds, with $k=3$; 10 epochs; the stochastic gradient descent algorithm; initial learning rate of 0.01, decaying by 0.75 every two epochs; and, loss function as cross-entropy. Similar experiments were described by [29]. Also, the CNN models were applied using the transfer learning strategy [82], considering each network pre-trained on the ImageNet dataset [29,32]. Thus, each dataset with a reduced number of examples was investigated using each model after a fine-tuned process, mapping the last corresponding layer of each architecture with the groups available in each H&E dataset. The final connections with their weights were updated based on the total number of classes in each context, ensuring appropriate results without overfitting. In addition, the input images were normalized according to the mean and standard deviation values of the ImageNet dataset [107]. The accuracy values the networks provide are shown in

Table 7. Accuracy provided by CNN architectures after fine-tuning process: best results are in bold.

CNN Architecture	CR	OED	LA	LG	NHL
ResNet-50	96.73%	96.00%	91.48%	97.78%	74.27%
VGG-19	98.67%	92.67%	77.96%	91.98%	65.07%
Inception v3	85.64%	90.50%	85.37%	92.59%	67.07%
DenseNet-121	97.45%	93.17%	87.41%	93.98%	72.40%
EfficientNet-B2	96.12%	94.67%	90.00%	98.89%	70.27%

The best accuracy rates for each dataset are highlighted in bold.

.

To understand the differences between the distinction rates of the models proposed here and those obtained via networks applied directly, we consider the average accuracy values achieved in each H&E-stained histological dataset, as summarized in Table 6. Thus, it is possible to verify that the accuracy values via the proposed approach overcome those provided by the ResNet-50, VGG-19, Inception v3, DenseNet-121, and EfficientNet-B2 networks. The classification rates with the convolutional networks ranged from 74.27% to 98.89%. Therefore, the gains in accuracy ranged from 0.6% to 16%, approximately. The smallest gain occurred in the CR dataset (0.58%) and the largest (16.45%) occurred in the NHL set. For example, we increased the classification rate in the NHL dataset, which involves three classes, from 74.27% to 90.72%, illustrating an additional contribution of this study.

Regarding the noted differences, the Friedman test was applied to verify if these solutions are statistically relevant. The Friedman test is a non-parametric statistical method capable of ranking the solutions under investigation, where the best option is set at the first position [108]. This type of test allows us to observe the variance of repeated measures and to analyze whether the existing differences are statistically significant via p-values. The smaller the p-value, the greater evidence that the difference is statistically relevant. It is possible to indicate that there is some relevant difference when the p-value is less than $0.05$. In the experiments carried out here, the resulting p-value was 0.0004, indicating that the differences between the solutions are statistically significant.

In addition, the Friedman test ranks the solutions as a table. The result involving the experiments is displayed in

Table 8. Friedman test: ranking of methods.

Ranking	Method	Friedman’s Score
1	HFs+DL + (ReliefF+bGWO)	2.40
1	HFs+DL + (ReliefF+GA)	2.40
2	HFs + (ReliefF+PSO)	3.60
2	HFs + (ReliefF+bGWO)	3.60
3	DL + (ReliefF+bGWO)	4.00
4	ResNet-50	6.20
5	EfficientNet-B2	6.80
6	DenseNet-121	7.20
7	VGG-19	8.60
8	Inception v3	9.40

, with Friedman’s score indicated. The solutions obtained in this study are the most relevant for each dataset. It is important to note that when applying the Friedman test, each dataset represents a different sample (row) in relation to the corresponding solution. Each performance obtained through a solution in an H&E set has a rank value assigned based on the order of the best solutions. In the case of a tie, average ranks were assigned to the solutions. In each column, Friedman’s score was calculated as the average of the ranks of the samples, providing a final score for each solution [108].

Finally, we believe that the heterogeneous ensemble of classifiers was another relevant factor in achieving the results listed previously. We were able to define a combination of algorithms that supported the pattern recognition process of different types of H&E images with a more robust and reliable system capable of covering the weaknesses that may exist in a single classifier. In addition, we believe that the bias and variance have been reduced, minimizing the overfitting. More comparisons or algorithms could be considered to indicate the possible limits of each solution or even whether the main combinations are maintained from more descriptors or selection methods. However, the set of techniques with their associations and experiments described here provided an important overview of the potential and discriminative capacity regarding H&E-stained histology images.

Observations Based on Related Work

Considering the best results (Table 6), an overview is presented in relation to the literature involving each H&E-stained histological dataset. It is essential to note that even though the images share the same type, the datasets are not similar in the number of examples, classes, metrics used, validations, or distinct samples. Therefore, the purpose of this validation is not to compare performance rates directly. This situation requires equal conditions between the models: for example, by reproducing each strategy in the H&E-stained histological datasets explored here. This task is difficult or even infeasible. Hence, it is expected to observe whether the obtained solutions provide results among those available in the specialized literature. This illustrative overview is displayed in Table 9, Table 10, Table 11 and Table 12 for the CR, OED, LA and LG, and NHL datasets, respectively.

From the values collected for this type of observation, it is verified that the solutions obtained here could provide highly competitive accuracy values, especially in three (CR, OED, and LG) of the five H&E datasets. In relation to the LA dataset (Table 11), the proposed approach provided an accuracy very close to the best solutions available in the literature: 99.24% (our solution) against 99.62% [29], which is a difference of only 0.38%. Taking into account the NHL dataset, even though the obtained solution indicated an accuracy among those achieved in the proposals of [62,109], a better-defined difference is verified in relation to the results of [20,29,56], showing that the strategies explored here deserve more attention for this type of image. Nevertheless, this overview has shown the discriminative capacity of the solutions obtained for different types of histological images, considering strategies not yet explored in the literature. Moreover, the ensemble learning approach provided relevant solutions with important information about the best combinations of descriptors and selection methods, using a reduced set of descriptors and revealing their main occurrences to recognize possible patterns in four types of histology tissues.

Table 9. Accuracy values provided by different strategies for CR image classification.

Author	Method	Accuracy
Proposed Method	DL + (ReliefF + bGWO)	100 %
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	99.39%
[15]	ResNet50 (activation_48_relu layer), ReliefF and 35 deep-learned features	98.00%
[56]	8 CNN models, handcrafted descriptors	97.60%
[20]	9 CNN models, handcrafted descriptors	97.50%
[62]	Le-Net, multidimensional and multiscale fractal features, Haralick and LBPs	91.06%
[14]	ResNet deep-tuning (DL)	86.67%

Table 9. Accuracy values provided by different strategies for CR image classification.

Author	Method	Accuracy
Proposed Method	DL + (ReliefF + bGWO)	100 %
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	99.39%
[15]	ResNet50 (activation_48_relu layer), ReliefF and 35 deep-learned features	98.00%
[56]	8 CNN models, handcrafted descriptors	97.60%
[20]	9 CNN models, handcrafted descriptors	97.50%
[62]	Le-Net, multidimensional and multiscale fractal features, Haralick and LBPs	91.06%
[14]	ResNet deep-tuning (DL)	86.67%

Table 10. Accuracy rates achieved in different methods for OED image classification.

Author	Method	Accuracy
Proposed Method	HFs + (ReliefF + PSO)	100%
[17]	OralNet: Fused Optimal Deep Features	99.50%
[110]	Neural architecture search and handcrafted descriptors (morphological and non-morphological)	95.20%
[111]	Handcrafted descriptors (SIFT, SURF, ORB)	92.80%
[100]	Handcrafted descriptors (morphological and non-morphological)	92.40%
[16]	Densenet121	91.91%

Table 10. Accuracy rates achieved in different methods for OED image classification.

Author	Method	Accuracy
Proposed Method	HFs + (ReliefF + PSO)	100%
[17]	OralNet: Fused Optimal Deep Features	99.50%
[110]	Neural architecture search and handcrafted descriptors (morphological and non-morphological)	95.20%
[111]	Handcrafted descriptors (SIFT, SURF, ORB)	92.80%
[100]	Handcrafted descriptors (morphological and non-morphological)	92.40%
[16]	Densenet121	91.91%

Table 11. Accuracy rates indicated in different approaches for distinguishing LA and LG images.

Author	Method	Accuracy (LA)	Accuracy (LG)
Proposed Method	LA: HFs + (ReliefF + bGWO); LG: HFs+DL + (ReliefF + bGWO)	99.24%	100%
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	99.62%	99.62%
[15]	ResNet50 (activation_48_relu layer), ReliefF and 5 deep-learned features	–	99.32%
[57]	Inception-V3, Fractal Dimension and Lacunarity (DL+HFs)	–	99.25%
[109]	CNN for texture	99.10%	98.20%
[112]	GIST handcrafted descriptor	88.40%	93.70%

Table 11. Accuracy rates indicated in different approaches for distinguishing LA and LG images.

Author	Method	Accuracy (LA)	Accuracy (LG)
Proposed Method	LA: HFs + (ReliefF + bGWO); LG: HFs+DL + (ReliefF + bGWO)	99.24%	100%
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	99.62%	99.62%
[15]	ResNet50 (activation_48_relu layer), ReliefF and 5 deep-learned features	–	99.32%
[57]	Inception-V3, Fractal Dimension and Lacunarity (DL+HFs)	–	99.25%
[109]	CNN for texture	99.10%	98.20%
[112]	GIST handcrafted descriptor	88.40%	93.70%

Table 12. Accuracy values defined by different approaches for NHL image classification.

Author	Method	Accuracy
[20]	9 CNN models, handcrafted descriptors	97.33%
[56]	8 CNN models, handcrafted descriptors	97.33%
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	95.55%
Proposed Method	HFs+DL + (ReliefF + GA)	92.25%
[62]	Le-Net, multidimensional and multiscale fractal features Haralick, LBPs	82.01%
[109]	CNN for texture	65.10%

Table 12. Accuracy values defined by different approaches for NHL image classification.

Author	Method	Accuracy
[20]	9 CNN models, handcrafted descriptors	97.33%
[56]	8 CNN models, handcrafted descriptors	97.33%
[29]	ResNet50 with fine-tuning, multidimensional and multiscale fractal features	95.55%
Proposed Method	HFs+DL + (ReliefF + GA)	92.25%
[62]	Le-Net, multidimensional and multiscale fractal features Haralick, LBPs	82.01%
[109]	CNN for texture	65.10%

4. Conclusions

In this work, an ensemble learning method was elaborated through multiple descriptors (handcrafted and deep-learned features), a two-stage feature selection, and a classification process with five algorithms (heterogeneous ensemble). The approach was utilized to categorize H&E histological images that are representative of various datasets, such as colorectal cancer, liver tissue, oral dysplasia, and non-Hodgkin’s lymphomas.

The best ensembles indicated average accuracy values ranging from 90.72% (NHL) to 100% (CR). Since the initial feature set was composed of 11,086 values (462 handcrafted descriptors and 10,624 deep-learned features), the best solutions used a maximum of 53 features, with the following scenarios being noteworthy: CR with only 21 descriptors via bGWO; OED with only 11 descriptors via PSO; LA with 40 descriptors via bGWO; LG with only 29 attributes, through the bGWO; NHL with 53 descriptors, via GA. A breakdown of the main descriptors was also presented. It was observed that deep-learned descriptors predominated in relation to handcrafted ones, especially in the solutions for the LG and NHL datasets, with occurrences of 65.17% and 97.16%, respectively. On the other hand, the best solution for the LA dataset involved only handcrafted attributes. Another interesting behavior regarding handcrafted attributes is that the improved version of the Higuchi method outperformed the occurrences of important fractal techniques, specifically $D{F}_n$ and $D{F}_p$, indicating the potential of the descriptor in multiple H&E-stained histological datasets. In addition, the handcrafted features with the highest occurrences were lacunarity (45.80%, LA dataset), percolation (46.36%, OED dataset), LBPs (39.60%, LG dataset) and Haralick (100%, NHL dataset). The indications of solutions, attributes, and occurrences represent important contributions of this study, since the composition of each model and the conditions involved are available to specialists interested in these issues.

When comparing the optimal outcomes with those achieved through CNN architectures applied directly to the H&E-stained histological datasets, it is noted that the proposed approach presented a superior performance in all conditions explored here. Moreover, regarding the performances available in the specialized literature for the same image contexts, the proposal provided the best solutions in three (CR, OED, and LG) of the five datasets, exploring from 11 (OED) to 29 (LG) features. Therefore, these results confirm the proposal as a robust baseline approach capable of providing models without overfitting, offering valuable insights for the assessment and enhancement of CAD systems tailored explicitly for H&E samples, particularly those representing CR, OED, and LG.

Finally, some issues concerning the proposed approach deserve attention. For instance, the effectiveness of parameter tuning, algorithm inclusion, and attribute selection methods may heavily depend on the dataset explored. The solutions may not generalize well to other types of histological images. Also, the success of applying metaheuristics and other algorithms relies on their suitability for the given problem. Biases might arise if specific algorithms are more effective due to the nature of the data, potentially favoring certain types of classifiers. Finally, the use of cutoff points for attribute selection via the ReliefF algorithm introduces a subjective element. The chosen cutoff points could impact the definition of best attributes, leading to potential biases based on the selected thresholds.

In future work, we intend to investigate the following: the limits and impacts on the best ensembles after applying parameter tuning methods for metaheuristics, including other algorithms; a scheme that aims to understand why the features were selected, in addition to which of them are most important for the classification process; influences of cutoff points to define the best attributes via the ReliefF algorithm (first stage of selection); the discriminative power of handcrafted attributes and corresponding ensembles based on quantifications of explainable artificial intelligence representations, specifically gradient-weighted class activation mapping and locally interpretable model-agnostic explanations; the discriminative capacity of these combinations and conditions in other H&E-stained histological images; comparisons of the main results with other existing methods or algorithms commonly used in the analysis of histological images.