# A Comparative Study on Discrete Shmaliy Moments and Their Texture-Based Applications

Authors: Germán González; Rodrigo Nava; Boris Escalante-Ramírez
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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 through5-fold cross-validation, and the results are computed using random forest, support vector machines, naïve Bayes, and k-nearest neighbors classifiers.


---

## Body

## 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 Section2, 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 f(x,y) on an orthogonal polynomial basis rpq(x,y). Formally, they are defined as follows: Mpq=∫x∫yrpq(x,y)f(x,y)dxdy. 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 N 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].Thepth-order Shmaliy polynomial is defined as follows: (1)hpn,N=-1pp+1n-ppN-npp!Np+1×F23-p,n+1,1-N+n;n-p,1-N-p+n;1where N is the length of the data, (x)p=Γ(x+p) is the Pochhammer symbol with Γ(x)=(x-1)!, p,n={0,1,…,N-1∣N>0}, and the hypergeometric function Fxy is given by(2)Fxya1,…,ax;b1,…,by;z=∑k=0∞a1ka2k⋯axkb1kb2k⋯bykzkk!.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: (3)hp^n,N=hpn,Nwn,Nρp,Nwhere w(n,N)=2n/N(N-1) is the weight function and ρ(p,N)=(p+1)(N-p-1)p/N(N)p+1 defines the norm.The weighted discrete Shmaliy polynomials satisfy the following orthogonality condition:(4)∑n=0N-1hp^n,Nhq^n,N=δpq.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 Hnm of order p<N of the image f(x,y) are defined as (5)Hnm=∑x=1N∑y=1Mh^nx,N+1h^my,M+1fx,ywhere N and M are the size of the image on X- and Y-axes, respectively. 1D and 2D discrete Shmaliy polynomials are shown in Figures 1(a) and 1(b), respectively.Figure 1
According to Shakibaei Asli and Flusser [21], the 2D discrete Shmaliy polynomials are obtained as the product of two unidimensional polynomials. (a) 1D and (b) 2D weighted discrete Shmaliy polynomials from order 0 to 4.
(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:(6)tnx,N=1-Nn3F2-n,-x,1+n;1,1-N;1wherex,n={0,1,…,N-1}.The property of orthogonality of the discrete Tchebichef polynomials is satisfied by(7)∑x=0N-1tnxtmx=ρn,Nδnmwhereδnm is the Kronecker delta and the normalizing factor ρ(n,N) is defined as (8)ρn,N=NN2-12N2-22⋯N2-n22n+1=2n!N+n2n+1.For a practical and stable implementation, the normalized discrete Tchebichef polynomials (see Figure2(a)) are used to build the following basis: (9)t^nx,N=tnx,Nρn,N.Using (9), the discrete Tchebichef moments (DTMs) of the image f(x,y) are computed as (10)τnm=∑x=0N-1∑y=0M-1t^nxt^myfx,ywhereN and M are the image size on X- and Y-axes, respectively.(b)
Discrete Krawtchouk Polynomials [32] constitute an orthogonal polynomial basis related to the binomial distribution. The classicnth-order Krawtchouk polynomial is defined as (11)Knx;p,N=∑k=0Nak,n,pxk=F12-n,-x;-N;1pwherex,n={0,1,2,…,N∣N>0} and pε(0,1).For stable computation, the weighted discrete Krawtchouk polynomials determine the basis as follows:(12)K^nx;p,N=Knx;p,Nwx;p,Nρn;p,Nwherew(x;p,N)=Nxpx(1-p)N-x is the weight function and ρ(n;p,N)=(-1)n1-p/pnn!/(-N)n is the norm.The weighted discrete Krawtchouk polynomials satisfy the orthogonality condition.(13)∑x=0NK^nx;p,NK^mx;p,N=δnmFrom (12), the discrete Krawtchouk moments (DKMs) of the image f(x,y) are defined as (see Figure 2(b)) (14)Qnm=∑x=0N-1∑y=0MK^nx;p1,N-1K^my;p2,M-1fx,ywhereN and M are the size of the image on X- and Y-axes, respectively.DKMs are able to extract local features when the parameterp varies between 0 to 1.Figure 2
Unidimensional discrete orthogonal polynomials. (a) Tchebichef and (b) Krawtchouk polynomial bases from order 0 to 4.
(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 then-recursion proposal by Asli and Flusser [21]; the normalized discrete Tchebichef polynomials are computed with thex-recursion formula proposed by Mukundan [5, 17]; and the weighted discrete Krawtchouk polynomials are calculated using then-recursion definition proposed by Yap et al. [18].Figure3 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.Figure 3
Mean computation time in milliseconds of the discrete Shmaliy, Tchebichef, and Krawtchouk polynomials.Asli and Flusser [21] performed a time comparison test between Tchebichef and Shmaliy moments. They also compared the hypergeometric andn- 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.


## 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:(6)tnx,N=1-Nn3F2-n,-x,1+n;1,1-N;1wherex,n={0,1,…,N-1}.The property of orthogonality of the discrete Tchebichef polynomials is satisfied by(7)∑x=0N-1tnxtmx=ρn,Nδnmwhereδnm is the Kronecker delta and the normalizing factor ρ(n,N) is defined as (8)ρn,N=NN2-12N2-22⋯N2-n22n+1=2n!N+n2n+1.For a practical and stable implementation, the normalized discrete Tchebichef polynomials (see Figure2(a)) are used to build the following basis: (9)t^nx,N=tnx,Nρn,N.Using (9), the discrete Tchebichef moments (DTMs) of the image f(x,y) are computed as (10)τnm=∑x=0N-1∑y=0M-1t^nxt^myfx,ywhereN and M are the image size on X- and Y-axes, respectively.(b)
Discrete Krawtchouk Polynomials [32] constitute an orthogonal polynomial basis related to the binomial distribution. The classicnth-order Krawtchouk polynomial is defined as (11)Knx;p,N=∑k=0Nak,n,pxk=F12-n,-x;-N;1pwherex,n={0,1,2,…,N∣N>0} and pε(0,1).For stable computation, the weighted discrete Krawtchouk polynomials determine the basis as follows:(12)K^nx;p,N=Knx;p,Nwx;p,Nρn;p,Nwherew(x;p,N)=Nxpx(1-p)N-x is the weight function and ρ(n;p,N)=(-1)n1-p/pnn!/(-N)n is the norm.The weighted discrete Krawtchouk polynomials satisfy the orthogonality condition.(13)∑x=0NK^nx;p,NK^mx;p,N=δnmFrom (12), the discrete Krawtchouk moments (DKMs) of the image f(x,y) are defined as (see Figure 2(b)) (14)Qnm=∑x=0N-1∑y=0MK^nx;p1,N-1K^my;p2,M-1fx,ywhereN and M are the size of the image on X- and Y-axes, respectively.DKMs are able to extract local features when the parameterp varies between 0 to 1.Figure 2
Unidimensional discrete orthogonal polynomials. (a) Tchebichef and (b) Krawtchouk polynomial bases from order 0 to 4.
(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 then-recursion proposal by Asli and Flusser [21]; the normalized discrete Tchebichef polynomials are computed with thex-recursion formula proposed by Mukundan [5, 17]; and the weighted discrete Krawtchouk polynomials are calculated using then-recursion definition proposed by Yap et al. [18].Figure3 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.Figure 3
Mean computation time in milliseconds of the discrete Shmaliy, Tchebichef, and Krawtchouk polynomials.Asli and Flusser [21] performed a time comparison test between Tchebichef and Shmaliy moments. They also compared the hypergeometric andn- 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 Section4. 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 Section2, 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 f(x,y) and the polynomial basis rnm(x,y), which is calculated using the weighted discrete Shmaliy polynomials.The polynomial basisrnm(x,y) acts as a filter for the DOMs. The value of the moment is larger when the variations of the image are similar to rnm(x,y) 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 thesth-order moment is evaluated. The description is related to the frequency responses of the polynomial basis. Thus, the textural features M(s) based on the discrete Shmaliy moments are computed as (15)Ms=∑s=m+nHmnwhere s=0,1,…,2N-2 and Hmn 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 Section2, 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 Figure4 that shows the overlapping window description).Figure 4
The description of the textures is based on sliding windows that move top-to-down and left-to-right. On every window position, the corresponding orthogonal moments are computed.The statistical textural feature,t, based on the textural features, Mi(s), is calculated for each window where i 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 vectort is built as follows: (16)t=μMi0,σMi0,κMi0,…,μMi2N-2,σMi2N-2,κMi2N-2where μ, σ, 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.Figure 5
General scheme of our proposal to classify texture images. After computing local DOMs using sliding windows, the statistical textural feature vectors are computed and classified.


### 3.3. Linear Discriminant Analysis

The computation of the statistical textural features involves large vectors because the textural features (Section3.1) have length of 2N-1 bins, where N 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: X1={x11,x21,…,xl11},…,XC={x1C,x2C,…,xlCC} are the samples of C classes; then the Fisher’s projection is given by vector w that maximizes (17)Jw=wTSBwwTSWw.The between-scatter matrix is(18)SB=∑i=1Cliμi-μμi-μT,where μi=1/li∑j=1lixji and μ is the general average.On the other hand, the within class scatter matrix is defined by(19)SW=∑i=1CSi,where Si=∑j=1li(xji-μi)(xji-μi)T.Finally, the optimal projection matrixw∗ is C-1 more relevant eigenvalues of SW-1SB.


## 3.1. Textural Features

The discrete orthogonal moments, described previously in Section2, 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 f(x,y) and the polynomial basis rnm(x,y), which is calculated using the weighted discrete Shmaliy polynomials.The polynomial basisrnm(x,y) acts as a filter for the DOMs. The value of the moment is larger when the variations of the image are similar to rnm(x,y) 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 thesth-order moment is evaluated. The description is related to the frequency responses of the polynomial basis. Thus, the textural features M(s) based on the discrete Shmaliy moments are computed as (15)Ms=∑s=m+nHmnwhere s=0,1,…,2N-2 and Hmn 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 Section2, 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 Figure4 that shows the overlapping window description).Figure 4
The description of the textures is based on sliding windows that move top-to-down and left-to-right. On every window position, the corresponding orthogonal moments are computed.The statistical textural feature,t, based on the textural features, Mi(s), is calculated for each window where i 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 vectort is built as follows: (16)t=μMi0,σMi0,κMi0,…,μMi2N-2,σMi2N-2,κMi2N-2where μ, σ, 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.Figure 5
General scheme of our proposal to classify texture images. After computing local DOMs using sliding windows, the statistical textural feature vectors are computed and classified.


## 3.3. Linear Discriminant Analysis

The computation of the statistical textural features involves large vectors because the textural features (Section3.1) have length of 2N-1 bins, where N 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: X1={x11,x21,…,xl11},…,XC={x1C,x2C,…,xlCC} are the samples of C classes; then the Fisher’s projection is given by vector w that maximizes (17)Jw=wTSBwwTSWw.The between-scatter matrix is(18)SB=∑i=1Cliμi-μμi-μT,where μi=1/li∑j=1lixji and μ is the general average.On the other hand, the within class scatter matrix is defined by(19)SW=∑i=1CSi,where Si=∑j=1li(xji-μi)(xji-μi)T.Finally, the optimal projection matrixw∗ is C-1 more relevant eigenvalues of SW-1SB.


## 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 25×25 to 1388×1040 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.Table 1
Characteristics of the datasets included in the IICBU-2008 benchmark.
Database Classes Num. of img. Format Microscopy Binucleate 2 40 1280 × 1024 16bit Fluorescence Celegans 4 252 1600 × 1200 16bit Fluorescence Liver aging 4 850 1388 × 1040 32bit color Brightfield CHO 5 340 512 × 382 16bit Fluorescence HeLa 10 860 382 × 382 16bit Fluorescence Lymphoma 3 375 1388 × 1040 32bit color Brightfield


### 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 toDrosophila melanogaster and are shown in Figure 6.Figure 6
The two classes defined in theBinucleate dataset are (a) binucleate and (b) nonbinucleate.
(a) (b)


### 4.2. Celegans

The full name of the dataset isCaenorhabditis 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.Figure 7
The chronological ages of theCelegans dataset: (a) 1, (b) 2, (c) 4, and (d) 8 days.
(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 eachCHO class.HeLa was published in 2001 [48] as an extension ofCHO dataset and supplies new structures for identification task. Figure 9 contains one slice of eachHeLa 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.Figure 8
Classes of theCHO dataset: (a) Anti-giantin, (b) Hoechst 33258 (DNA), (c) Anti-lamp2, (d) Anti-nop4, and (e) Anti-tubulin.
(a) (b) (c) (d) (e)Figure 9
Classes ofHeLa dataset: (a) Actin, (b) DNA (Nuclei), (c) TfR (Endosome), (d) ER (Endoplasmic reticulum), (e) Giantin (cis/medial Golgi), (f) GPP130 (cis Golgi), (g) Lamp2 (Lysosomes), (h) Tubulin, (i) Mitochondria, and (j) Nucleolin (Nucleoli).
(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, withad 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.Figure 10
The images belong toad libitum diet from withLiver Aging database: (a) 1, (b) 6, (c) 16, and (d) 24 months.
(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 theLymphoma dataset.Figure 11
Classes ofLymphoma database: (a) CLL, (b) FL, and (c) MCL.
(a) (b) (c)


## 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 toDrosophila melanogaster and are shown in Figure 6.Figure 6
The two classes defined in theBinucleate dataset are (a) binucleate and (b) nonbinucleate.
(a) (b)


## 4.2. Celegans

The full name of the dataset isCaenorhabditis 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.Figure 7
The chronological ages of theCelegans dataset: (a) 1, (b) 2, (c) 4, and (d) 8 days.
(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 eachCHO class.HeLa was published in 2001 [48] as an extension ofCHO dataset and supplies new structures for identification task. Figure 9 contains one slice of eachHeLa 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.Figure 8
Classes of theCHO dataset: (a) Anti-giantin, (b) Hoechst 33258 (DNA), (c) Anti-lamp2, (d) Anti-nop4, and (e) Anti-tubulin.
(a) (b) (c) (d) (e)Figure 9
Classes ofHeLa dataset: (a) Actin, (b) DNA (Nuclei), (c) TfR (Endosome), (d) ER (Endoplasmic reticulum), (e) Giantin (cis/medial Golgi), (f) GPP130 (cis Golgi), (g) Lamp2 (Lysosomes), (h) Tubulin, (i) Mitochondria, and (j) Nucleolin (Nucleoli).
(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, withad 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.Figure 10
The images belong toad libitum diet from withLiver Aging database: (a) 1, (b) 6, (c) 16, and (d) 24 months.
(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 theLymphoma dataset.Figure 11
Classes ofLymphoma database: (a) CLL, (b) FL, and (c) MCL.
(a) (b) (c)


## 5. Experimental Results

During the preprocessing stage, color images inLymphoma andLiver aging datasets were transformed into gray scale and the size of the images was reduced by 50%.In addition, the backgrounds of the images inC. elegans,CHO, andHeLa (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 from20×20 to 50×50 pixels, the step of window size is 10×10 pixels, and the overlap is 50%. This process is shown in Figure 4. For Krawtchouk moments, p1=p2=0.5; see (14).The overlapping windows are processed as Figure5 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 Section3.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 C-1 length, where C is the number of classes. The new feature vectors are shown in Figure 12. Note thatCHO andHeLa datasets are partially shown because their dimensionality exceeds the three-dimensional representation with 4 and 9 features, respectively.Figure 12
The statistical textural features are transformed with linear discriminant analysis. Note that the projections obtained improve the quality of the clusters and the boundaries among them.The classification of the datasets is performed through5-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, andLymphoma in Figure 18.Figure 13
Mean accuracy results forBinucleate dataset. This dataset consists of two classes. The dataset has the particularity that LDA projects and separates the classes in such a way that the classification becomes trivial. However, when SVM-rbf is used, there is a small error caused by the influence of the parameters of the Gaussian kernel.Figure 14
Mean accuracy results forCelegans dataset. The dataset consists of four classes. The best results were obtained when the size of the local window is bigger than 20.Figure 15
Mean accuracy results forCHO dataset. The dataset consists of five classes of cellular structures of hamster ovary. On this dataset, Shmaliy and Tchebichef features show a similar and uniform performance when the window size is bigger than 30.Figure 16
Mean accuracy results forHeLa dataset. The dataset, as an extension ofCHO database, consists of ten classes of cellular structures of hamster ovary. The best results were obtained when the size of the local window is bigger than 40.Figure 17
Mean accuracy results forLiver Aging dataset. This dataset consists of four classes. The dataset has the particularity that classifiers perform similarly among them, and the accuracy increases with larger windows.Figure 18
Mean accuracy results forLymphoma dataset. The dataset consists of three classes of hematological diseases. The best results were obtained when the size of the local window is bigger than 40.Every subfigure in Figures13–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,H0=μ1=μ2=⋯=μk, where μ is the mean of group k. In other words, the null hypothesis is that any group election does not have relevance. If this hypothesis is rejected, the alternative hypothesis, H1, 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 meansN=3.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 thep-values; onep-value is computed for each nominal variable and the interactions among them. In general, if thep-value is less than 5%, the null hypothesis H0 is rejected. Table 2 shows thep-values, in percentages, of independent variables and their interactions.Table 2
p-values of N-way ANOVA.
WSize & WSize & Class. Database WSize Class. Mom. Class Mom. Mom. Celegans 0 0.4 0 97.3 0 99.4 CHO 0 0.5 3.2 98.6 2.7 98.7 HeLa 0 0 0 97.5 39.7 79.9 Liver aging 0 0 0.6 99.7 0 98.2 Lymphoma 0 0 0 97.3 3.3 99Based on Table2,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 forCelegans andHeLa were obtained with Krawtchouk moments. ForCHO the best results were obtained with Tchebichef moments, and forLiver aging andLymphoma 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 thatBinucleate 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, includingCHO andHeLa 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 forBinucleate,Celegans,CHO, andHeLa 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) forLiver Aging andLymphoma datasets.Table3 summarizes and compares the results of our proposal and the ones reported in the state-of-the-art.Table 3
Comparison results between discrete orthogonal moments and state-of-the-art methods.
Shmaliy Tchebichef Krawtchouk WND- Shamir(2) Siji Meng Database Moments Moments Moments CHRM et al. et al. et al. Binucleate 100 100 100 100 - 100 - Celegans 100 100 100 53 - 75 - Liver aging 99.81 100 99.05 51 - - 95.09 CHO 100 100 100 93 94 99 - HeLa 98.25 97.79 98.15 84 83 74 - Lymphoma 100 99.74 99.74 85 - - 84.79


## 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,H0=μ1=μ2=⋯=μk, where μ is the mean of group k. In other words, the null hypothesis is that any group election does not have relevance. If this hypothesis is rejected, the alternative hypothesis, H1, 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 meansN=3.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 thep-values; onep-value is computed for each nominal variable and the interactions among them. In general, if thep-value is less than 5%, the null hypothesis H0 is rejected. Table 2 shows thep-values, in percentages, of independent variables and their interactions.Table 2
p-values of N-way ANOVA.
WSize & WSize & Class. Database WSize Class. Mom. Class Mom. Mom. Celegans 0 0.4 0 97.3 0 99.4 CHO 0 0.5 3.2 98.6 2.7 98.7 HeLa 0 0 0 97.5 39.7 79.9 Liver aging 0 0 0.6 99.7 0 98.2 Lymphoma 0 0 0 97.3 3.3 99Based on Table2,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 forCelegans andHeLa were obtained with Krawtchouk moments. ForCHO the best results were obtained with Tchebichef moments, and forLiver aging andLymphoma 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 thatBinucleate 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, includingCHO andHeLa 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 forBinucleate,Celegans,CHO, andHeLa 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) forLiver Aging andLymphoma datasets.Table3 summarizes and compares the results of our proposal and the ones reported in the state-of-the-art.Table 3
Comparison results between discrete orthogonal moments and state-of-the-art methods.
Shmaliy Tchebichef Krawtchouk WND- Shamir(2) Siji Meng Database Moments Moments Moments CHRM et al. et al. et al. Binucleate 100 100 100 100 - 100 - Celegans 100 100 100 53 - 75 - Liver aging 99.81 100 99.05 51 - - 95.09 CHO 100 100 100 93 94 99 - HeLa 98.25 97.79 98.15 84 83 74 - Lymphoma 100 99.74 99.74 85 - - 84.79


## 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.



---

*Source: 1673283-2018-08-09.xml*