Abstract
In recent years, discrete orthogonal moments have attracted the attention of the scientific community because they are a suitable tool for feature extraction. However, the numerical instability that arises because of the computation of high-order moments is the main drawback that limits their wider application. In this article, we propose an image classification method that avoids numerical errors based on discrete Shmaliy moments, which are a new family of moments derived from Shmaliy polynomials. Shmaliy polynomials have two important characteristics: one-parameter definition that implies a simpler definition than popular polynomial bases such as Krawtchouk, Hahn, and Racah; a linear weight function that eases the computation of the polynomial coefficients. We use IICBU-2008 database to validate our proposal and include Tchebichef and Krawtchouk moments for comparison purposes. The experiments are carried out through 5-fold cross-validation, and the results are computed using random forest, support vector machines, naïve Bayes, and k-nearest neighbors classifiers.
1. Introduction
Over the past few years, discrete orthogonal moments (DOMs) have attracted the attention of the image analysis community because they possess the attribute of describing local and global features in images efficiently. The DOMs are optimal in the sense that they represent images with minimal information redundancy. The characteristic is originated by the orthogonality condition of the polynomial basis where the moments are computed [1]. Furthermore, contrary to the continuous moments, the DOMs are defined in the discrete domain. Hence, their computation does not require spatial quantization [2].
Families of DOMs such as Racah [3], dual Hahn [4], and Tchebichef [5] have been used successfully in many applications, for instance, feature analysis [6, 7], face recognition [8, 9], and image retrieval [10, 11] to name a few. Particularly, medical imaging has exploited orthogonal moments to characterize different types of biological tissue and assess the severity of several diseases [12–16].
An important impulse that has placed them in the spotlight is the discussion of the numerical instability problem that occurs in high-order moments. Since the DOMs are calculated as the projection of the image on a weighted kernel or set of orthogonal polynomials, their computation is usually linked to the size of the image. The higher the order of the moment, the higher the numerical error.
Traditionally, the computation of the DOMs is carried out by the use of recursive equations. However, the methodology tends to accumulate and propagate errors that degrade the representation of the image. Other methodologies have tackled the issue. In [17], Mukundan presented a seminal paper where the symmetry and renormalization of the Tchebichef polynomials are used to reduce the accumulation of numerical errors; the discrete Krawtchouk moments [18] are defined on a set of weighted Krawtchouk polynomials to improve the numerical stability; and Bayraktar et al. [19] proposed computing an offline lookup table of Tchebichef and Krawtchouk polynomial coefficients to overcome calculation errors.
In addition to the aforementioned normalization methods, other approaches such as partitioning and sliding window have been also proposed. The former divides the image into smaller nonoverlapped sub-images where low-order moments are calculated independently, while the latter uses a rectangular window that slides usually from left to right and from top to down. On every window region low-order DOMs are computed.
A new class of discrete orthogonal moments derived from Shmaliy polynomials [20] was recently proposed. Shakibaei Asli and Flusser called this new class of moments “discrete Shmaliy moments” (DSMs) [21]. We highlight that the Shmaliy polynomial basis has two important characteristics: one-parameter definition that implies a simpler definition than popular polynomial bases such as Krawtchouk, Hahn, Dual Hahn, and Racah; a linear weight function that eases the computation of the polynomial coefficients. The Shmaliy polynomials are opposite to the Tchebichef polynomials, which use a symmetric nonlinear weight function.
So far, the experiments conducted by Shakibaei Asli and Flusser have shown the capability of the DSMs as feature descriptors in one dimension. In this paper, we explore the use of the discrete Shmaliy moments as 2D texture descriptors and present an extensive comparative study using the IICBU-2008 database [22]. A quantitative analysis based on N-way ANOVA is also conducted to support our proposal.
In Section 2, the mathematical theory of DOMs is presented. We introduce DSMs as the main part of this paper. Tchebichef and Krawtchouk moments are also briefly described for comparison purposes. The rest of this paper is organized as follows: in Section 3, we review the texture model and the general classification scheme; the image database is described in Section 4; and in Section 5, the experimental results are summarized. Finally, Section 6 contains concluding statements.
2. Discrete Shmaliy Moments
In 1980, Teague defined the orthogonal moments in the continuous space [23]. Since then, many variations have been published and used successfully. For example, Zernike and Legendre moments have been used as detectors of invariant features [24–26]. However, the implementation of a continuous basis involves a discrete approximation that affects the properties of invariance and orthogonality. In addition, in some cases, the coordinated space of the image must be transformed in order to compute the continuous moments [27, 28].
Generally speaking, the orthogonal moments are a set of scalar quantities that are not correlated among them. They are useful for characterizing local and global features in images [29]. The orthogonal moments are computed as the projection of the image on an orthogonal polynomial basis . Formally, they are defined as follows: . One way of interpreting the projection is as the correlation measure between the image and the polynomial basis [15]. On the other hand, the discrete orthogonal moments overcome the aforementioned limitations because they still satisfy the property of orthogonality and are also defined in the same image domain.
Discrete Shmaliy polynomials were developed as unbiased finite impulse responses for predictive FIR filters [30] and have been applied on blind fitting in finite-length data [20]. The discrete Shmaliy polynomials have two important characteristics. Their definition is simple because it requires only one parameter that represents the length of the data, and they use a linear weight function, which eases the computation. In addition, the discrete Shmaliy polynomials represent the discrete version of the radial Mellin polynomials when the length of the signal approaches infinity [31].
The pth-order Shmaliy polynomial is defined as follows: where is the length of the data, is the Pochhammer symbol with , , and the hypergeometric function is given by
In practice, the weighted discrete Shmaliy polynomials [21] are used to minimize the instability at high-order moments. The orthogonal basis is built as follows: where is the weight function and defines the norm.
The weighted discrete Shmaliy polynomials satisfy the following orthogonality condition:
According to Shakibaei Asli and Flusser, it is possible to define multivariate polynomials (2D and 3D) as products of the univariate polynomials (1D) in each dimension.
Therefore, using (3), the discrete Shmaliy moments of order of the image are defined as where and are the size of the image on - and -axes, respectively. 1D and 2D discrete Shmaliy polynomials are shown in Figures 1(a) and 1(b), respectively.
(a)
(b)
2.1. Other Discrete Orthogonal Moments
In order to conduct a comparative study, we also use two well-known discrete orthogonal moments to compare with. A brief introduction of the basic theory of the discrete Tchebichef moments [5] and the discrete Krawtchouk moments [32] is presented in this section.(a)The discrete Tchebichef polynomials are defined in terms of hypergeometric functions as follows: where . The property of orthogonality of the discrete Tchebichef polynomials is satisfied by where is the Kronecker delta and the normalizing factor is defined as For a practical and stable implementation, the normalized discrete Tchebichef polynomials (see Figure 2(a)) are used to build the following basis: Using (9), the discrete Tchebichef moments (DTMs) of the image are computed as where and are the image size on - and -axes, respectively.(b)Discrete Krawtchouk Polynomials [32] constitute an orthogonal polynomial basis related to the binomial distribution. The classic nth-order Krawtchouk polynomial is defined as where and . For stable computation, the weighted discrete Krawtchouk polynomials determine the basis as follows: where is the weight function and is the norm. The weighted discrete Krawtchouk polynomials satisfy the orthogonality condition. From (12), the discrete Krawtchouk moments (DKMs) of the image are defined as (see Figure 2(b)) where and are the size of the image on - and -axes, respectively. DKMs are able to extract local features when the parameter varies between to .
(a)
(b)
2.2. Computational Aspects of the Discrete Orthogonal Moments
We have mentioned that a key issue of the orthogonal polynomials is the numerical instability when the order of the polynomial becomes large. However, another important problem is the computational cost that tends to be high and limits the use of the orthogonal moments when an on-line computation is required [33].
There are many different ways to compute the coefficients of the orthogonal polynomials, i.e., by the hypergeometric definition [34], in terms of the Rodrigues formula [35], using the closed form representation, or by iterating the recurrence equations [36].
Morales-Mendoza et al. [20] showed that the orthogonal Shmaliy polynomials satisfy the three-term recurrence relation. This relation can be used for the numerical implementation of the polynomials. The three terms are the length of the signal, the polynomial order, and the signal index. As we mentioned in Section 2, the discrete Shmaliy polynomials have two main advantages in comparison to Tchebichef and Krawtchouk polynomials: one-parameter recursion definition and linear weight function.
We conducted time complexity experiments to compare the recursive implementation of Shmaliy, Tchebichef, and Krawtchouk polynomials. Specifically, the implementation of the weighted discrete Shmaliy polynomials is based on the n-recursion proposal by Asli and Flusser [21]; the normalized discrete Tchebichef polynomials are computed with the x-recursion formula proposed by Mukundan [5, 17]; and the weighted discrete Krawtchouk polynomials are calculated using the n-recursion definition proposed by Yap et al. [18].
Figure 3 shows the mean time computation of the discrete polynomials up to the 100th order; each polynomial basis was calculated 1000 times. The computation of the Shmaliy polynomials is faster than Tchebichef and Krawtchouk polynomials. Specifically, in comparison to Krawtchouk, Shmaliy polynomials have the same type of recursion but Krawtchouk polynomials depend on two parameters that may affect the computation time. Tchebichef polynomials use different recursion implementation; therefore, it is not adequate to make a direct comparison.
Asli and Flusser [21] performed a time comparison test between Tchebichef and Shmaliy moments. They also compared the hypergeometric and n- recursion implementations and showed that the recursive method is faster than the hypergeometric definition. They concluded that there is no significant differences between the time complexity of Shmaliy and Tchebichef moments. The recursive computation of DOMs and their optimization is a state-of-the-art topic; i.e., Shu et al. [33] studied the computational and time complexity of the recursive computation of Tchebichef moments. They also compared and contrasted their proposal with Kotoulas and Andreadis [37], Papakostas et al. [38], and Wang and Wang [39] recursive methods.
3. Texture Description Based on Discrete Shmaliy Moments
Texture classification involves two main assignments: feature extraction and similarity measure. This paper is focused on feature extraction; however, the target is the correct classification of the image database described in Section 4. A complete scheme of classification is explained below and is based on statistical descriptors through local texture analysis using overlapping square windows.
3.1. Textural Features
The discrete orthogonal moments, described previously in Section 2, are scalar quantities that characterize a specific function. For example, in accordance with (5), the value of the DSMs is determined by the correlation between the image and the polynomial basis , which is calculated using the weighted discrete Shmaliy polynomials.
The polynomial basis acts as a filter for the DOMs. The value of the moment is larger when the variations of the image are similar to in X and Y directions. This characteristic is important for the analysis of the textures because the texture is defined as the spatial repetition of gray-level patterns in a region [40]. Therefore, it is possible to obtain a description of the texture when the magnitude dependency of the sth-order moment is evaluated. The description is related to the frequency responses of the polynomial basis. Thus, the textural features based on the discrete Shmaliy moments are computed as where and is defined in (5).
Marcos and Cristóbal [29] performed a validation process with DTMs as texture descriptors. They compared their proposal against standard methods such as Haralick cooccurrence matrix [41], Gaussian filtering [42], and local binary patterns [43]. The results showed that DTMs collect essential texture information with similar results compared to conventional methods.
The textural features could be computed using any of the polynomial basis mentioned on Section 2, and they are used in the experiments of Section 5.
3.2. Statistical Textural Features
In [29], the authors calculated one textural feature vector per image that implies the computation of high-order moments. According to Mukundan [17], this way of computing moments may produce overflow calculations due to unstable oscillations of the polynomials. The numerical instability issue may appear in Krawtchouk and Tchebichef polynomial bases and is one of the most important reasons why DOMs have not been applied widely on texture description tasks.
Our proposal includes a modification based on overlapping square windows that shift over the image. This modification avoids the computation of high-order moments and assures numerical stability. Since overlapping windows generate an over-description of the image, we calculate the corresponding first three central moments (see Figure 4 that shows the overlapping window description).
The statistical textural feature, , based on the textural features, , is calculated for each window where is the window position, usually left-to-right and top-to-down.
This description scheme is useful in spite of the differences of size of inter- and intra-classes. The vector is built as follows: where , , and are the mean, standard deviation, and kurtosis, respectively.
We can say that (16) describes the texture. The classification scheme shows how the statistical textural features are computed for each image (Figure 5). Once the features are calculated, they are used to determine the classification rule.
3.3. Linear Discriminant Analysis
The computation of the statistical textural features involves large vectors because the textural features (Section 3.1) have length of bins, where is the window size. Therefore, an overfitting may occur, which consists on an excessive adjustment of the classification algorithm leading to inaccurate results when new data are processed [44, 45].
In order to avoid that issue, Fisher proposed the discriminant analysis [46] for binary problems. This method is one of the most used for dimensionality reduction. Linear discriminant analysis (LDA) projects feature vectors into clusters in such a way that preserves the data geometry; LDA is defined as follows: are the samples of classes; then the Fisher’s projection is given by vector that maximizes
The between-scatter matrix is where and is the general average.
On the other hand, the within class scatter matrix is defined by where .
Finally, the optimal projection matrix is more relevant eigenvalues of .
4. IICBU-2008 Benchmark
IICBU-2008 benchmark was presented by Shamir [22]. It is composed of 11 datasets, which represent a wide range of biological imaging problems. The image format of the datasets is TIFF (Tagged Image File Format) of 8 or 16 bits, and the image sizes vary from to pixels. The datasets utilized on this study are shown in Table 1. Next, we present a short review of the datasets used in this study.
4.1. Binucleate
The dataset was acquired for classification purposes of binucleate cellular phenotype and normal cells. Binucleate phenotype is associated with failures in cell division and is a target when specialists look for genes that have effects on cell division. Many agents used on chemotherapy also alter cell division. The cells included in this dataset belong to Drosophila melanogaster and are shown in Figure 6.
(a)
(b)
4.2. Celegans
The full name of the dataset is Caenorhabditis elegans muscle aging and its principal purpose is to estimate the age of nematodes using images of muscle tissue. The images were captured with a light microscope at different ages. The images were assigned to the respective chronological age, rather than the morphological age. These ages are not fully correlated because the physical development among individuals of the same species is not the same even though they may have the same chronological age.
Celegans has a different purpose than the rest of the datasets included in the benchmark IICBU-2008. The principal goal is to find a variable that correlates morphological and chronological ages. The ages vary in the range of 1, 2, 4, and 8 days. A sample of each class is shown in Figure 7.
(a)
(b)
(c)
(d)
4.3. CHO and HeLa
The Chinese Hamster Ovary (CHO) dataset was published in 1998 [47]. CHO is focused on the identification of subcellular organelles and cellular structures. Figure 8 shows one slice of each CHO class. HeLa was published in 2001 [48] as an extension of CHO dataset and supplies new structures for identification task. Figure 9 contains one slice of each HeLa class. These images were originally dedicated to find the location of proteins in cells. The location of one protein in relation to other proteins or functional structures points out the protein purpose.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
4.4. Liver Aging
The Atlas of Gene Expression in Mouse Aging Project (AGEMAP) is a study performed by the National Institute of Aging and the National Institute of Health in the United States. This study involves 48 mice, males and females, with ad libitum (with no restriction) and caloric restriction diet and different ages: 1, 6, 12, and 24 months. The livers were extracted from sacrificed mice; they were stained with Hematoxylin/Eosin (H+E). The images were acquired through brightfield microscopy technique. Figure 10 shows slices of different liver ages. The AGEMAP images could be analyzed on different manners: age, gender, and diet.
(a)
(b)
(c)
(d)
4.5. Lymphoma
The hematologic diseases include different types of leukemia and lymphomas. Although these diseases are not usual, they have aroused medical interest in recent years. For that reason, this dataset is one of the most popular databases because it has been used on several studies for classification [49, 50] and registration [51]. The whole image dataset is stained with Hematoxylin/Eosin (H+E). Lymphoma database is composed of three types of hematologic diseases: chronic lymphocytic leukemia (CLL) with 133 samples, follicular lymphoma (FL) with 139, and mantle cell lymphoma (MCL) with 122. Figure 11 shows slices of the Lymphoma dataset.
(a)
(b)
(c)
5. Experimental Results
During the preprocessing stage, color images in Lymphoma and Liver aging datasets were transformed into gray scale and the size of the images was reduced by 50%.
In addition, the backgrounds of the images in C. elegans, CHO, and HeLa (Figures 7, 8, and 9, respectively) have irrelevant information. For this reason, the background content was discarded through the analysis of the overlapping windows in the following manner: if at least 50% of the pixels within the window contain information, it means the window is not background and the statistical textural features are computed. Otherwise, the window is dismissed.
For local analysis, every image is divided into overlapping square windows. We varied the size of the windows from to pixels, the step of window size is pixels, and the overlap is 50%. This process is shown in Figure 4. For Krawtchouk moments, ; see (14).
The overlapping windows are processed as Figure 5 depicts. The textural features, see (15), are calculated on every overlapping window. The proposed statistical textural features are computed according to their general definition given in (16).
Linear discriminant analysis, described in Section 3.3, is applied to the statistical textural features before the classification and validation process. Afterwards, the feature vectors are mean-centered (zero mean and unit standard the error). Keep in mind that the resulting vectors from LDA are length, where is the number of classes. The new feature vectors are shown in Figure 12. Note that CHO and HeLa datasets are partially shown because their dimensionality exceeds the three-dimensional representation with 4 and 9 features, respectively.
The classification of the datasets is performed through 5-Fold cross-validation. We use five classifiers: k-nearest neighbors (KNN) [52], support vector machines (SVM) with linear and Gaussian kernels [53, 54], naïve Bayes [55], and random forest (RF) [56].
The mean accuracy results for each dataset are shown in the following figures: Binucleate results in Figure 13, Celegans in Figure 14, CHO in Figure 15, HeLa in Figure 16, Liver aging in Figure 17, and Lymphoma in Figure 18.
Every subfigure in Figures 13–18 shows the results for a specific window size. Red bars correspond to Shmaliy moments, blue bars to Tchebichef moments, and green bars to Krawtchouk moments. The vertical lines that appear at the center of the bars represent the standard deviations computed by averaging the cross-validated results.
5.1. ANOVA Analysis
N-way ANOVA [57], or factorial ANOVA, is performed in order to find out if the choice of the window size, classifier, and DOM is statistically significant because the classification results are very similar to each other.
ANOVA compares the means between groups and resolves if any of them has significant statistical difference from each other. ANOVA tests the null hypothesis, , where is the mean of group . In other words, the null hypothesis is that any group election does not have relevance. If this hypothesis is rejected, the alternative hypothesis, , is accepted and implies that at least two group means are significantly different from each other.
In our experiments, the accuracy is the dependent variable. The window size, classifier, and discrete orthogonal moment are the three nominal or independent variables; it means . N-way ANOVA can also examine the interaction between the independent variables because the interactions show that the differences between the groups are not uniform in the categories of the independent variables. The result after applying ANOVA is the p-values; one p-value is computed for each nominal variable and the interactions among them. In general, if the p-value is less than 5%, the null hypothesis is rejected. Table 2 shows the p-values, in percentages, of independent variables and their interactions.
Based on Table 2, N-way ANOVA shows that the independent variables are statistically significant, separately, for all image databases.
First, the window size is important because the larger the size, the better the overall accuracy classification results. The classifier is also statistically relevant; it shows that better results are obtained, in general, with both naïve Bayes and SVM (linear and Gaussian) classifiers. DOMs have also a statistical significance. The best classification results for Celegans and HeLa were obtained with Krawtchouk moments. For CHO the best results were obtained with Tchebichef moments, and for Liver aging and Lymphoma discrete Shmaliy moments scored the best results. This situation may occur because the datasets are more correlated to a certain polynomial basis.
It is important to mention that Binucleate dataset was not included in ANOVA because, for all the parameter combinations, the classification results were very high.
The best results from our experiments are compared with the state-of-the-art works: Shamir et al. [22, 58], Siji et al. [59], and Meng et al. [60] classified the same datasets from IICBU-2008 benchmark.
Shamir et al. presented IICBU-2008 benchmark in [22]. They also used a popular tool called WND-CHRM [61, 62] to extract and classify different feature sets, for example, Radon transform, Gabor coefficients, multiscale histograms, and Tamura texture features. Afterwards, Shamir et al. used Fourier, discrete Tchebichef moments, and Wavelet transforms up to fourth order to extract texture features; WND-CHRM was used again to perform the classification. They reported results with a subset of data, including CHO and HeLa datasets.
Siji et al. proposed to apply principal component analysis (PCA) and Fisher analysis score based on feature selection for image retrieval applications. They used the same features computed by WND-CHRM for Binucleate, Celegans, CHO, and HeLa databases. They used naïve Bayes and SVMs as classifiers.
Meng et al. proposed a classification model called collateral representative subspace projection modeling (C-RSPM) for Liver Aging and Lymphoma datasets.
Table 3 summarizes and compares the results of our proposal and the ones reported in the state-of-the-art.
6. Discussion and Conclusions
Since their recent introduction, Shmaliy moments have not been treated or tested on image texture analysis. In this paper, we show for the first time the capability of Shmaliy moments for computing texture features and compare them with the most popular families of discrete orthogonal moments: Tchebichef and Krawtchouk moments. Moreover, the results are compared with a set of methods that use the IICBU-2008 benchmark.
The proposed classification method is sufficiently general and useful in many different types of images. We included a preprocessing stage that consists in color space transformation and, in some cases, resizing the datasets to accelerate the description and classification processes. Our proposal is robust and it may be useful in many different cases that not necessarily belong to the biomedical field.
The classification results obtained by discrete Shmaliy moments are similar to the ones achieved with Tchebichef and Krawtchouk moments; therefore, a 3-way ANOVA was applied in order to verify if the parameter selection (descriptor, window size, and moment) was significant. ANOVA suggests that the descriptor combination is the most important parameter because the null hypothesis was rejected for all databases while the discrete orthogonal moment and the window size elections show equal low significance since the null hypothesis was rejected only on three of five databases.
On the other hand, comparison among discrete orthogonal moments and other techniques found in the literature shows competitive results. Accordingly, Shmaliy moments have shown the same capability for texture description as other classical discrete polynomial bases; therefore, its consideration is pertinent for texture analysis, and probably a combination of Shmaliy texture features and other discrete moments could be the topic for a future work.
Data Availability
To download the IICBU-2008 benchmark, use https://ome.grc.nia.nih.gov/iicbu2008/.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work has been sponsored by UNAM Grants PAPIIT IN116917 and SECITI 110/2015. Germán González thanks CONACYT263921 Scholarship.